[GenABEL-dev] Issues with implementing robust standard errors for cox proportional hazard regression

[Alisa Manning]

Dear ProbABEL developers,

I've been working on implemening robust standard errors for the cox
proportional hazards analysis. The problem that I've run into is that the
function call that I've implemented (copied from the R survival package
much like the coxfit2 function) returns different values from what I
expect.  I've implemented the call in R with what I think is the same input
data.  Is anyone available to spend a few minutes looking over this code
and let me know if you can spot my mistake?

I modified the branched version of 0.1-3-- If I should implement these
changes in another version, let me know.  I chose this version because the
Cox Proportional Hazard analysis was working straight off; I had plans to
try implementing the changes (and working out the cox bugs) in the more
recent versions AFTER I understood the necessary changes in 0.1-3.

I downloaded the source code with the command:
$  svn checkout svn://svn.r-forge.r-project.org/svnroot/genabel/

>From the genabel directory, I created the patch (attached)
$ svn diff > ../svn.diff.addcoxscore.patch

#################
# ProbABEL command
~/amanning/ProbABEL/genabel/branches/ProbABEL-pacox/v.0.1-3/examples 531 $
 ../bin/pacoxph -p ./coxph_data.txt -d ./test.mldose -i ./test.mlinfo -m
./test.map -c 19 --robust -o coxph.robust >
../../../../../coxrobust.out.22Aug.dev.txt

#################
# ProbABEL output - see attached "coxrobust.out.22Aug.dev.txt"

##################
## R output

> resid[1:50]

 [1]  0.0000000000  0.1086258272  0.0011085883 -0.0014120667  0.0008283735
 [6]  0.0008186706  0.1066213091 -0.0019426329  0.1084221333 -0.0041228497
[11] -0.0052619695  0.0032972320 -0.8810385183  0.0031559494  0.0034423290
[16]  0.0031134893  0.0024421635  0.1001538285 -0.8829088458  0.0046289440
[21]  0.0041215655  0.1022440869 -0.8828842583  0.0068676096 -0.0083150569
[26]  0.0077873828  0.0090103153  0.0074532482  0.0950587064  0.0070292693
[31]  0.0114207938  0.0922110596 -0.0063786832  0.0124870186  0.0093464958
[36] -0.0095624990  0.0927836541  0.0142498550  0.0119589572 -0.0124419529
[41]  0.0156961998  0.0066048746  0.0097690163 -0.0126152895  0.0089864928
[46]  0.0129145714 -0.0122369770  0.0076375518  0.0121946660  0.0806112482

########################
## R code to be called from the
genabel/branches/ProbABEL-pacox/v.0.1-3/examples directory

in.pheno <- read.table("coxph_data.txt",header=T,as.is=T)
in.pheno[1:5,]
attach(in.pheno)

dose <- read.table("./test.mldose",header=F,as.is=T)
dose[1:5,]
table(in.pheno[,"id"]==
sapply(dose[,1],function(x){strsplit(x,split="->")[[1]][2]}))

# Look at first SNP

i <- 3

library(survival)

test.coxph <- coxph(Surv(fupt_chd,chd) ~ sex + age + othercov +
dose[,i],x=TRUE,y=TRUE )

n <- length(test.coxph$residuals)
#  rr <- object$residuals
y <- test.coxph$y
x <- test.coxph[['x']]  # avoid matching object$xlevels
vv <- test.coxph$var

weights <- rep(1,n)

ny <- ncol(y)
status <- y[,ny,drop=TRUE]
nvar <- ncol(x)

ord <- order(y[,ny-1], -status)
newstrat <- rep(0,n)
newstrat[n] <- 1

# sort the data
x <- x[ord,]
y <- y[ord,]
score <- exp(test.coxph$linear.predictors)[ord]

method<-c("erfron")

resid <- .C("coxscore", as.integer(n),
            as.integer(nvar),
            as.double(y),
            x=as.double(x),
            as.integer(newstrat),
            as.double(score),
            as.double(weights[ord]),
            as.integer(method=='efron'),
            resid= double(n*nvar),
            double(2*nvar))$resid

tmp.rr <-  matrix(0, n, nvar)

tmp.rr <- matrix(resid, ncol=nvar)

rr <- matrix(0, n, nvar)
rr[ord,] <- matrix(resid, ncol=nvar)
dimnames(rr) <- list(names(test.coxph$residuals),

names(test.coxph$coefficients))
 rr <- rr %*% vv
 rr <- rr * weights   ## At this point, rr is what is returned by a call to
residuals(test.coxph,type="dfbeta")


robust.variances <- t(rr) %*% rr
test.coxph.robust$var


## compare to call to get robust variances
test.coxph.robust <- coxph(Surv(fupt_chd,chd) ~ sex + age + othercov +
dose[,i],x=TRUE,y=TRUE,robust=T )




--

Alisa Knodle Manning, Ph.D.
Postdoctoral Research Fellow
Broad Institute
617-714-7662
amanning at broadinstitute.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20120822/42406f62/attachment-0001.html>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: coxrobust.out.22Aug.dev.txt
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20120822/42406f62/attachment-0001.txt>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: svn.diff.addcoxscore.patch
Type: application/octet-stream
Size: 17239 bytes
Desc: not available
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20120822/42406f62/attachment-0001.obj>