[GenABEL-dev] Issues with implementing robust standard errors for cox proportional hazard regression
L.C. Karssen
lennart at karssen.org
Sun Sep 30 16:56:38 CEST 2012
Hi Alisa,
I've created a branch as per your request:
svn co svn checkout
svn://svn.r-forge.r-project.org/svnroot/genabel/branches/ProbABEL-pacox-robust
I've applied your patch and added the R-script you used to test the code
as well (r.964).
Best,
Lennart.
On 29/09/12 18:52, Alisa Manning wrote:
> 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
>> ------------------------------------------------------------------
>>
>
>
>
--
-----------------------------------------------------------------
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/20120930/11046e47/attachment.sig>
More information about the genabel-devel
mailing list