Dear All,

Please can anyone help me to solve this simple problem that I cannot solve. I tried to read several books but cannot figure it out. I added the STATA code at the end of this post

The problem is the following: I am trying to simulate some survival data from a weibull model with a treatment effect that decreases with time. The "gamma" specification of 1 in survsim is in fact the exponential but I put this just for simplicity:

Code:
survsim pfs_time pfs_event, distribution(weibull) lambdas(.00440942) gammas(1) covariates(trt -.80) tde(trt 0.20) maxtime(365)
I then wanted to see the results of the weibull model with a time varying coefficient so I stsplit the data and then added to the model the variable "gen lntimextrt=ln(_t)*trt" so I could confirm the estimate I inputed in survsim "tde(trt 0.20)"

Code:
stset pfs_time, fail(pfs_event) id(personid)
stsplit time, every(1)
keep personid trt _d _t _t0
rename _t pfs_time
rename _d pfs_event
rename _t0 origintime
stset pfs_time, fail(pfs_event) id(personid) enter(origintime) origin(origintime)
gen lntimextrt=ln(_t)*trt
streg i.trt c.lntimextrt, distribution(weibull)
All of the above seems to be correct

However, I would like to manually compare the kaplan meier plot by treatment arm from simulated data with the weibull model parameters I used in survsim. However, I cannot figure this out.

Basically, I would like to derive the survival function using the following parameters I used in survsim:

"lambdas(.00440942) gammas(1) covariates(trt -.80) tde(trt 0.20)"

This should be a simple thing to do but I cannot figure out how to do it. I tried to understand how to do this by deriving manually "predict surv, surv" and "predict csurv, csurv"

I know that "predict csurv, csurv" is just the running product of "predict surv, surv" when there are multiple records for the same subject. However, I do not know how to obtain "predict surv, surv" and therefore I cannot figure out how to derive the theoretical survival weibull curve under the assumed estimates above.

Please could anyone help me with this?

Please see the code below to understand what I mean

Code:
drop _all
set seed 83474
set obs 2660 //####
gen personid=_n
gen trt = rbinomial(1,0.5)
//#########################
survsim pfs_time pfs_event, distribution(weibull) lambdas(.00440942) gammas(1) covariates(trt -.80) tde(trt 0.20)  maxtime(365) //####
//#########################
stset pfs_time, fail(pfs_event) id(personid)

stsplit time, every(1)

keep personid trt _d _t _t0

rename _t pfs_time
rename _d pfs_event
rename _t0 origintime

stset pfs_time, fail(pfs_event) id(personid) enter(origintime) origin(origintime)

gen lntimextrt=ln(_t)*trt

streg i.trt c.lntimextrt, distribution(weibull)

// csurv is the running product within each subject of surv
predict double surv_wei, surv

bysort personid: generate double runprod = surv_wei[1]
by personid: replace runprod = runprod[_n-1]*surv_wei if _n>1

predict double csurv_wei, csurv

list if runprod!=csurv_wei

gen double H=exp(_b[_cons]+_b[1.trt]*trt+_b[lntimextrt]*lntimextrt)*(_t^exp(_b[/ln_p]))

gen double csurv_fromH=exp(-H)

// it seems that csurv_fromH matches csurv_wei when trt==0 but not when trt==1? why???
tab trt if csurv_fromH!=csurv_wei

// I would like to obtain this but comparing the observed simulated data with the theorical model
// with assumed survsim paramaters. HOW DO I DO THIS?
sts graph, by(trt) addplot((line csurv_wei _t if trt==0, sort) (line csurv_wei _t if trt==1, sort))
Many thanks in advance for your help!

Andrew