Why is Logistic Regression Wrong?
If you had observed all subjects for the same amount of time, logistic regression would be fine. The censoring -- not knowing what happens after 180 days -- is the problem.
Think of this in an extreme case. Suppose we're interested in examining the risk of death. However, our sample consists mostly of children (who do not typically die young, and often go on to live into adulthood). The short observation time of the children means they will likely contribute an excess of 0's to our estimate of the risk of death. Hence, that estimate is inappropriate to apply to the entire population; the children simply have not enough time to accrue risk of death by ... much of anything. Adults, by contrast, have had more opportunities to die from crossing roads, smoking cigarettes, having plaque accumulate in their arteries, and so on and so on. Observing all subjects for the same amount of time means they have had (assuming all else equal) the same time to accumulate the same risk.
How Can We Analyze Such Data?
You could use a cox model, but you could also do a Poisson regression with robust covariance estimation and use exposure time as an offset, as is explained here. The {survival}
R library has a function called pyears
which calculates the length of exposure. The docs for pyears
show an example of how to use the exposure time as an offset, but they do not show how to estimate robust covariance in that example.
Let's see how to do this with an example. First, let's simulate some data. Assume we've randomly assigned patients to two groups. Group membership is indicated by trt
.
N <- 1000
trt <- rbinom(N, 1, 0.5)
raw_tm <- rexp(N, 1/(100 - 50*trt))
cens_tm <- 365
tm <- ceiling(pmin(raw_tm, cens_tm))
event <- 1.0*(raw_tm < cens_tm)
d <- data.frame(tm, event, trt)
Now, let's use pyears
to count the outcomes in each group and the total person years contributed.
library(survival)
p <- pyears(Surv(tm/365, event)~trt, data=d, scale=1, data.frame=T)
p$data
> p$data
trt pyears n event
1 0 145.32877 519 503
2 1 67.09041 481 481
There were 503 events in the 145 person years contributed by the trtr=0
group. That's a rate of 503/145.33 ~ 3.46
events per person year. This means if I followed 10 people for a year I would expect 10 people x 1 year x 3.46 = 34.6
events. Same if I had followed 5 people for 2 years, or 20 people for 6 months.
Now, onto the regression
fit <- glm(event ~ trt + offset(log(pyears)),
data=p$data,
family = poisson)
coeftest(fit, vcov. = sandwich)
z test of coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2416e+00 2.5167e-32 4.9334e+31 < 2.2e-16 ***
trt 7.2823e-01 8.2724e-16 8.8031e+14 < 2.2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The base rate of events in the trt=0
per 1 person year is exp(1.24) ~ 3.46
-- which is exactly what we saw. We can also see that the rate of events per 1 person year in the trt=1
group is exp(0.73) = 2.07
times as large as compared to the trt=0
group. That also makes sense from the pyear
tabulation we saw above.
So far as the 180 day cliff is concerned, your data are censored at 180 days. That is still informative, because we know that those individuals do not get the outcome within 180 days. Their person years and 0 events still count towards the regression.