06-11-2015 10:25 AM
I am doing simulation to compare a couple of methods. I simulate survival data from the one-level parametric frailty model with Weibull baseline hazard, normally distributed clusters. I generate 500 of these datasets at the same time in the dataset "Survival_MultiLevel". Each of these 500 datasets is recognizable by "sample_id". Clusters are, also, determined by "cluster_id". I have two covariates as well.
Since I used to get errors when using the codes posted on SAS website to model parameteric frailty models in NLMIXED, I wrote the loglikelihood in this way. These initial values also help no to get over-flow values. "qpoints = 30" also helps to get more precise estimated. I used to qpoints = 5, qpoints = 10 and ... and I used to get worst estimates. In a word, many "NOTE"s and "ERROR"s led me to use the following codes. These codes work perfectly and I get very appealing results for estimation of b1 and b2. However, estimations of logsig, rho and b0 are very biased.
ods output on;
proc nlmixed data=Survival_MultiLevel qpoints = 30 noad;
bounds rho > 0 ;
parms rho = 2 b1 = 0 b2 = 0 logsig = .5 b0 = 0;
eta = b0 + b1 * covariate1 + b2 * covariate2 + W;
lambda = exp(eta);
log_S_t = -lambda * (survival_t ** rho);
log_h = (censor_flag&j = 1) * (eta + log(rho) + (rho - 1) * log(survival_t));
ll = log_h + log_S_t;
model survival_t ~ general(ll);
random W ~ normal(0, exp(2 * logsig)) subject = cluster_id ;
ods select off;
ods output ParameterEstimates = data&j;