R Textbook Examples
Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence
by Judith D. Singer and John B. Willett
Chapter 12: Extending the Discrete-Time Hazard Model

The examples on the page were written in R version 2.4.1. 


Table 12.2, page 413. To recreate this table, we will run each model separately and then combine the model statistics.

tenure<-read.table("http://www.ats.ucla.edu/stat/examples/alda/tenure_pp.csv", sep=",", header=T)
attach(tenure)
PERIOD2 <- PERIOD^2
PERIOD3 <- PERIOD^3
PERIOD4 <- PERIOD^4
PERIOD5 <- PERIOD^5
PERIOD6 <- PERIOD^6
general <- glm(EVENT ~ D1 + D2 + D3 + D4 + D5 + D6 + D7 + D8 + D9, family = "binomial")
order0<-glm(EVENT~1, family="binomial")
order1<-glm(EVENT~PERIOD, family="binomial")
order2<-glm(EVENT~PERIOD + PERIOD2, family="binomial")
order3<-glm(EVENT~PERIOD + PERIOD2 + PERIOD3, family="binomial")
order4<-glm(EVENT~PERIOD + PERIOD2 + PERIOD3 + PERIOD4, family="binomial")
order5<-glm(EVENT~PERIOD + PERIOD2 + PERIOD3 + PERIOD4 + PERIOD5, family="binomial")
dev <- c(order0$deviance, order1$deviance, order2$deviance, order3$deviance, order4$deviance, order5$deviance, general$deviance)
dev.diff.p <- c(0, dev[1:5] - dev[2:6])
dev.diff.gen <- c(dev - dev[7])
aic <- c(order0$aic, order1$aic, order2$aic, order3$aic, order4$aic, order5$aic, general$aic)
table12.2 <- cbind(dev, dev.diff.p, dev.diff.gen, aic)
table12.2

           dev   dev.diff.p dev.diff.gen       aic
[1,] 1037.5652   0.00000000   206.361395 1039.5652
[2,]  867.4619 170.10332478    36.258070  871.4619
[3,]  836.3041  31.15779577     5.100274  842.3041
[4,]  833.1725   3.13158698     1.968687  841.1725
[5,]  832.7427   0.42981046     1.538877  842.7427
[6,]  832.7322   0.01051811     1.528358  844.7322
[7,]  831.2038   0.00000000     0.000000  849.2038
 

Figure 12.1, page 414.  First, we will rerun the general model without an intercept for greater ease in using the model coefficients later. 


general <- glm(EVENT ~ D1 + D2 + D3 + D4 + D5 + D6 + D7 + D8 + D9 - 1, family = "binomial")

In the loop below, we will use our model estimates to generate a fitted value for each of the nine time points from each of the models.  These are the values that will be plotted in the first of the three plots. We also calculate the survival and hazard rates from the quadratic and general models to be plotted in the second and third graphs.

fits <- c()
survivor.quad = 1
survivor.gen = 1
for (i in 1:9){
  constant = order0$coef[1]
  linear = order1$coef[1] + order1$coef[2]*i
  quadratic = order2$coef[1] + order2$coef[2]*i + order2$coef[3]*i**2
  cubic = order3$coef[1] + order3$coef[2]*i + order3$coef[3]*i**2 + order3$coef[4]*i**3
  hazard.quad = 1/(1 + exp(-quadratic));
  survivor.quad = (1 - hazard.quad)*survivor.quad;
  generalval = general$coef[i]
  hazard.gen = 1/(1 + exp(-generalval));
  survivor.gen = (1 - hazard.gen)*survivor.gen;
  z <- c(i, constant, linear, quadratic, cubic, generalval, hazard.quad, survivor.quad, hazard.gen, survivor.gen)
  fits <- rbind(fits, z)}

Each of the three graphs is created by generating the plot area while plotting one line, and then overlaying any additional lines on the same plot.

par(mfrow=c(1,1))
plot(fits[,1], fits[,2], type = "l", lty = 1, xlim = c(0,9), ylim = c(-6,0), xlab = "Years after hire", ylab = "Fitted logit(hazard)")
points(fits[,1], fits[,3], type = "l", lty = 2)
points(fits[,1], fits[,4], type = "l", lty = 3)
points(fits[,1], fits[,5], type = "l", lty = 4)
points(fits[,1], fits[,6], type = "l", lty = 5)
legend("bottomright", c("constant", "linear", "quadratic", "cubic", "general"), lty = c(1, 2, 3, 4, 5)) 


par(mfrow=c(1,2))
plot(fits[,1], fits[,7], type = "l", lty = 1, xlim = c(0,9), ylim = c(0,.4), xlab = "Years after hire", ylab = "Fitted hazard")
points(fits[,1], fits[,9], type = "l", lty = 2)
legend("bottomright", c("quadratic","general"), lty = c(1, 2)) 


