I am in doubt whether it is possible to get BLUP in this model. As far as I know, it is not possible to get any unbiased estimator of Odds ratios in a fixed effect model. And since the predictions of the random effects should be adjusted for the estimated Odds ratios of the fixed effects it should also not be possible to give unbiased predictions of random effects. In practice, however, maksimum likelihood estimates, as GLIMMIX calculates, are good if just the sample size are not extreme small. ... View more

Here you have the program that calculate Aalen Johansen estimates. It was not written as a macro, so bit editing will perhaps be neccessary. I hope you find it useful. ... View more

You can not calculate all the survival curves in only one run of phreg. You can of course run phreg for each outcome. The Fine and Gray method which proc phreg use only take account for one outcome at a time. I am not sure that if you run phreg for each outcome that the total sum will add up to one (exactly).  If you like, I can give you a macro which calculate all curves in one time using the Aalen Johansen estimator. ... View more

Of course you can use Cox regression also when all individuals have an event. As in this simulated example where all individuals have an observed event. The cox regression estimate the rate-ratio fairly close to the true hazard ratio. data simulation; do i=1 to 1000; exposure=(i<=500); t=rand('exponential',exp(-0.5*exposure)); output; end; run; proc phreg data=simulation; model t=exposure; run;    Estimation problems will only occur if all individuals in one level of exposure are censored, or if all indivudals from one level come have events before any of the individuals form other levels. ... View more

I think you did it correct.   The difference must be that the cox-regression allows the baseline rate to depend on the underlying time. In contrast, Poissonregression have a constant baseline hazard, which you estimate with the intercept.   If you are very eager to see if this the difference is due to time, then you can divide your observations into pieces, and make a Poissonregression that include a timedependent variable, "time", that says what timeinterval that are observed. Remember that only the last interval should have can have an event.     ... View more

I agree that Cox regression and Poisson regression usually gives almost same result. That is because the difference between the two models are very small. Its likely that you did something wrong. ... View more

I see no other explanation than the Hazard Ratios is calculated with "6" as a reference value. Therefore the HR at 6 is 1, and the confidence intervals at 6 is just a point [1-1] (or undefined). ... View more

In case your age variable is time dependent (age increase as time increase) then you have to divide each waiting time into smaller pieces, where the first intervals will end up with a censoring (count=0) and only the last interval will end up with a event (or censoring if the original waiting time was censored). The offset should be log(length of interval). If you just want to have age at beginning of the interval, you just put that into the model as a covariate, without splitting up the interval.   On the left side in the model statement you have the count which is either 0 or 1.   If you want to speed up the calculation time, you can summarize your data and have sum-of-events on left side and use log(sum-of-length of-waiting time) as offset. ... View more

Hi Lei, It is described very detailed in my favourite book: "Survival and Event History of Analysis" of Odd Aalen, page 223.   But the quick argument is like this:   Let T be exponential distributed with rate exp(β).  Make the likelihood function for β. Then lets say X is Poisson distributed with intensity t*exp(β). Make the likelihod function for β. Observe that these two likelihood functions are the same if you put log(t) as offset in the latter. ... View more

You will not get biased result because of different numbers of control-persons. But of course, as more control person the more statistical power you will have. Just remember to put the variable that defines the matching group into the strata statement in proc logistic (it could forexample be the case-ID-number). ... View more

Waiting time can be modelled with poisson regression, but not with negative binomial regression. This is becuase the likelihood for exponential distributed data (or even just data with piecewice constant rate)  is the same as the likelihood from poisson distributed data.   In this example I model an event which has a constant value of 1 with poisson regression, to estimate the rate which true value is 1/5. data mydata; event=1; do i=1 to 1000; t=rand('exponential',5); logt=log(t); output; end; run; proc genmod data=mydata; model event=/dist=poisson link=log offset=logt; estimate 'rate' intercept 1/exp; run; But it is not meaningfull to talk about overdispersion, because that is a term which use assumption about the distribution. The data is clearly not poisson-distributed, but as the likelihood is the same as from poisson distrubuted data, poisson regression can be used to estimate parameters.   ... View more

You can not separately get an odds ratio of the variables that are involved in some interactionsterms. But, you can easily get the oddsratios for a variable for each level of the other variable it interacts with: just use the "at" option in the odds ratio statement. Like in this example which calculate the oddsratio for the variable A at each level of B, and the oddsratio for B for at level of A.   data mydata; n=1000; do a=0 to 1; do b=0 to 1; odds=(1.2**a)*(0.75**b)*(0.8**(a*b)); p=odds/(1+odds); y=rand('binomial',p,n); output; end; end; run; proc logistic data=mydata; class a(ref="0") b(ref="0")/param=glm ; model y/n=a b a*b; oddsratio a/at(b=all) diff=ref; oddsratio b/at(a=all) diff=ref; run;   ... View more

Proc phreg does not calculate the expected lifetime directly. But PHREG can calculate the survival function, which then can be used to calculate the expected lifetime. It is such that the integrated survival function gives the expected lifetime. One should be carefull in practice, since the survival function can be difficult to estimate in the tail.   Here is an example, where the datastep after PHREG do the integration: data mydata; do i=1 to 10000; predictor=mod(i,2); time=rand('gamma',5*exp(log(2)*predictor)); censurtime=rand('gamma',10); event=(time<=censurtime); time=min(time,censurtime); output; end; keep time predictor event; run; data covariates; do predictor=0 to 1; output; end; run; proc phreg data=mydata; class predictor /param=glm; model time*event(0)=predictor; baseline out=kaplanmeier survival=survival covariates=covariates/method=pl; run; data meansurvival; set kaplanmeier; by predictor; retain lasttime; if first.predictor then integral=0; else do; integral+survival*(time-lasttime); end; lasttime=time; if last.predictor; put integral; keep integral predictor; run;   Unfortunately, it is not possible to estimate confidence limits of the expected survival in this way. If you want this, you should better use a failuretime model (proc lifereg). ... View more