This is a brief introduction to performing longitudinal analyses with R. The method is described in various sources, perhaps the most approachable is Singer & Willett (2003) from which the first example was derived.
If you haven't installed packages then see Appendix I (below).
To use these examples you will need to install my BES package and several packages from CRAN.
If your library statements fail with a "no such package" error then you should try to install the package.
There are two packages in R with which you can perform mixed model analyses, the lme function from the nlme package or the lmer function with the lme4 package.
lme4 is generally preferred because it's presently under active development and tends to be faster.
However, for various pedagogical reasons (Bates: lmer, p-values and all that), lmer does not report degrees of freedom or p-values for fixed effects (though as you'll see below there's a way around this).
All of the examples in Singer & Willett are coded in Mplus, MLwiN, HLM, SAS, Stata, R, and SPSS at UCLA.
I'll reproduce an example from Chapter 5 (Treating time more flexibly) where they examine changes in CES-D by employment status (model C, p 166ff).
They use nlme, but I prefer lme4.
library(lme4)
unemployment <- read.table("http://www.ats.ucla.edu/stat/r/examples/alda/data/unemployment_pp.txt", header = T, sep = ",")
mer = lmer(cesd ~ months * unemp + (months | id), data = unemployment, REML = F)
summary(mer)
Linear mixed model fit by maximum likelihood ['lmerMod']
Formula: cesd ~ months * unemp + (months | id)
Data: unemployment
AIC BIC logLik deviance
5119 5155 -2552 5103
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 93.713 9.681
months 0.451 0.672 -0.60
Residual 62.031 7.876
Number of obs: 674, groups: id, 254
Fixed effects:
Estimate Std. Error t value
(Intercept) 9.617 1.889 5.09
months 0.162 0.194 0.84
unemp 8.529 1.878 4.54
months:unemp -0.465 0.217 -2.14
Correlation of Fixed Effects:
(Intr) months unemp
months -0.888
unemp -0.911 0.863
months:unmp 0.755 -0.878 -0.852
Part of the convenience of using R is that model coefficients are readily accessible and can be used to generate plots.
Singer & Willett illustrate a brute force method for plotting the results (brute force in the sense that you need to recode the plotting commands if the model changes -- later I'll show a different method that's less sensitive to model alterations).
months.seq <- seq(0, 14)
fixef.c <- fixef(mer)
unemp.c1 <- fixef.c[[1]] + fixef.c[[2]] * months.seq + fixef.c[[3]] + fixef.c[[4]] * months.seq
unemp.c0 <- fixef.c[[1]] + fixef.c[[2]] * months.seq
plot(months.seq, unemp.c1, type = "l", ylim = c(5, 20), ylab = "CES-D.hat", xlab = "Months since job loss")
lines(months.seq, unemp.c0, lty = 3)
legend(7, 20, c("Unemp = 1", "Unemp = 0"), lty = c(1, 3))
title("Interaction Between Unemp and Time")
For simple models we can generally "see" the effects and the directions of association. For more complex model, particularly with higher order interactions (especially with time), we've found that plotting the "predicted" values is the only way to visualize the effects.
These data are from repeated administrations of the CVLT. We'll look at changes over time (linear and quadratic) and differences in rates of change between women and men. BLSA does not code sex explicitly because men have ID<5000.
In this analysis I centered age at 65 and coded explicitly for the quadratic term (unnecessary but easier to comprehend).
(load(url("http://handls.nih.gov/zBES/2013-10-21-mixedDemo1.rdata")))
[1] "cvl"
cvl$Sex = (cvl$BLSid < 5000)
cvl$Sex = factor(cvl$Sex, labels = zQ(Women, Men))
cvl$Age = cvl$Age - 65
cvl$AgeSq = cvl$Age^2
table(cvl$Sex)
Women Men
3557 4008
library("lme4")
mer = lmer(CVLtca ~ (Age + AgeSq) * Sex + (Age | BLSid), data = cvl, na.action = na.omit)
summary(mer)
Linear mixed model fit by REML ['lmerMod']
Formula: CVLtca ~ (Age + AgeSq) * Sex + (Age | BLSid)
Data: cvl
REML criterion at convergence: 56643
Random effects:
Groups Name Variance Std.Dev. Corr
BLSid (Intercept) 79.812 8.934
Age 0.148 0.385 0.64
Residual 69.390 8.330
Number of obs: 7529, groups: BLSid, 1854
Fixed effects:
Estimate Std. Error t value
(Intercept) 56.102681 0.406401 138.0
Age -0.478675 0.025301 -18.9
AgeSq -0.010124 0.001080 -9.4
SexMen -5.439701 0.568141 -9.6
Age:SexMen -0.022228 0.034702 -0.6
AgeSq:SexMen -0.000687 0.001525 -0.5
Correlation of Fixed Effects:
(Intr) Age AgeSq SexMen Ag:SxM
Age 0.266
AgeSq -0.343 0.318
SexMen -0.715 -0.190 0.245
Age:SexMen -0.194 -0.729 -0.232 0.211
AgeSq:SexMn 0.243 -0.226 -0.708 -0.364 0.222
b = fixef(mer)
pAge = seq(50, 90)
xAge = pAge - 65
xAgeSq = xAge^2
pW = b[1] + b[2] * xAge + b[3] * xAgeSq
pM = b[1] + b[2] * xAge + b[3] * xAgeSq + b[4] + b[5] * xAge + b[6] * xAgeSq
plot(pAge, pM, type = "l", col = "blue", ylim = c(25, 65), ylab = "CVLtca", xlab = "Age")
lines(pAge, pW, lty = 1, col = "red")
legend(75, 65, zQ(Women, Men), lty = c(1, 1), col = zQ(red, blue))
This method is also "brute force." The code below yields the same plot using zMixHat, a utility function included in the BES package. I also illustrate a few other commands that make the code a little more accessible and the plot a little prettier.
library(BES)
pAge = seq(50, 90)
xAge = pAge - 65
xAgeSq = xAge^2
cvlHatCI = zMixHat(cvl, mer, vary = "Age=xAge, Sex=zQ(Women,Men)", derivCov = "AgeSq=Age**2")
head(cvlHatCI)
Age Sex CVLtca AgeSq hat
1 -15 Women 0 225 61.00
2 -14 Women 0 196 60.82
3 -13 Women 0 169 60.61
4 -12 Women 0 144 60.39
5 -11 Women 0 121 60.14
6 -10 Women 0 100 59.88
par(las = 1, lwd = 2)
with(cvlHatCI[cvlHatCI$Sex == "Men", ], plot(pAge, hat, type = "l", col = "blue", ylim = c(25, 65), ylab = "CVLtca", xlab = "Age"))
with(cvlHatCI[cvlHatCI$Sex == "Women", ], lines(pAge, hat, lty = 1, col = "red"))
legend(75, 65, zQ(Women, Men), lty = c(1, 1), col = zQ(red, blue))
Perhaps more relevant but with only one repeated measure we also have longitudinal data in HANDLS.
In this example I examine change over approximately 5 years by poverty status, race, and sex for BVR total errors.
I'll start with a complete model including the (dreaded) 4-way interaction.
I also apply the cftest function from the multcomp package to generate t-tests and p-values for the fixed effects.
library(lme4)
library(multcomp)
library(zStat)
(load(url("http://handls.nih.gov/zBES/2013-10-22-mixedDemo2.rdata")))
[1] "BVR"
describe(BVR)
BVR
7 Variables 5193 Observations
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
HNDid : HANDLS identification
n missing unique Mean .05 .10 .25 .50 .75 .90 .95
5193 0 3009 8160081647 8111208302 8113218601 8132109401 8161727401 8192756601 8211969461 8221847802
lowest : 8031079801 8031107801 8031107802 8031117901 8031117902, highest: 8224505601 8224515701 8224520701 8224521901 8224521902
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
HNDwave : HANDLS Wave
n missing unique Mean
5193 0 2 1.925
1 (2791, 54%), 3 (2402, 46%)
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Age : Initial age
n missing unique Mean .05 .10 .25 .50 .75 .90 .95
5193 0 2097 50.37 34.00 37.00 43.29 51.00 58.00 62.75 64.82
lowest : 30.00 31.00 32.00 32.86 33.00, highest: 70.30 70.32 70.40 70.57 70.78
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
BVRtot : NeuPsy: BVR Total errors
n missing unique Mean .05 .10 .25 .50 .75 .90 .95
5193 0 30 7.527 0 0 3 7 11 15 18
lowest : 0 1 2 3 4, highest: 25 26 27 29 30
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Sex
n missing unique
5193 0 2
Women (2974, 57%), Men (2219, 43%)
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Race
n missing unique
5193 0 2
White (2073, 40%), AfrAm (3120, 60%)
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
PovStat : Poverty status
n missing unique
5193 0 2
Above (3036, 58%), Below (2157, 42%)
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
merBVR1 = lmer(BVRtot ~ PovStat * Race * Sex * Age + (Age | HNDid), data = BVR, na.action = na.omit)
summary(merBVR1)
Linear mixed model fit by REML ['lmerMod']
Formula: BVRtot ~ PovStat * Race * Sex * Age + (Age | HNDid)
Data: BVR
REML criterion at convergence: 31089
Random effects:
Groups Name Variance Std.Dev. Corr
HNDid (Intercept) 2.32609 1.5252
Age 0.00179 0.0423 1.00
Residual 12.98787 3.6039
Number of obs: 5193, groups: HNDid, 3009
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.6314 1.0677 -0.59
PovStatBelow 0.3545 1.7923 0.20
RaceAfrAm 0.4816 1.4633 0.33
SexMen -2.6248 1.5842 -1.66
Age 0.1373 0.0214 6.41
PovStatBelow:RaceAfrAm -1.9043 2.2881 -0.83
PovStatBelow:SexMen 2.9199 2.8882 1.01
RaceAfrAm:SexMen -1.7018 2.1896 -0.78
PovStatBelow:Age 0.0329 0.0362 0.91
RaceAfrAm:Age 0.0254 0.0291 0.87
SexMen:Age 0.0461 0.0317 1.45
PovStatBelow:RaceAfrAm:SexMen 0.5546 3.6243 0.15
PovStatBelow:RaceAfrAm:Age 0.0167 0.0459 0.36
PovStatBelow:SexMen:Age -0.0580 0.0581 -1.00
RaceAfrAm:SexMen:Age 0.0287 0.0436 0.66
PovStatBelow:RaceAfrAm:SexMen:Age -0.0283 0.0728 -0.39
Correlation of Fixed Effects:
(Intr) PvSttB RcAfrA SexMen Age PvSB:RAA PvSB:SM RcAA:SM PvSB:A RcAA:A SxMn:A PvSB:RAA:SM PSB:RAA:A PSB:SM: RAA:SM:
PovStatBelw -0.596
RaceAfrAm -0.730 0.435
SexMen -0.674 0.402 0.492
Age -0.979 0.583 0.714 0.660
PvSttBl:RAA 0.467 -0.783 -0.640 -0.315 -0.457
PvSttBlw:SM 0.370 -0.621 -0.270 -0.548 -0.362 0.486
RcAfrAm:SxM 0.488 -0.291 -0.668 -0.723 -0.477 0.427 0.397
PvSttBlw:Ag 0.580 -0.979 -0.423 -0.391 -0.592 0.767 0.608 0.283
RacAfrAm:Ag 0.721 -0.429 -0.979 -0.486 -0.736 0.626 0.266 0.654 0.436
SexMen:Age 0.661 -0.394 -0.482 -0.979 -0.675 0.308 0.537 0.708 0.400 0.497
PvSB:RAA:SM -0.295 0.495 0.404 0.437 0.288 -0.631 -0.797 -0.604 -0.484 -0.395 -0.428
PvStB:RAA:A -0.456 0.770 0.620 0.308 0.466 -0.979 -0.478 -0.414 -0.787 -0.633 -0.315 0.618
PvSttB:SM:A -0.361 0.609 0.263 0.534 0.368 -0.477 -0.980 -0.387 -0.622 -0.271 -0.546 0.781 0.489
RcAfrA:SM:A -0.481 0.286 0.653 0.712 0.491 -0.418 -0.391 -0.979 -0.291 -0.667 -0.727 0.592 0.422 0.397
PSB:RAA:SM: 0.288 -0.486 -0.391 -0.426 -0.294 0.618 0.782 0.587 0.496 0.400 0.436 -0.980 -0.631 -0.798 -0.599
cftest(merBVR1)
Simultaneous Tests for General Linear Hypotheses
Fit: lmer(formula = BVRtot ~ PovStat * Race * Sex * Age + (Age | HNDid),
data = BVR, na.action = na.omit)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
(Intercept) == 0 -0.6314 1.0677 -0.59 0.554
PovStatBelow == 0 0.3545 1.7923 0.20 0.843
RaceAfrAm == 0 0.4816 1.4633 0.33 0.742
SexMen == 0 -2.6248 1.5842 -1.66 0.098
Age == 0 0.1373 0.0214 6.41 0.00000000015
PovStatBelow:RaceAfrAm == 0 -1.9043 2.2881 -0.83 0.405
PovStatBelow:SexMen == 0 2.9199 2.8882 1.01 0.312
RaceAfrAm:SexMen == 0 -1.7018 2.1896 -0.78 0.437
PovStatBelow:Age == 0 0.0329 0.0362 0.91 0.364
RaceAfrAm:Age == 0 0.0254 0.0291 0.87 0.382
SexMen:Age == 0 0.0461 0.0317 1.45 0.146
PovStatBelow:RaceAfrAm:SexMen == 0 0.5546 3.6243 0.15 0.878
PovStatBelow:RaceAfrAm:Age == 0 0.0167 0.0459 0.36 0.717
PovStatBelow:SexMen:Age == 0 -0.0580 0.0581 -1.00 0.319
RaceAfrAm:SexMen:Age == 0 0.0287 0.0436 0.66 0.511
PovStatBelow:RaceAfrAm:SexMen:Age == 0 -0.0283 0.0728 -0.39 0.698
(Univariate p values reported)
Using backward elimination we can remove the 4-way interaction because it was nonsignificant. Here's a convenient syntactical trick for specifying all of the interactions up to a specific number of terms (note pluses within parentheses and the ^3 indicates the maximum "ways").
merBVR2 = lmer(BVRtot ~ (PovStat + Race + Sex + Age)^3 + (Age | HNDid), data = BVR, na.action = na.omit)
summary(merBVR2)
Linear mixed model fit by REML ['lmerMod']
Formula: BVRtot ~ (PovStat + Race + Sex + Age)^3 + (Age | HNDid)
Data: BVR
REML criterion at convergence: 31086
Random effects:
Groups Name Variance Std.Dev. Corr
HNDid (Intercept) 2.34510 1.5314
Age 0.00178 0.0422 1.00
Residual 12.98732 3.6038
Number of obs: 5193, groups: HNDid, 3009
Fixed effects:
Estimate Std. Error t value
(Intercept) -0.5113 1.0225 -0.50
PovStatBelow 0.0155 1.5664 0.01
RaceAfrAm 0.2582 1.3468 0.19
SexMen -2.8875 1.4329 -2.02
Age 0.1348 0.0205 6.58
PovStatBelow:RaceAfrAm -1.3548 1.7997 -0.75
PovStatBelow:SexMen 3.7967 1.8001 2.11
PovStatBelow:Age 0.0398 0.0314 1.27
RaceAfrAm:SexMen -1.2028 1.7732 -0.68
RaceAfrAm:Age 0.0300 0.0267 1.12
SexMen:Age 0.0515 0.0286 1.80
PovStatBelow:RaceAfrAm:SexMen -0.8234 0.7267 -1.13
PovStatBelow:RaceAfrAm:Age 0.0054 0.0357 0.15
PovStatBelow:SexMen:Age -0.0760 0.0350 -2.17
RaceAfrAm:SexMen:Age 0.0185 0.0349 0.53
Correlation of Fixed Effects:
(Intr) PvSttB RcAfrA SexMen Age PvSB:RAA PvSB:SM PvSB:A RcAA:SM RcAA:A SxMn:A PSB:RAA:S PSB:RAA:A PSB:SM:
PovStatBelw -0.545
RaceAfrAm -0.700 0.304
SexMen -0.636 0.246 0.390
Age -0.977 0.527 0.681 0.618
PvSttBl:RAA 0.383 -0.703 -0.550 -0.072 -0.366
PvSttBlw:SM 0.242 -0.442 0.063 -0.381 -0.221 0.006
PvSttBlw:Ag 0.526 -0.972 -0.287 -0.228 -0.538 0.674 0.405
RcAfrAm:SxM 0.411 -0.008 -0.589 -0.646 -0.394 0.102 -0.123 -0.012
RacAfrAm:Ag 0.690 -0.294 -0.975 -0.380 -0.706 0.526 -0.080 0.299 0.565
SexMen:Age 0.621 -0.232 -0.377 -0.975 -0.636 0.056 0.350 0.235 0.621 0.392
PvSB:RAA:SM -0.065 0.105 0.112 0.106 0.001 -0.167 -0.246 0.013 -0.181 -0.021 -0.007
PvStB:RAA:A -0.370 0.684 0.523 0.055 0.379 -0.966 0.031 -0.704 -0.070 -0.536 -0.057 0.001
PvSttB:SM:A -0.227 0.420 -0.088 0.356 0.232 0.034 -0.947 -0.432 0.167 0.086 -0.365 -0.010 -0.030
RcAfrA:SM:A -0.402 -0.007 0.568 0.630 0.411 -0.076 0.156 0.009 -0.968 -0.583 -0.647 0.030 0.072 -0.169
cftest(merBVR2)
Simultaneous Tests for General Linear Hypotheses
Fit: lmer(formula = BVRtot ~ (PovStat + Race + Sex + Age)^3 + (Age |
HNDid), data = BVR, na.action = na.omit)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
(Intercept) == 0 -0.5113 1.0225 -0.50 0.617
PovStatBelow == 0 0.0155 1.5664 0.01 0.992
RaceAfrAm == 0 0.2582 1.3468 0.19 0.848
SexMen == 0 -2.8875 1.4329 -2.02 0.044
Age == 0 0.1348 0.0205 6.58 0.000000000046
PovStatBelow:RaceAfrAm == 0 -1.3548 1.7997 -0.75 0.452
PovStatBelow:SexMen == 0 3.7967 1.8001 2.11 0.035
PovStatBelow:Age == 0 0.0398 0.0314 1.27 0.204
RaceAfrAm:SexMen == 0 -1.2028 1.7732 -0.68 0.498
RaceAfrAm:Age == 0 0.0300 0.0267 1.12 0.261
SexMen:Age == 0 0.0515 0.0286 1.80 0.071
PovStatBelow:RaceAfrAm:SexMen == 0 -0.8234 0.7267 -1.13 0.257
PovStatBelow:RaceAfrAm:Age == 0 0.0054 0.0357 0.15 0.879
PovStatBelow:SexMen:Age == 0 -0.0760 0.0350 -2.17 0.030
RaceAfrAm:SexMen:Age == 0 0.0185 0.0349 0.53 0.596
(Univariate p values reported)
Now this is getting a little more interesting. It's nearly impossible to get a feel for the direction and magnitude of the effects because the coefficients are not independent (uncorrelated). Plotting the predicted values by group is an easy way to grasp the overall directions and type of interactions.
pAge = seq(30, 70)
hatBVR1 = zMixHat(BVR, merBVR2, vary = "Age=pAge, PovStat=zQ(abovePovStat,belowPovStat),Race=zQ(White,AfrAm),Sex=zQ(Women,Men)")
par(las = 1, lwd = 2)
HNDcolors = HNDpltColors()
with(hatBVR1[hatBVR1$Sex == "Women" & hatBVR1$PovStat == "abovePovStat" & hatBVR1$Race == "White", ], plot(pAge, hat, lty = 1, col = "red", typ = "l", ylim = c(0, 14), ylab = "BVR", xlab = "Age"))
with(hatBVR1[hatBVR1$Sex == "Women" & hatBVR1$PovStat == "abovePovStat" & hatBVR1$Race == "AfrAm", ], lines(pAge, hat, lty = 2, col = "red"))
with(hatBVR1[hatBVR1$Sex == "Women" & hatBVR1$PovStat == "belowPovStat" & hatBVR1$Race == "White", ], lines(pAge, hat, lty = 3, col = "red"))
with(hatBVR1[hatBVR1$Sex == "Women" & hatBVR1$PovStat == "belowPovStat" & hatBVR1$Race == "AfrAm", ], lines(pAge, hat, lty = 4, col = "red"))
with(hatBVR1[hatBVR1$Sex == "Men" & hatBVR1$PovStat == "abovePovStat" & hatBVR1$Race == "White", ], lines(pAge, hat, lty = 1, col = "blue", typ = "l", ylim = c(0, 14), ylab = "BVR", xlab = "Age"))
with(hatBVR1[hatBVR1$Sex == "Men" & hatBVR1$PovStat == "abovePovStat" & hatBVR1$Race == "AfrAm", ], lines(pAge, hat, lty = 2, col = "blue"))
with(hatBVR1[hatBVR1$Sex == "Men" & hatBVR1$PovStat == "belowPovStat" & hatBVR1$Race == "White", ], lines(pAge, hat, lty = 3, col = "blue"))
with(hatBVR1[hatBVR1$Sex == "Men" & hatBVR1$PovStat == "belowPovStat" & hatBVR1$Race == "AfrAm", ], lines(pAge, hat, lty = 4, col = "blue"))
legend(30, 14, zQ(abovePovWhite, abovePovAfrAm, belowPovWhite, belowPovAfrAm), lty = 1:4, col = "black")
text(30, 10, "Women in red", adj = c(0, 0), col = "red")
text(30, 9.5, "Men in blue", adj = c(0, 0), col = "blue")
Although it's a little hard to discern, it appears as though differences associated with poverty status are larger than those associated with race differences.
We can reanalyze and re-plot, or we can recalculate the predicted scores holding constant race (these aren't necessarily equivalent approaches).
Just for illustration, I'll recalculate the predicted scores with zMixHat and re-plot.
par(las = 1, lwd = 2)
hatBVR2 = zMixHat(BVR, merBVR2, vary = "Age=pAge, PovStat=zQ(abovePovStat,belowPovStat),Sex=zQ(Women,Men)", fixedCov = "Race")
with(hatBVR2[hatBVR2$Sex == "Women" & hatBVR2$PovStat == "abovePovStat", ], plot(pAge, hat, lty = 1, col = "red", typ = "l", ylim = c(0, 14), ylab = "BVR", xlab = "Age"))
with(hatBVR2[hatBVR2$Sex == "Women" & hatBVR2$PovStat == "belowPovStat", ], lines(pAge, hat, lty = 2, col = "red"))
with(hatBVR2[hatBVR2$Sex == "Men" & hatBVR2$PovStat == "abovePovStat", ], lines(pAge, hat, lty = 1, col = "blue"))
with(hatBVR2[hatBVR2$Sex == "Men" & hatBVR2$PovStat == "belowPovStat", ], lines(pAge, hat, lty = 2, col = "blue"))
legend(30, 14, zQ(abovePovWomen, belowPovWomen, abovePovMen, belowPovMen), lty = c(1, 2), col = zQ(red, red, blue, blue))
This is a quick plot. I would probably clean up the legend if I were considering publishing it or using it in a more formal presentation
There is a large and growing larger library of packages in CRAN (http://cran.r-project.org/ -- Comprehensive R Archive Network) that you can install into R (select "packages" under Software on the left-hand menu).
To install lme4, issue the command:
install.packages("lme4")
To use the package, issue the command:
library(lme4)
I also created a small package of useful utility routines called BES. I'm sharing them with you because I use these functions all the time. It's easier to share the functions than it is to rewrite all my code.
Install the BES package by downloading the file from BES archive and issuing this install.package command:
install.packages('BES_0.01.tar.gz', repos=NULL, type='source')
Douglas Bates: lme4: Mixed-effects modeling with R
Venables: Exegeses on Linear Models
Singer, J. D., & Willett, J. B. (2003). Applied Longitudinal Data Analysis: Modeling Change and Event Occurrence (1st ed.). Oxford University Press, USA.