plot(fits[,1], fits[,8], type = "l", lty = 1, xlim = c(0,9), ylim = c(0,1), xlab = "Years after hire", ylab = "Fitted survival")
points(fits[,1], fits[,10], type = "l", lty = 2)
legend("bottomleft", c("quadratic","general"), lty = c(1, 2)) 


Figure 12.3, page 423.  After reading in the data, we will subset it by parenting transition (PT) values.  Then, we will run separate models on these datasets using both logit and cloglog transformations--four models in total.  Finally, we will plot the parameter estimates associated with each of the time points from each of the four models.

firstsex<-read.table("http://www.ats.ucla.edu/stat/examples/alda/firstsex_pp.csv", sep=",", header=T)
firstsex0  <- subset(firstsex, pt==0)
firstsex1  <- subset(firstsex, pt==1)
fs0.logit<-glm(event~d7+d8+d9+d10+d11+d12 - 1, family=binomial(link = "logit"), data = firstsex0)
fs1.logit<-glm(event~d7+d8+d9+d10+d11+d12 - 1, family=binomial(link = "logit"), data = firstsex1)
fs0.loglog<-glm(event~d7+d8+d9+d10+d11+d12 - 1, family=binomial(link = "cloglog"), data = firstsex0)
fs1.loglog<-glm(event~d7+d8+d9+d10+d11+d12 - 1, family=binomial(link = "cloglog"), data = firstsex1)

Having run these four models, we create a single dataset that contains the coefficients from each model and a column of time points. This is the dataset we will use to generate the plot. First, we create a plot with one of the lines. Then, we overlay the other three lines.

fig12.3 <- cbind(time = c(7, 8, 9 , 10, 11, 12), fs0.logit = fs0.logit$coef, fs0.loglog = fs0.loglog$coef, fs1.logit = fs1.logit$coef, fs1.loglog = fs1.loglog$coef)
par(mfrow=c(1,1))
plot(fig12.3[,1], fig12.3[,2], type = "l", ylab = "Transformed hazard probability", xlab = "Grade", ylim = c(-4,0))
points(fig12.3[,1], fig12.3[,3], type = "l", lty=2)
points(fig12.3[,1], fig12.3[,4], type = "l", lty=3)
points(fig12.3[,1], fig12.3[,5], type = "l", lty=4)
legend(7, 0, c("pt=0, logit", "pt=0, loglog", "pt=1, logit", "pt=1, loglog"), lty = c(1, 2, 3, 4)) 

Table 12.3, page 424. For this table, we will run two models with different link functions and then find the fitted baseline hazard using the corresponding model coefficients.


fs.loglog<-glm(event~d7+d8+d9+d10+d11+d12+pt - 1, family=binomial(link = "cloglog"), data = firstsex)
fs.logit<-glm(event~d7+d8+d9+d10+d11+d12+pt - 1, family=binomial(link = "logit"), data = firstsex)

est.logit <- fs.logit$coef
est.loglog <- fs.loglog$coef
fbh.logit <- 1/(1+exp(-est.logit))
fbh.loglog <- 1-exp(-exp(est.loglog))

tab12.3 <- cbind(est.loglog, est.logit, fbh.loglog, fbh.logit)
tab12.3 
    
    est.loglog  est.logit fbh.loglog  fbh.logit
d7  -2.9732973 -2.9943266 0.04984906 0.04768284
d8  -3.6592441 -3.7001239 0.02542322 0.02412410
d9  -2.3156320 -2.2811235 0.09398896 0.09269841
d10 -1.9001387 -1.8225985 0.13890279 0.13912236
d11 -1.7621385 -1.6542272 0.15774911 0.16053845
d12 -1.3426380 -1.1790570 0.22983911 0.23522180
pt   0.7854048  0.8736184 0.88845136 0.70549806

fs.loglog$deviance
[1] 634.5013

fs.logit$deviance
[1] 634.6618

Figure 12.4, page 432. To generate these graphs, we will first do the calculations required for the scatter plots. Then we will run the model required to add the curves.

depression<-read.table("http://www.ats.ucla.edu/stat/examples/alda/depression_pp.csv", sep=",", header=T)

For each age between 4 and 39, we are interested in the proportion of subjects that experienced the event (in this example, an initial depressive episode) within the group of subjects with divorced parents (PD = 1) and the group of subjects without divorced parents (PD = 0).  To calculate this, we can loop through these combinations of age and PD values, subsetting our data and calculating the proportion within each subset.  For the second graph, we will take the logit of this proportion. 

