Dear Göran I too have been working on this for the last couple of weeks and this is a great help thanks. I wonder can you advise on how the covariance matrix of estimated parameters in phreg relates to those estimated in survreg.
Im hoping you can just put your finger on this. Your eha package has been really useful so far and if I can get this working then i will be a very happy person indeed. Essentially, I am doing all this in order to run phreg, then convert and pass the results to a survreg object in order to use the predict function to create survival curves for each stratification with confidence intervals around them. This is clunky I know, but unless there is a similar predict method for phreg objects I am stuck, and I havent found one yet. What follow is a piece of code ive been working on to explore the relationships between phreg and survreg outputs using simulated data. # ------------------------------------------------------------------------------ # START OF SCRIPT EXAMPLE # ------------------------------------------------------------------------------ rm(list=ls()) # clear the memory graphics.off() # turn graphics off set.seed(1) # for repeatability, set the random number seed # Load the libraries library(survival) library(eha) # ------------------------------------------------------------------------------ n <- 50 # number of observations per stratification true.shape <- c(3,3,3,3) # common shape parameter for each category true.scale <- c(3,5,7,20) # separate scale parameter X1.code <- c("A","A","B","B") # X1 encoded as A or B X2.code <- c("P","Q","P","Q") # X2 encoded as P or Q # Create the response variables and two fixed factor variables Y <- NULL X1 <- NULL X2 <- NULL # Loop over each combinatin of Factors and simulate some data from a # weibull distribution for (i in 1:length(true.shape)){ # simulate the Y data tmp <- rweibull(n,true.shape[i],scale=true.scale[i]) Y <- c(Y,tmp) # Encode the fixed factors X1 <-c (X1,rep(X1.code[i],n)) X2 <-c (X2,rep(X2.code[i],n)) } # close loop over factors # X1 and X1 are factors X1 <- as.factor(X1) X2 <- as.factor(X2) # Everyone has status 1 here (i.e. they were all observed to death/failure) status <- rep(1,length(Y)) # ------------------------------------------------------------------------------ # plot the data dev.new() plot(survfit(Surv(Y,status)~X1+X2),"All Data") # run survreg m.surv <- survreg(Surv(Y,status)~X1*X2,dist="weibull") # run phreg m.ph <- phreg(Surv(Y,status)~X1*X2,dist="weibull") # Covariance matrix for survreg() method m.surv$var # Covariance matrix for phreg() method m.ph$var # ------------------------------------------------------------------------------ # And these are the type of plots i want to produce. new.data <- data.frame(X1=factor("B",levels=levels(X1)), X2=factor("Q",levels=levels(X2))) prob <- seq(0.0, 1.0, length=200) temp <- predict(m.surv, type='quantile', p=prob, se.fit=T, newdata=new.data) matlines(cbind(temp$fit, temp$fit +1.96*temp$se, temp$fit - 1.96*temp$se), 1-prob, lty=c(1,2,2), col=2, lwd=1) # ------------------------------------------------------------------------------ # END OF SCRIPT EXAMPLE # ------------------------------------------------------------------------------ As far as I can decipher the following are related: var(log(shape)) from phreg equates directly to var(log(scale)) from survreg var(log(scale)) = var(Intercept) -cov(log(shape),log(scale) = cov(Intercept,log(scale)) # note the minus in front of the phreg output I am stuck though now on how the following relate The variances of the covariates var(X1B) var(X2Q) var(X1B:X2Q) The covariance of the covariates cov(X1B,X2Q) cov(X1B:X2Q,X1B) cov(X1B:X2Q,X2Q) The covariance of the covariates with the shape and scale parameters cov(log(shape),X1B) cov(log(scale),X1B) .. and similary with X2Q and X1B:X2Q I look forward to hearing from you on this. best wishes Andrew Jackson -- Dr Andrew Jackson Lecturer School of Natural Sciences Zoology Building, Trinity College Dublin, Dublin 2, Ireland Tel. + 353 1 896 2728, Fax. + 353 1 677 8094, Email. a.jack...@tcd.ie http://www.tcd.ie/Zoology/research/research/theoretical/AndrewJackson.php -- View this message in context: http://r.789695.n4.nabble.com/survreg-vs-aftreg-eha-the-relationship-between-fitted-coefficients-tp3082263p3219204.html Sent from the R help mailing list archive at Nabble.com. [[alternative HTML version deleted]]
______________________________________________ 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.