[GenABEL-dev] Issues with implementing robust standard errors for cox proportional hazard regression
L.C. Karssen
lennart at karssen.org
Sat Sep 29 18:16:18 CEST 2012
Dear Alisa,
I finally found some time to take a look at the coxscore problem. Your
patch applies well on the latest v0.1.3 branch, however, I haven't been
able to tackle the problem yet. As far as I can see the variables you
send into the coxscore() function in pacoxph are OK.
The strange thing is that the output of coxscore (in the tempresiduals
matrix) is (virtually) identical to the output in R, but only for the
first six rows of the matrix. On row seven things are different all of a
sudden.
I replaced the -O3 option with -O0 in the Makefile to make sure no
optimisations were done, but no luck.
I also used valgrind to check for memory leaks. There are a few of them,
but they are either in coxfit2() and coxscore() (so not really our
concern), or in parts of pacoxph that I don't think are relevant for us
now.
Anyway, I take another look at it today or tomorrow. If someone else
wants to help us out: you're most welcome!
Best,
Lennart.
On 22/08/12 21:47, Alisa Manning wrote:
>
>
> 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/
> <http://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 <http://as.is>=T)
> in.pheno[1:5,]
> attach(in.pheno)
>
> dose <- read.table("./test.mldose",header=F,as.is <http://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 <mailto:amanning at broadinstitute.org>
>
>
> _______________________________________________
> genabel-devel mailing list
> genabel-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/genabel-devel
>
--
-----------------------------------------------------------------
L.C. Karssen
Utrecht
The Netherlands
lennart at karssen.org
http://blog.karssen.org
Stuur mij aub geen Word of Powerpoint bestanden!
Zie http://www.gnu.org/philosophy/no-word-attachments.nl.html
------------------------------------------------------------------
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 198 bytes
Desc: OpenPGP digital signature
URL: <http://lists.r-forge.r-project.org/pipermail/genabel-devel/attachments/20120929/97d4398c/attachment.sig>
More information about the genabel-devel
mailing list