percents <- c()
for (i in 4:39){
  for (j in 0:1){
    x <- subset(depression, PERIOD==i & PD==j)
    y <- table(x$EVENT)
    if (dim(y)==2){
      z <- c(i, j, y[[2]], y[[1]], (y[[2]]/(y[[1]]+y[[2]])))}
    if (dim(y)==1){
      z <- c(i, j, 0, y[[1]], NA)}
    percents <- rbind(percents, z)		
  }
}

percents <- cbind(percents, log(percents[,5]/(1-percents[,5])))
colnames(percents) <- c("age", "parent", "event", "nonevent", "percent", "logpercent")
percents.nm <- as.data.frame(na.omit(percents))
percent.pd0 <- subset(percents.nm, parent == 0)
percent.pd1 <- subset(percents.nm, parent == 1)

Next, we will fit a logistic model to the data.  Then, for each age between 4 and 39 and for each PD group, we are interested in the model-predicted proportion of subjects that experienced the event.  To calculate this, we can again loop through these combinations of age and PD values and use the model coefficients to generate predicted proportions.  We can also calculate the hazard rate for the second plot.

dmodel<-glm(EVENT ~ ONE + age_18 + age_18sq + age_18cub + PD - 1, family=binomial(link = "logit"), data = depression)

modelfit<- c()
for (i in 0:1){
  for (j in 4:39){
    fitx = dmodel$coef[1] + dmodel$coef[2]*(j-18) + dmodel$coef[3]*(j-18)^2 + dmodel$coef[4]*(j-18)^3 + dmodel$coef[5]*i
    fithazard = 1/(1 + exp(-fitx))
    modelfit <- rbind(modelfit, c(i, j, fitx, fithazard))
  }
}

Lastly, we subset our data and generate the plots, starting each with one scatter plot, then layering on the other scatter points and the two fitted curves.

modelfit.0 <- subset(data.frame(modelfit), modelfit[,1]==0)
modelfit.1 <- subset(data.frame(modelfit), modelfit[,1]==1)

par(mfrow=c(1,2))
plot(percent.pd0$age, percent.pd0$percent, pch = 19, ylim = c(0, .06), xlab = "Age", ylab = "Proportion experiencing event")
points(percent.pd1$age, percent.pd1$percent, pch = 22)
points(modelfit.0[,2], modelfit.0[,4], type = 'l', lty = 1)
points(modelfit.1[,2], modelfit.1[,4], type = 'l', lty = 2)
legend(5, 0.06, c("PD = 0", "PD = 1"), lty = c(1, 2))

plot(percent.pd0$age, percent.pd0$logpercent, pch = 19, ylim = c(-8, -2), xlab = "Age", ylab = "Logit(proportion experiencing event)")
points(percent.pd1$age, percent.pd1$logpercent, pch = 22)
points(modelfit.0[,2], modelfit.0[,3], type = 'l', lty = 1)
points(modelfit.1[,2], modelfit.1[,3], type = 'l', lty = 2)
legend(5, -2, c("PD = 0", "PD = 1"), lty = c(1, 2)) 

Figure 12.5, page 437. To generate these graphs, we will first fit a logistic regression model that predicts the event using three age terms, parental divorce, and sex.
depression<-read.table("http://www.ats.ucla.edu/stat/examples/alda/depression_pp.csv", sep=",", header=T)

dmodel<-glm(EVENT ~ ONE + age_18 + age_18sq + age_18cub + PD + FEMALE - 1, family=binomial(link = "logit"), data = depression)

Next, we can to use the model coefficients to generate predicted values for all of the time/PD/sex combinations we are interested in.  We will do this in a loop.  Within each PD/sex combination, we are also interested in the "survival" rate at each time point.  This will be the proportion of subjects that have not yet experienced the event at a given time point. Both the hazard and the survival rate will be calculated in the loop using the predicted values from the model. We will then subset the data by PD/sex combination, resulting in four datasets with one observation for each time point.

modelfit<- c()
for (i in 0:1){
  for (k in 0:1){
  survivor <- 1
    for (j in 4:39){
    fitx = dmodel$coef[1] + dmodel$coef[2]*(j-18) + dmodel$coef[3]*(j-18)^2 + dmodel$coef[4]*(j-18)^3 + dmodel$coef[5]*i + dmodel$coef[6]*k
    hazard = 1/(1 + exp(-fitx))
    survivor = (1-hazard)*survivor
    modelfit <- rbind(modelfit, c(i, j, k, fitx, hazard, survivor))}}}

colnames(modelfit) <- c("pd", "age", "female", "fit", "hazard", "survival")

