[GenABEL-dev] Issues with implementing robust standard errors for cox proportional hazard regression
Alisa Manning
amanning at broadinstitute.org
Sat Sep 29 18:52:58 CEST 2012
Hey--- Thanks for looking in to this..... it's the fact that the
output is the same up to row 6 that made me bang my head on my desk.
I'm still planning on trying to tackle this problem. Could you
create a new code branch for me to work with?
Thanks,
Alisa
On Sat, Sep 29, 2012 at 12:16 PM, L.C. Karssen <lennart at karssen.org> wrote:
> 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
> ------------------------------------------------------------------
>
--
Alisa Knodle Manning, Ph.D.
Postdoctoral Research Fellow
Broad Institute
617-714-7662
amanning at broadinstitute.org
More information about the genabel-devel
mailing list