<br><br><div>Dear ProbABEL developers,<br><br>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?<br>
<br>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.<br>
<br>I downloaded the source code with the command:<br>$ svn checkout svn://<a href="http://svn.r-forge.r-project.org/svnroot/genabel/">svn.r-forge.r-project.org/svnroot/genabel/</a><br><br>From the genabel directory, I created the patch (attached)<br>
$ svn diff > ../svn.diff.addcoxscore.patch<br><br>#################<br># ProbABEL command<br>~/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<br>
<br>#################<br># ProbABEL output - see attached "coxrobust.out.22Aug.dev.txt"<br><br><div>##################<br>## R output<br><br>> resid[1:50]<br><br> [1] 0.0000000000 0.1086258272 0.0011085883 -0.0014120667 0.0008283735<br>
[6] 0.0008186706 0.1066213091 -0.0019426329 0.1084221333 -0.0041228497<br>[11] -0.0052619695 0.0032972320 -0.8810385183 0.0031559494 0.0034423290<br>[16] 0.0031134893 0.0024421635 0.1001538285 -0.8829088458 0.0046289440<br>
[21] 0.0041215655 0.1022440869 -0.8828842583 0.0068676096 -0.0083150569<br>[26] 0.0077873828 0.0090103153 0.0074532482 0.0950587064 0.0070292693<br>[31] 0.0114207938 0.0922110596 -0.0063786832 0.0124870186 0.0093464958<br>
[36] -0.0095624990 0.0927836541 0.0142498550 0.0119589572 -0.0124419529<br>[41] 0.0156961998 0.0066048746 0.0097690163 -0.0126152895 0.0089864928<br>[46] 0.0129145714 -0.0122369770 0.0076375518 0.0121946660 0.0806112482<br>
<br>########################<br>## R code to be called from the genabel/branches/ProbABEL-pacox/v.0.1-3/examples directory<br><br>in.pheno <- read.table("coxph_data.txt",header=T,<a href="http://as.is">as.is</a>=T)<br>
in.pheno[1:5,]<br>attach(in.pheno)<br><br>dose <- read.table("./test.mldose",header=F,<a href="http://as.is">as.is</a>=T)<br>dose[1:5,]<br>table(in.pheno[,"id"]== sapply(dose[,1],function(x){strsplit(x,split="->")[[1]][2]}))<br>
<br># Look at first SNP<br><br>i <- 3<br><br>library(survival)<br><br>test.coxph <- coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + dose[,i],x=TRUE,y=TRUE )<br><br>n <- length(test.coxph$residuals)<br># rr <- object$residuals<br>
y <- test.coxph$y<br><div>x <- test.coxph[['x']] # avoid matching object$xlevels<br>vv <- test.coxph$var<br><br>weights <- rep(1,n)<br><br>ny <- ncol(y)<br>status <- y[,ny,drop=TRUE]<br>nvar <- ncol(x)<br>
<br>ord <- order(y[,ny-1], -status)<br>newstrat <- rep(0,n)<br>newstrat[n] <- 1<br><br># sort the data<br>x <- x[ord,]<br>y <- y[ord,]<br>score <- exp(test.coxph$linear.predictors)[ord]<br><br>method<-c("erfron")<br>
<br>resid <- .C("coxscore", as.integer(n),<br> as.integer(nvar),<br> as.double(y),<br> x=as.double(x),<br> as.integer(newstrat),<br> as.double(score),<br>
as.double(weights[ord]),<br> as.integer(method=='efron'),<br> resid= double(n*nvar),<br> double(2*nvar))$resid<br><br>tmp.rr <- matrix(0, n, nvar)<br><br>tmp.rr <- matrix(resid, ncol=nvar)<br>
<br>rr <- matrix(0, n, nvar)<br>rr[ord,] <- matrix(resid, ncol=nvar)<br>dimnames(rr) <- list(names(test.coxph$residuals), <br><br>names(test.coxph$coefficients))<br> rr <- rr %*% vv<br> rr <- rr * weights ## At this point, rr is what is returned by a call to residuals(test.coxph,type="dfbeta")<br>
<br><br>robust.variances <- t(rr) %*% rr<br>test.coxph.robust$var<br><br><br>## compare to call to get robust variances<br>test.coxph.robust <- coxph(Surv(fupt_chd,chd) ~ sex + age + othercov + dose[,i],x=TRUE,y=TRUE,robust=T )<br>
<br><br><br><br>--<br><br>Alisa Knodle Manning, Ph.D.<br>Postdoctoral Research Fellow<br>Broad Institute<br>617-714-7662<br><a href="mailto:amanning@broadinstitute.org">amanning@broadinstitute.org</a>
</div></div></div>