modelfit.0male <- subset(data.frame(modelfit), modelfit[,1]==0 & modelfit[,3]==0)
modelfit.0female <- subset(data.frame(modelfit), modelfit[,1]==0 & modelfit[,3]==1)
modelfit.1male <- subset(data.frame(modelfit), modelfit[,1]==1 & modelfit[,3]==0)
modelfit.1female <- subset(data.frame(modelfit), modelfit[,1]==1 & modelfit[,3]==1)

To plot the four curves in each plot, we start with one plot, and then overlay the other three lines in the same plot area.

par(mfrow=c(1,2))

plot(modelfit.0male$age, modelfit.0male$hazard, type = 'l', lty = 1, ylim = c(0, .04), xlab = "Age", ylab = "Fitted hazard")
points(modelfit.0female$age, modelfit.0female$hazard, type = 'l', lty = 2)
points(modelfit.1male$age, modelfit.1male$hazard, type = 'l', lty = 3)
points(modelfit.1female$age, modelfit.1female$hazard, type = 'l', lty = 4)
legend(5, 0.04, c("PD = 0, Male", "PD = 0, Female", "PD = 1, Male", "PD = 1, Female"), lty = c(1, 2, 3, 4))


plot(modelfit.0male$age, modelfit.0male$survival, type = 'l', lty = 1, ylim = c(0, 1), xlab = "Age", ylab = "Fitted survival")
points(modelfit.0female$age, modelfit.0female$survival, type = 'l', lty = 2)
points(modelfit.1male$age, modelfit.1male$survival, type = 'l', lty = 3)
points(modelfit.1female$age, modelfit.1female$survival, type = 'l', lty = 4)
legend(0, 0.2, c("PD = 0, Male", "PD = 0, Female", "PD = 1, Male", "PD = 1, Female"), lty = c(1, 2, 3, 4))


Figure 12.6, page 445. This set of plots is similar to previous examples: we will first calculate the proportion of a subset of subjects who experienced the event at a each time point, and then create plots of the proportion (or some transformation) against time with different lines representing different subject characteristics.


fa<-read.table("http://www.ats.ucla.edu/stat/examples/alda/firstarrest_pp.csv", sep=",", header=T)

We loop through all of the combinations of time, race, and abused status in our dataset, subsetting on each and calculating the proportion of subjects experiencing the event within each subset.  We then take the logit of this proportion. After completing the loop, we subset the percent and logit data into four datasets, one for each race/abuse combination.

percents <- c()
for (i in 8:18){
  for (j in 0:1){
    for (k in 0:1){
      x <- subset(fa, PERIOD==i & BLACK==j & ABUSED==k)
      if (sum(x$EVENT) > 0){
        y <- mean(x$EVENT)}
      if (sum(x$EVENT)==0){
        y <- NA}
      logity <- log(y/(1-y))
      z <- c(i, j, k, y, logity)
      percents <- rbind(percents, z)}}}

colnames(percents) <- c("age", "black", "abused", "pct", "hazard")
percents.nm <- as.data.frame(na.omit(percents))

percents.w.a <- subset(percents.nm, black==0 & abused==1)
percents.b.a <- subset(percents.nm, black==1 & abused==1)
percents.w.na <- subset(percents.nm, black==0 & abused==0)
percents.b.na <- subset(percents.nm, black==1 & abused==0)

Next, we create two plots: one of the white subjects and one of the black subjects.  In each, we first plot a line for one of the abuse categories, and then overlay a line for the other.

par(mfrow=c(1,2))
plot(percents.w.a$age, percents.w.a$hazard, type = "l", lty = 1, ylim = c(-7, -2), xlim = c(8,18),
main = "White", xlab = "Age", ylab = "Sample logit(hazard)")
points(percents.w.na$age, percents.w.na$hazard, type = "l", lty = 2)
legend(8, -2, c("Abused", "Not Abused"), lty = c(1, 2))

plot(percents.b.a$age, percents.b.a$hazard, type = "l", lty = 1, ylim = c(-7, -2), xlim = c(8,18),
main = "Black", xlab = "Age", ylab = "Sample logit(hazard)")
points(percents.b.na$age, percents.b.na$hazard, type = "l", lty = 2)
legend(8, -2, c("Abused", "Not Abused"), lty = c(1, 2))

The first two graphs were created using raw proportions.  The third graph is based on a logistic model fitted to the data, using time, race, abuse status, and an abuse-race interaction term as predictors.  We will run this model and then use the coefficients to predict hazard rates at each time point for each race/abuse combination.  These predicted values will be generated using a loop, as we have done in previous examples. 

famodel<-glm(EVENT ~ D8 + D9 + D10 + D11 + D12 + D13 + D14 + D15 + D16 + D17 + D18 + ABUSED + BLACK + ABLACK - 1,
family=binomial(link = "logit"), data = fa)

