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.
