Equivalence of Poisson and Weibull PH regression in a survival setting

Question

Suppose we have the following dataset that records individual survival times (dur) and a covariate z:

id| dur | z
-------------
1 |  1  | -1
2 |  2  |  1
3 |  3  | -1
4 |  4  |  1

I want to model the duration as a function of z. I may specify dur ~ weibull() and parameterize the scale as a function of z. This model can be fitted easily, e.g. with phreg or survreg in R.

If I am now interested to incorporate time-dependent covariates, I need to transform the dataset into something like this:

id | event  | time  | z  | orig.row
_____________________________________
1  | 1      | 1     | -1 | 1
2  | 0      | 1     | 1  | 2
2.1| 1      | 2     | 1  | 2
3  | 0      | 1     |-1  | 3
3.1| 0      | 2     |-1  | 3
3.2| 1      | 3     |-1  | 3
4  | 0      | 1     | 1  | 4
4.1| 0      | 2     | 1  | 4
4.2| 0      | 3     | 1  | 4
4.3| 1      | 4     | 1  | 4

Apparently, a to say event ~ poisson() with a "poisson mean" $$\log(u_i) = \beta_0 + \alpha \log(t_i) + \beta_1 z_i$$ is equivalent to a the Weibull model above (Lindsey 1995). But when I run the analysis in R, I get two largely different values for the scale:

# Generate the data: 
library(eha)
enter <- rep(0, 4)
exit <- 1:4
event <- rep(1, 4)
z <- rep(c(-1, 1), 2)
dat <- data.frame(enter, exit, event, z)
binDat <- toBinary(dat)
binDat <- binDat[order(rownames(binDat)),]

# Run the model: 
summary(phreg(Surv(enter,exit, event) ~ z, data = dat, dist="weibull"))
summary(glm(event ~ z + log(risktime), data = binDat, family = poisson("log")))

The shape parmeter (glm: log(risktime)) is sort of similar as well as the beta for z. But the scale parameter is different (glm: Intercept). What am I doing wrong?

## Results
summary(phreg(Surv(enter,exit, event) ~ z, data = dat, dist="weibull"))
Call:
phreg(formula = Surv(enter, exit, event) ~ z, data = dat, dist = "weibull")

Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
z                   0.200    -0.432     0.649     0.519     0.406 

log(scale)                    1.020     2.774     0.197     0.000 
log(shape)                    0.985     2.677     0.421     0.019 

Events                    4 
Total time at risk            10 
Max. log. likelihood      -5.6567 
LR test statistic         0.68 
Degrees of freedom        1 
Overall p-value           0.410367

summary(glm(event ~ z + log(risktime), data = binDat, family = poisson("log")))

Call:
glm(formula = event ~ z + log(risktime), family = poisson("log"), 
    data = binDat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0696  -0.7695  -0.5391   0.3482   1.0671  

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)    -1.6122     1.0019  -1.609    0.108
z              -0.3166     0.5149  -0.615    0.539
log(risktime)   1.0634     1.0928   0.973    0.330

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 7.3303  on 9  degrees of freedom
Residual deviance: 6.1388  on 7  degrees of freedom
AIC: 20.139

Number of Fisher Scoring iterations: 5

Answer 1

Poisson regression is not equivalent to a a Weibull survival model. Instead, it's assuming an exponential distribution, where the baseline hazard is not only proportional, but constant.

The Weibull model relaxes this assumption somewhat. While the two models will give you the same answer if indeed the underlying survival distribution is exponential, as the exponential distribution is a special case of the weibull distribution, they need not under other circumstances.

It appears one such circumstance is your data.

Answer 2

In the usual setup for Poisson regression, the "risktime" variable would be offered as an offset.

glm(event ~ z +offset(log(risktime)), family ="poisson")
     selected output
Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
(Intercept)  -1.5623     0.5000  -3.125  0.00178 **
z            -0.3095     0.5000  -0.619  0.53589   

The Weibull with a shape parameter of 1 and the glm(... , family="poisson") call and the CPH models with (truly) constant hazard should give similar results, since they are all basically the same. I am not familiar with pkg:eha, but here is a warning from Therneau in his survreg.distributions page: "The location-scale parameterizaion of a Weibull distribution found in survreg is not the same as the parameterization of rweibull." And the comment in the example code: "# survreg parameters are scale=1/shape, intercept=log(scale)". On the other hand this post to r-help says pkg:eha uses "standard R parameterization"

Edit:

On viewing this question again, I see that I did not grasp the essential question. Yes. the log(time) covariate should create a similar model to to the Weibull. A value of 0 for the estimated coefficient would correspond to the exponential model with constant effect over time while values above 0 would yield models with increasing risk over time. In R regression packages there are generally predict functions. That does not seem to be available in G:oran Brostr:om's eha package. He does include a vignette with the package which has further information on how he uses the parameters. See page 4 of http://cran.r-project.org/web/packages/eha/vignettes/parametric.pdf to see that he uses the log of the parameters that are typical.

Then looking at the output we see that the log(scale) parameter (which is what would correspond to the time-associated increase in risk) in the phreg estimate is 1.02 while the glm( ~log(time), "poisson") estimate is 1.06. I would say those are basically the same.