modelfit <- c()
for (j in 0:1){
  for (k in 0:1){
  survivor <- 1
    for (i in 8:18){
      logitfit <- famodel$coef[i-7] + famodel$coef[12]*j + famodel$coef[13]*k + famodel$coef[14]*j*k
      hazard = 1/(1 + exp(-logitfit))
      survivor = (1-hazard)*survivor
      z <- c(i, j, k, j*k, logitfit, hazard, survivor)
      modelfit <- rbind(modelfit, z)}}}

After completing the loop, we will subset our data as we did for the first two plots. Then, we create one graph by plotting the line for one of the four race/abuse combinations and then overlaying the other three lines. 


colnames(modelfit) <- c("age", "abused", "black", "ablack", "logitfit", "hazard", "survival")
modelfit <- as.data.frame(modelfit)
modelfit.w.a <- subset(modelfit, black==0 & abused==1)
modelfit.b.a <- subset(modelfit, black==1 & abused==1)
modelfit.w.na <- subset(modelfit, black==0 & abused==0)
modelfit.b.na <- subset(modelfit, black==1 & abused==0)

par(mfrow=c(1,1))
plot(modelfit.w.a$age, modelfit.w.a$logitfit, type = "l", lty = 1, ylim = c(-8, -2), xlim = c(8,18), xlab = "Age", ylab = "Fitted logit(hazard)")
points(modelfit.w.na$age, modelfit.w.na$logitfit, type = "l", lty = 2)
points(modelfit.b.a$age, modelfit.b.a$logitfit, type = "l", lty = 3)
points(modelfit.b.na$age, modelfit.b.na$logitfit, type = "l", lty = 4)
legend(8, -2, c("White Abused", "White Not Abused", "Black Abused", "Black Not Abused"), lty = c(1, 2, 3, 4))

Table 12.4, page 449. The parameter estimates and fit statistics in this table can be found in the output from each of the three models below. 

depression<-read.table("http://www.ats.ucla.edu/stat/examples/alda/depression_pp.csv", sep=",", header=T)

modelA<-glm(EVENT ~ ONE + age_18 + age_18sq + age_18cub + PD + FEMALE + NSIBS- 1,
  family=binomial(link = "logit"), data = depression)
summary(modelA)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
ONE       -4.359e+00  1.216e-01 -35.845  < 2e-16 ***
age_18     6.106e-02  1.166e-02   5.236 1.64e-07 ***
age_18sq  -7.308e-03  1.224e-03  -5.972 2.34e-09 ***
age_18cub  1.816e-04  7.897e-05   2.300 0.021447 * 
PD         3.726e-01  1.624e-01   2.295 0.021754 * 
FEMALE     5.587e-01  1.095e-01   5.103 3.34e-07 ***
NSIBS     -8.141e-02  2.227e-02  -3.655 0.000257 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 51288.7  on 36997  degrees of freedom
Residual deviance:  4124.3  on 36990  degrees of freedom
AIC: 4138.3

modelB<-glm(EVENT ~ ONE + age_18 + age_18sq + age_18cub + PD + FEMALE + SIBS12 + SIBS34 + SIBS56 + SIBS78 + SIBS9PLUS - 1,
  family=binomial(link = "logit"), data = depression)
summary(modelB)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
ONE       -4.500e+00  2.067e-01 -21.773  < 2e-16 ***
age_18     6.145e-02  1.166e-02   5.269 1.37e-07 ***
age_18sq  -7.289e-03  1.224e-03  -5.958 2.56e-09 ***
age_18cub  1.813e-04  7.897e-05   2.296   0.0217 * 
PD         3.727e-01  1.625e-01   2.294   0.0218 * 
FEMALE     5.596e-01  1.095e-01   5.109 3.24e-07 ***
SIBS12     2.085e-02  1.976e-01   0.106   0.9160   
SIBS34     1.076e-02  2.100e-01   0.051   0.9591   
SIBS56    -4.942e-01  2.545e-01  -1.942   0.0522 . 
SIBS78    -7.754e-01  3.437e-01  -2.256   0.0241 * 
SIBS9PLUS -6.585e-01  3.440e-01  -1.914   0.0556 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 51289  on 36997  degrees of freedom
Residual deviance:  4118  on 36986  degrees of freedom
AIC: 4140

modelC<-glm(EVENT ~ ONE + age_18 + age_18sq + age_18cub + PD + FEMALE + BIGFAMILY - 1,
  family=binomial(link = "logit"), data = depression)
