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 ~~~~~~~~~~ - ([EMAIL PROTECTED]) FAX: (+45) 35327907 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.