[R] Variance-covariance matrix
Peter Dalgaard
P.Dalgaard at biostat.ku.dk
Wed Aug 27 12:09:06 CEST 2008
Laura Bonnett wrote:
> Here is the exact code I have written which does the standard vs nt1
> and standard vs nt2 and also gives me the hazard ratio for nt1 vs nt2.
>
> with <- read.table("allwiths.txt",
> header=TRUE)
> fix(arm)
> function (data)
> {
> dummy <- rep(0,2437)
> for(i in 1:2437){
> if(data$Arm[i]=="CBZ")
> dummy[i] <- i
> }
> return(data[-dummy,])
> }
> fix(x1)
> function (data)
> {
> x1 <- rep(0,716)
> for(i in 1:716){
> if(data$Treat[i]=="TPM")
> x1[i] <- 0
> if(data$Treat[i]=="VPS")
> x1[i] <- 0
> if(data$Treat[i]=="LTG")
> x1[i] <- 1
> }
> return(x1)
> }
> fix(x2)
> function (data)
> {
> x2 <- rep(0,716)
> for(i in 1:716){
> if(data$Treat[i]=="TPM")
> x2[i] <- 1
> if(data$Treat[i]=="VPS")
> x2[i] <- 0
> if(data$Treat[i]=="LTG")
> x2[i] <- 0
> }
> return(x2)
> }
> fit1 = coxph(Surv(Withtime,Wcens)~x1(armb),data=armb,method="breslow")
> exp(fit1$coefficients)
> exp(confint(fit1))
> fit2 = coxph(Surv(Withtime,Wcens)~x2(armb),data=armb,method="breslow")
> exp(fit2$coefficients)
> exp(confint(fit2))
> exp(fit2$coefficients)/exp(fit1$coefficients)
>
> From that, how do I get the necessary variance-covaraince matrix.
>
> Sorry if I appear dense. It really isn't my intention.
>
You need to read up on vectorized arithmetic, (and don't use fix() to
define functions, x1 <- function(....) {....} will do, and what's with
the unused arm() function??).
However, if I get your data structure correctly, what I meant was that
you could just look at
fit3 <- coxph(Surv(Withtime,Wcens)~Treat, data=armb,method="breslow")
(possibly factor(Treat) instead). If VPS is your control, I think you'll
even get the desired coefficient directly from the output because the
group comparisons will be against LTG.
--
O__ ---- Peter Dalgaard Øster Farimagsgade 5, Entr.B
c/ /'_ --- Dept. of Biostatistics PO Box 2099, 1014 Cph. K
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
More information about the R-help
mailing list