summary(modelC)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   
ONE       -4.483e+00  1.087e-01 -41.237  < 2e-16 ***
age_18     6.141e-02  1.166e-02   5.265 1.40e-07 ***
age_18sq  -7.291e-03  1.224e-03  -5.959 2.54e-09 ***
age_18cub  1.815e-04  7.897e-05   2.298   0.0215 * 
PD         3.710e-01  1.623e-01   2.286   0.0222 * 
FEMALE     5.581e-01  1.095e-01   5.098 3.44e-07 ***
BIGFAMILY -6.108e-01  1.446e-01  -4.225 2.39e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 51288.7  on 36997  degrees of freedom
Residual deviance:  4118.8  on 36990  degrees of freedom
AIC: 4132.8


Figure 12.8, page 458. This example compares the hazard rates of male and female subjects across time for unfitted values, and then for fitted values from several different models. Note that it would be possible to generate all of these fitted and unfitted values in a single loop after running the three models.  We have opted to keep these steps separate for easier reading. 


math<-read.table("http://www.ats.ucla.edu/stat/examples/alda/mathdropout_pp.csv", sep=",", header=T)

First, for each combination of time and sex, we use a loop to generate the proportion of subjects experiencing the event and calculate the logit of the proportion. Then we subset this proportion and logit data by sex.

percents <- c()
for (i in 1:5){
  for (j in 0:1){
    x <- subset(math, PERIOD==i & FEMALE==j)
    if (sum(x$EVENT) > 0){
      y <- mean(x$EVENT)}
    if (sum(x$EVENT)==0){
      y <- NA}
    logity <- log(y/(1-y))
    z <- c(i, j, y, logity)
    percents <- rbind(percents, z)}}

colnames(percents) <- c("term", "female", "pct", "logit")
percents <- as.data.frame(na.omit(percents))

percents.m <- subset(percents, female==0)
percents.f <- subset(percents, female==1)

Next, we run the first of three models, Model A.  Then, using a loop, we generate predicted proportion values for each combination of time and sex with the coefficients from Model A and calculate the logit of the fitted proportion.  Then we subset this proportion and logit data by sex.

modelA<-glm(EVENT ~ HS11 + HS12 + COLL1 + COLL2 + COLL3 + FEMALE - 1,
  family=binomial(link = "logit"), data = math)


modelfitA <- c()
for (i in 1:5){
  for (j in 0:1){
    logitfit <- modelA$coef[i] + modelA$coef[6]*j
    hazard = 1/(1 + exp(-logitfit))
    z <- c(i, j, logitfit, hazard)
    modelfitA <- rbind(modelfitA, z)}}

colnames(modelfitA) <- c("term", "female", "logitfit", "hazard")
modelfitA <- as.data.frame(modelfitA)
modelfitA.m <- subset(modelfitA, female==0)
modelfitA.f <- subset(modelfitA, female==1)

Next, we run the second of three models, Model B, and go through the same steps as for Model A. 

modelB <- glm(EVENT ~ HS11 + HS12 + COLL1 + COLL2 + COLL3 + FHS11 + FHS12 + FCOLL1 + FCOLL2 + FCOLL3 - 1,
  family=binomial(link = "logit"), data = math)

modelfitB <- c()
for (i in 1:5){
  for (j in 0:1){
    logitfit <- modelB$coef[i] + modelB$coef[i+5]*j
    hazard = 1/(1 + exp(-logitfit))
    z <- c(i, j, logitfit, hazard)
    modelfitB <- rbind(modelfitB, z)}}

colnames(modelfitB) <- c("term", "female", "logitfit", "hazard")
modelfitB <- as.data.frame(modelfitB)
modelfitB.m <- subset(modelfitB, female==0)
modelfitB.f <- subset(modelfitB, female==1)

Next, we run the last of three models, Model C, and go through the same steps as for Model A and Model B. 

modelC <- glm(EVENT ~ HS11 + HS12 + COLL1 + COLL2 + COLL3 + FEMALE + FLTIME - 1,
  family=binomial(link = "logit"), data = math)

modelfitC <- c()
for (i in 1:5){
  for (j in 0:1){
    logitfit <- modelC$coef[i] + modelC$coef[6]*j + modelC$coef[7]*i*j
    hazard = 1/(1 + exp(-logitfit))
    z <- c(i, j, logitfit, hazard)
    modelfitC <- rbind(modelfitC, z)}}

colnames(modelfitC) <- c("term", "female", "logitfit", "hazard")
modelfitC <- as.data.frame(modelfitC)
modelfitC.m <- subset(modelfitC, female==0)
modelfitC.f <- subset(modelfitC, female==1)

Lastly, we create four plots: one for the unfitted hazard and one for each of the three fitted hazards. For each graph, we start by plotting one line for males and then overlay a line for females.  

par(mfrow=c(2,2))

plot(percents.m$term, percents.m$logit, type = "l", lty = 1, ylim = c(-2.5, 0),
main = "Within-group sample hazard functions", xlab = "Term", ylab = "Sample logit(hazard)")
points(percents.f$term, percents.f$logit, type = "l", lty = 2)
legend("bottomright", c("Male", "Female"), lty = c(1, 2))


plot(modelfitA.m$term, modelfitA.m$logitfit, type = "l", lty = 1, ylim = c(-2.5, 0),
main = "Model A", xlab = "Term", ylab = "Fitted logit(hazard)")
points(modelfitA.f$term, modelfitA.f$logitfit, type = "l", lty = 2)
legend("bottomright", c("Male", "Female"), lty = c(1, 2))

plot(modelfitB.m$term, modelfitB.m$logitfit, type = "l", lty = 1, ylim = c(-2.5, 0),
main = "Model B", xlab = "Term", ylab = "Fitted logit(hazard)")
points(modelfitB.f$term, modelfitB.f$logitfit, type = "l", lty = 2)
legend("bottomright", c("Male", "Female"), lty = c(1, 2))


plot(modelfitC.m$term, modelfitC.m$logitfit, type = "l", lty = 1, ylim = c(-2.5, 0),
main = "Model C", xlab = "Term", ylab = "Fitted logit(hazard)")
points(modelfitC.f$term, modelfitC.f$logitfit, type = "l", lty = 2)
legend("bottomright", c("Male", "Female"), lty = c(1, 2))

Table 12.5, page 459. The parameter estimates and fit statistics in this table can be found in the output from each of the three models below. 


modelA<-glm(EVENT ~ HS11 + HS12 + COLL1 + COLL2 + COLL3 + FEMALE - 1,
family=binomial(link = "logit"), data = math)
summary(modelA)

Coefficients:
       Estimate Std. Error z value Pr(>|z|)   
HS11   -2.13080    0.05674 -37.555  < 2e-16 ***
HS12   -0.94248    0.04789 -19.681  < 2e-16 ***
COLL1  -1.44948    0.06345 -22.845  < 2e-16 ***
COLL2  -0.61765    0.07571  -8.158 3.42e-16 ***
COLL3  -0.77160    0.14280  -5.403 6.54e-08 ***
FEMALE  0.37865    0.05015   7.551 4.33e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13250.2  on 9558  degrees of freedom
Residual deviance:  9804.3  on 9552  degrees of freedom
AIC: 9816.3
modelB <- glm(EVENT ~ HS11 + HS12 + COLL1 + COLL2 + COLL3 + FHS11 + FHS12 + FCOLL1 + FCOLL2 + FCOLL3 - 1,
family=binomial(link = "logit"), data = math)
summary(modelB)

Coefficients:
       Estimate Std. Error z value Pr(>|z|)   
HS11   -2.00767    0.07148 -28.087  < 2e-16 ***
HS12   -0.96425    0.05854 -16.470  < 2e-16 ***
COLL1  -1.48240    0.08471 -17.499  < 2e-16 ***
COLL2  -0.71001    0.10073  -7.048 1.81e-12 ***
COLL3  -0.86904    0.19077  -4.555 5.23e-06 ***
FHS11   0.15680    0.09777   1.604 0.108791   
FHS12   0.41866    0.07924   5.284 1.27e-07 ***
FCOLL1  0.44068    0.11580   3.805 0.000142 ***
FCOLL2  0.57075    0.14455   3.948 7.86e-05 ***
FCOLL3  0.60077    0.28573   2.103 0.035502 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13250.2  on 9558  degrees of freedom
Residual deviance:  9796.3  on 9548  degrees of freedom
AIC: 9816.3

modelC <- glm(EVENT ~ HS11 + HS12 + COLL1 + COLL2 + COLL3 + FEMALE + FLTIME - 1,
family=binomial(link = "logit"), data = math)
summary(modelC)

Coefficients:
       Estimate Std. Error z value Pr(>|z|)   
HS11   -2.04592    0.06464 -31.649  < 2e-16 ***
HS12   -0.92551    0.04821 -19.197  < 2e-16 ***
COLL1  -1.49655    0.06649 -22.507  < 2e-16 ***
COLL2  -0.71781    0.08605  -8.341  < 2e-16 ***
COLL3  -0.91655    0.15570  -5.887 3.94e-09 ***
FEMALE  0.22749    0.07744   2.937  0.00331 **
FLTIME  0.11977    0.04698   2.549  0.01080 * 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13250.2  on 9558  degrees of freedom
Residual deviance:  9797.8  on 9551  degrees of freedom
AIC: 9811.8

Table 12.6, page 465. This table shows some summary information for a set of eight subjects. To generate this table in R, we wrote a function that takes an ID number as an argument and returns a vector containing the information seen in a given row of this table.


firstsex<-read.table("http://www.ats.ucla.edu/stat/examples/alda/firstsex_pp.csv", sep=",", header=T)

First, we fit a logistic model to the data predicting the event using time, PT, and PAS. From this model, we then generated a deviance residual for each observation.

model12.6<-glm(event ~ d7 + d8 + d9 + d10 + d11 + d12 + pt+pas- 1, family=binomial(link = "logit"), data = firstsex)
firstsex.r <- cbind(firstsex, dev.res = residuals(model12.6, type="deviance"))

Next, we defined a function that limits the dataset to a given subject, and then extracts the subjects last (greatest time point value) values for PT, PAS, grade, and censor. The function then creates a vector containing this information, the deviance residuals from all the time points for which the subject has data, and an SS deviance for the subject.

get12.6 <- function(id.num){
  x <- subset(firstsex.r, id==id.num)
  pt <- max(x$pt)
  pas <- max(x$pas)
  grade <- max(x$period)
  censor <- abs(max(x$event)-1)
  gr <- x$dev.res
  gr.l <- length(gr)
  if (gr.l < 6){
    gr <- c(gr, c(rep(NA, (6-gr.l))))}
  ss.dev <- sum(na.omit(gr*gr))
  z <- c(id.num, pt, pas, grade, censor, gr, ss.dev)
  return(z)}

The table was created by calling the above function with the eight IDs that appear in the table and stacking the eight vectors together.

tab12.6 <- rbind(get12.6(22), get12.6(112), get12.6(166), get12.6(89), get12.6(102), get12.6(87), get12.6(67), get12.6(212))
colnames(tab12.6) <- c("id", "pt", "pas", "grade", "censor", "gr7", "gr8", "gr9", "gr10", "gr11", "gr12", "ss.dev")

tab12.6


      id pt        pas grade censor        gr7        gr8        gr9       gr10       gr11       gr12   ss.dev
[1,]  22  1 -0.6496482    12      0 -0.4117311 -0.2944395 -0.5840395 -0.7176492 -0.7747741  1.4144633 3.713321
[2,] 112  1 -0.6609311    12      1 -0.4110716 -0.2939581 -0.5831423 -0.7165919 -0.7736556 -0.9562858 2.621976
[3,] 166  1  2.7814130    11      0 -0.6614533 -0.4806816 -0.9108230 -1.0903113  1.1913978         NA 4.106381
[4,]  89  0 -0.0751572    11      0 -0.3248377 -0.2313901 -0.4644566 -0.5751775  1.8623690         NA 4.174028
[5,] 102  1  0.6049317     8      0 -0.4913390  2.3694820         NA         NA         NA         NA 5.855859
[6,]  87  1  2.6779010     7      0  1.8176297         NA         NA         NA         NA         NA 3.303778
[7,]  67  1  2.2746520    12      0 -0.6180257 -0.4476515 -0.8559141 -1.0294308 -1.1007321  1.0429686 4.674059
[8,] 212  0 -0.9617912    12      1 -0.2857033 -0.2032141 -0.4097550 -0.5090181 -0.5524318 -0.6958426 1.339299

Figure 12.8, page 467. These two plots are generated using the deviance residuals from a model.  The first plot depicts the residuals from all of the observations, while the second depicts the sum of squared residuals for each ID. We first need to fit the model and generate the residuals.


firstsex<-read.table("http://www.ats.ucla.edu/stat/examples/alda/firstsex_pp.csv", sep=",", header=T)
model12.6<-glm(event ~ d7 + d8 + d9 + d10 + d11 + d12 + pt+pas- 1, family=binomial(link = "logit"), data = firstsex)
firstsex.r <- cbind(firstsex, dev.res = residuals(model12.6, type="deviance"))

Next, we have defined a function that accepts an ID as an argument and calculates the sum of the squared deviance residuals for observations associated with the given ID.

get12.9 <- function(id.num){
  x <- subset(firstsex.r, id==id.num)
  gr <- x$dev.res
  ss.dev <- sum((gr*gr))
  return(ss.dev)}

Having defined the above function, we loop through the values of ID, calling the above function to calculate the sum of squared residuals for each ID. This gives us a dataset containing one observation per ID.

ss.devs <- c()
  for (i in 1:216){
  z <- c(i, get12.9(i))
  if (z[2]!=0){
    ss.devs <- rbind(ss.devs, z)}}

Finally, we create two plots: one of all of the deviance residuals against ID and one of the sum of squared residuals against ID.

par(mfrow=c(2,1))
plot(firstsex.r$id, firstsex.r$dev.res, ylab = "Deviance residual", xlab = "ID")
plot(ss.devs[,1], ss.devs[,2], ylab = "SS Deviance residual", xlab = "ID")

How to cite this page

Report an error on this page or leave a comment

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.