The binaryMM package allows users to fit marginalized transition and latent variables (mTLV) models for binary longitudinal data. The aim of this vignette is to provide an overview of the models together with example code and analyses.

The Marginalized Transition and Latent Variables Model

Let N be the total number of subjects, Y_i be the n_i−vector of binary responses for subject i, X_i be the n_i × p matrix of covariates and U_i ∼ N(0, 1). The marginalized transition and latent variable (mTLV) model described in Schildcrout and Heagerty (2007) can be defined by two equations:

logit(μ_ij^m) = β^TX_i logit(μ_ij^c) = Δ_ij(X_i) + γ(X_i)Y_{i(j − 1)} + σ(X_i)U_i

where μ_ij^m = E[Y_ij|X_i] and μ_ij^c = E[Y_ij|X_i, Y_{i(j − 1)}, U_i].

The first equation describes the marginal mean model and the relationship between the outcome Y_i and the covariates X_i. The second equation describes the conditional mean model (also named the dependence model) and the relationship between the outcome Y_i measured over time for each subject i. In particular, the conditional model includes a short-term transition component γ(X_i)Y_{i(j − 1)}, and a random intercept term, σ(X_i)U_i, describing long-term non-diminishing dependence.

Δ_ij(X_i) is a function of the marginal mean, μ_ij^m, and the conditional mean, μ_ij^c, such that the two model above are cohesive. In particular, Δ_ij(X_i) is the value that satisfies the convolution equation:

μ_ij^m = E_{U_i, Y_{i(j − 1)}}(μ_ij^c) = E_{Z_i}[E_{Y_{i(j − 1)}}[logit⁻¹(Δ_ij(X_i) + γ(X_i)Y_{i(j − 1)} + σ(X_i)U_i)]] Δ_ij(X_i) in mTLV is analytically intractable and its value is computed iteratively with a Newton-Raphson method.

Detailed information on marginalized models with transition and/or latent terms can be found in Heagerty (2002), Heagerty and Zeger (1999) and Schildcrout and Heagerty (2007).

Basic Examples

The next two sections explain how different specifications of mTLV models can be fitted using the binaryMM package. The data used are part of the package.

library(binaryMM)

The Madras Longitudinal Schizophrenia Study

madras contains a subset of the data from the Madras Longitudinal Schizophrenia Study Diggle et al. (2002), which collected monthly symptom data on 86 schizophrenia patients after their initial hospitalization. The dataframe has 922 observations on 86 patients and includes the variables:

though. An indicator for thought disorders
age. An indicator for age-at-onset ≥ 20 years
gender. An indicator for female gender
month. Months since hospitalization
id. A unique patient identifiers

The primary question of interest is whether subjects with an older age-at-onset tend to recover more or less quickly, and whether female patients recover more or less quickly. Recovery is measured by a reduction in the presentation of symptoms.

data(madras)
str(madras)
#> 'data.frame':    922 obs. of  5 variables:
#>  $ id     : int  1 1 1 1 1 1 1 1 1 1 ...
#>  $ thought: int  1 1 1 1 1 0 0 0 0 0 ...
#>  $ month  : int  0 1 2 3 4 5 6 7 8 9 ...
#>  $ gender : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ age    : num  1 1 1 1 1 1 1 1 1 1 ...

The marginal mean model is defined as:

logit(μ_ij^m) = β₀ + β₁month_ij + β₂age_i + β₃gender_i + β₄age_i × month_ij + β₅gender_i × month_ij

Multiple dependence models are explored to demonstrate how the mm function can be used. The different dependence models are declared by changing the t.formula and lv.formula arguments. Note that by default formula both are initially assigned NULL and if neither association models are specified, then an error is returned.

Case 1. A dependence model with a transition term only:

logit(μ_ij^c) = Δ_ij(X_i) + γY_{i(j − 1)}

mod.mt <- mm(thought ~ month*gender + month*age, t.formula = ~1,  
             data = madras, id = id)
summary(mod.mt)
#> 
#> Class:
#> MMLong
#> 
#> Call:
#> mm(mean.formula = thought ~ month * gender + month * age, t.formula = ~1, 
#>     id = id, data = madras)
#> 
#> Information Criterion:
#>       AIC        BIC     logLik   Deviance  
#>  688.3789   705.5594  -337.1895   674.3789  
#> 
#> Marginal Mean Parameters:
#>               Estimate  Model SE Chi Square  Pr(>Chi)
#> (Intercept)   1.183683  0.444318     7.0971  0.007721
#> month        -0.342857  0.081841    17.5501 2.798e-05
#> gender       -0.141884  0.416152     0.1162  0.733147
#> age          -0.649770  0.449183     2.0925  0.148021
#> month:gender -0.143788  0.081853     3.0859  0.078975
#> month:age     0.111555  0.085896     1.6867  0.194040
#> 
#> Dependence Model Parameters:
#>                   Estimate Model SE Chi Square  Pr(>Chi)
#> gamma:(Intercept)  3.16583  0.23014     189.23 < 2.2e-16
#> 
#> Number of clusters:             86 
#> Maximum cluster size:           12 
#> Convergence status (nlm code):  1 
#> Number of iterations:           22

Case 2. A dependence model with a latent term only:

logit(μ_ij^c) = Δ_ij(X_i) + σU_i

mod.mlv <- mm(thought ~ month*gender + month*age, lv.formula = ~1, 
              data = madras, id = id)
summary(mod.mlv)
#> 
#> Class:
#> MMLong
#> 
#> Call:
#> mm(mean.formula = thought ~ month * gender + month * age, lv.formula = ~1, 
#>     id = id, data = madras)
#> 
#> Information Criterion:
#>       AIC        BIC     logLik   Deviance  
#>  750.5767   767.7571  -368.2883   736.5767  
#> 
#> Marginal Mean Parameters:
#>               Estimate  Model SE Chi Square  Pr(>Chi)
#> (Intercept)   1.569015  0.399400    15.4326 8.550e-05
#> month        -0.399007  0.062845    40.3107 2.166e-10
#> gender       -0.539932  0.355893     2.3016   0.12924
#> age          -0.911165  0.393487     5.3621   0.02058
#> month:gender -0.081899  0.060179     1.8521   0.17354
#> month:age     0.140215  0.063107     4.9366   0.02629
#> 
#> Dependence Model Parameters:
#>                        Estimate Model SE Chi Square  Pr(>Chi)
#> log(sigma):(Intercept)  0.81289  0.12764     40.559 1.908e-10
#> 
#> Number of clusters:             86 
#> Maximum cluster size:           12 
#> Convergence status (nlm code):  1 
#> Number of iterations:           42

Case 3. A dependence model with a transition and a latent term:

logit(μ_ij^c) = Δ_ij(X_i) + γY_{i(j − 1)} + σU_i

mod.mtlv <- mm(thought ~ month*gender + month*age,
               t.formula = ~1, lv.formula = ~1, 
               data = madras, id = id)
summary(mod.mtlv)
#> 
#> Class:
#> MMLong
#> 
#> Call:
#> mm(mean.formula = thought ~ month * gender + month * age, lv.formula = ~1, 
#>     t.formula = ~1, id = id, data = madras)
#> 
#> Information Criterion:
#>       AIC        BIC     logLik   Deviance  
#>  680.1283   699.7631  -332.0642   664.1283  
#> 
#> Marginal Mean Parameters:
#>               Estimate  Model SE Chi Square  Pr(>Chi)
#> (Intercept)   1.327440  0.434588     9.3299  0.002255
#> month        -0.367564  0.077443    22.5268 2.072e-06
#> gender       -0.282497  0.402015     0.4938  0.482241
#> age          -0.732176  0.436180     2.8177  0.093228
#> month:gender -0.111740  0.078306     2.0362  0.153591
#> month:age     0.117705  0.080870     2.1184  0.145536
#> 
#> Dependence Model Parameters:
#>                        Estimate Model SE Chi Square Pr(>Chi)
#> gamma:(Intercept)      2.511119 0.303963    68.2486   <2e-16
#> log(sigma):(Intercept) 0.074494 0.244870     0.0925    0.761
#> 
#> Number of clusters:             86 
#> Maximum cluster size:           12 
#> Convergence status (nlm code):  1 
#> Number of iterations:           50

Case 4. A dependence model with a transition term that is modified by gender.

logit(μ_ij^c) = Δ_ij(X_i) + (γ₀ + γ₁gender_i)Y_{i(j − 1)}

mod.mtgender <- mm(thought ~ month*gender + month*age,
                   t.formula = ~gender, data = madras, id = id)
summary(mod.mtgender)
#> 
#> Class:
#> MMLong
#> 
#> Call:
#> mm(mean.formula = thought ~ month * gender + month * age, t.formula = ~gender, 
#>     id = id, data = madras)
#> 
#> Information Criterion:
#>       AIC        BIC     logLik   Deviance  
#>  690.2915   709.9263  -337.1458   674.2915  
#> 
#> Marginal Mean Parameters:
#>               Estimate  Model SE Chi Square  Pr(>Chi)
#> (Intercept)   1.175157  0.444200     6.9990  0.008156
#> month        -0.342431  0.081702    17.5662 2.775e-05
#> gender       -0.137202  0.416448     0.1085  0.741809
#> age          -0.635415  0.452122     1.9752  0.159901
#> month:gender -0.143961  0.082080     3.0763  0.079443
#> month:age     0.110312  0.085973     1.6464  0.199456
#> 
#> Dependence Model Parameters:
#>                   Estimate Model SE Chi Square Pr(>Chi)
#> gamma:(Intercept)  3.11501  0.28599    118.634   <2e-16
#> gamma:gender       0.14321  0.48550      0.087    0.768
#> 
#> Number of clusters:             86 
#> Maximum cluster size:           12 
#> Convergence status (nlm code):  1 
#> Number of iterations:           34

Case 5. A dependence model with a latent term that is modified by gender. Note that because σ is a positive quantity, to fit a mTLV model where the latent term is modified by gender, we need to specify two indicator variables: I0 for gender = 0, and I1 for gender = 1. The model to be specified in lv.fomula will take the form: ~0+I0+I1.

logit(μ_ij^c) = Δ_ij(X_i) + [σ₀I(gender_i = = 0) + σ₁I(gender_i = = 1)]U_i

# set-up two new indicator variables for gender
madras$g0    <- ifelse(madras$gender == 0, 1, 0)
madras$g1    <- ifelse(madras$gender == 1, 1, 0)
mod.mlvgender <- mm(thought ~ month*gender + month*age,
                   lv.formula = ~0+g0+g1, data = madras, id = id)
summary(mod.mlvgender)
#> 
#> Class:
#> MMLong
#> 
#> Call:
#> mm(mean.formula = thought ~ month * gender + month * age, lv.formula = ~0 + 
#>     g0 + g1, id = id, data = madras)
#> 
#> Information Criterion:
#>       AIC        BIC     logLik   Deviance  
#>  752.4992   772.1340  -368.2496   736.4992  
#> 
#> Marginal Mean Parameters:
#>               Estimate  Model SE Chi Square  Pr(>Chi)
#> (Intercept)   1.566764  0.400938    15.2705 9.316e-05
#> month        -0.395507  0.064395    37.7227 8.155e-10
#> gender       -0.515176  0.366086     1.9804   0.15935
#> age          -0.921331  0.395394     5.4296   0.01980
#> month:gender -0.091526  0.069640     1.7273   0.18875
#> month:age     0.140228  0.063187     4.9251   0.02647
#> 
#> Dependence Model Parameters:
#>               Estimate Model SE Chi Square  Pr(>Chi)
#> log(sigma):g0  0.84429  0.17128     24.297 8.257e-07
#> log(sigma):g1  0.77199  0.19523     15.636 7.678e-05
#> 
#> Number of clusters:             86 
#> Maximum cluster size:           12 
#> Convergence status (nlm code):  1 
#> Number of iterations:           50

The parameters from the marginal mean model have the same interpretation regardless of the dependence model used. Overall, older individuals tend to have slower recovery time than younger subjects, while females recover quicker than males.

Weighted Likelihood

The binaryMM package allows user to add sampling weights and estimates the parameters of interest in those cases where the available sample might not be representative of the target population (i.e., survey data). This section shows how the sampling weights can be added in the mm syntax using the datarand dataframe.

The dataframe has 24,999 observation on 2,500 subjects and includes the variables:

id. A unique patient identifier
Y. A binary longitudinal outcome
time. A continuous time-varying covariate indicating time of each follow-up
binary. A binary time-fixed covariate indicating whether a patient was assigned to a treatment arm (1) or a control arm (0)

data(datrand)
str(datrand)
#> 'data.frame':    24999 obs. of  4 variables:
#>  $ id    : int  1 1 1 1 1 1 1 1 1 2 ...
#>  $ Y     : int  0 0 1 1 0 0 0 0 0 0 ...
#>  $ time  : num  0 1 2 3 4 5 6 7 8 0 ...
#>  $ binary: num  0 0 0 0 0 0 0 0 0 0 ...

From datarand a biased sampled can be created by assuming that complete data are available only for 1) every one who experienced the event Y at least once, and 2) 20% of the subjects who never experienced the event Y.

# create the sampling scheme
Ymean     <- tapply(datrand$Y, FUN = mean, INDEX = datrand$id)
some.id   <- names(Ymean[Ymean != 0])
none.id   <- names(Ymean)[!(names(Ymean) %in% some.id)]
samp.some <- some.id[rbinom(length(none.id), 1, 1) == 1]
samp.none <- none.id[rbinom(length(none.id), 1, 0.20) == 1]

# sample subjects and create a weight vector
datrand$sampled <- ifelse(datrand$id %in% c(samp.none, samp.some), 1, 0)
dat.small       <- subset(datrand, sampled == 1)
wt              <- ifelse(dat.small$id %in% samp.none, 1/1, 1/0.2)

# fit the mTLV model
mod.wt          <- mm(Y ~ time*binary, t.formula = ~1, data = dat.small, 
                      id = id, weight = wt)
summary(mod.wt)
#> Warning in summary.MMLong(mod.wt): When performing a weighted likelihood
#> analysis (by specifying the weight argument), robust standard errors are
#> reported. Model based standard errors will not be correct and should not be
#> used.
#> 
#> Class:
#> MMLong
#> 
#> Call:
#> mm(mean.formula = Y ~ time * binary, t.formula = ~1, id = id, 
#>     data = dat.small, weight = wt)
#> 
#> Information Criterion:
#>       AIC        BIC     logLik   Deviance  
#>  74718.33   74745.52  -37354.17   74708.33  
#> 
#> Marginal Mean Parameters:
#>              Estimate Robust SE Chi Square  Pr(>Chi)
#> (Intercept) -1.011814  0.045344    497.930 < 2.2e-16
#> time        -0.161071  0.010257    246.591 < 2.2e-16
#> binary       0.319141  0.072757     19.241 1.152e-05
#> time:binary  0.158303  0.013989    128.055 < 2.2e-16
#> 
#> Dependence Model Parameters:
#>                   Estimate Robust SE Chi Square  Pr(>Chi)
#> gamma:(Intercept)  1.04156   0.04267     595.83 < 2.2e-16
#> 
#> Number of clusters:             1697 
#> Maximum cluster size:           15 
#> Convergence status (nlm code):  1 
#> Number of iterations:           27

Note that when the weight argument is specified, model-based standard error will not be correct and should not be reported. Thus, the software will return robust standard errors only together with a warning message.

Functions Available in the Package

The two examples above showed how different mTLV model can be used using simulated data as well as data from the Madras Longitudinal Schizophrenia Study. The table below summarizes the functions in mm available to the user.

Function	Description
`GenBinaryY`	Generate binary response variable under a user-specified mTLV model. The outcome is generated from a Bernoulli distribution where the probability of success is computed as the inverse-logit of the conditional mean. The function requires the user to specify the mean model formula (`mean.formula`) in which a binary covariate is regressed on covariates, one or both components of the dependence model (the latent variable component `lv.formula` or the transition term component `t.formula`), the vector of cluster identifiers (`id`), a vector of values for the parameters of the mean model (`beta`), a vector of values for the parameters of the transition component of the dependence model (`gamma`), a vector of values for the latent component of the dependence model (`sigma`), a dataframe (`data`) with the mean model covariates (ordered by id and time) and a string of the mane of the new binary variable (`Yname`). The function returns the entire `data` object with an additional column `Yname` of the binary longitudinal outcome
`mm`	Fit mTLV model. The function requires the user to specify the mean model formula (`mean.formula`) in which a binary covariate is regressed on covariates, one or both components of the dependence model (the latent variable component `lv.formula` or the transition term component `t.formula`), the vector of cluster identifiers (`id`), and the dataframe to use (`data`). Users can additionally specify the sampling weights (`weight`) to estimate the parameters using weighted likelihood.
`summary`	Summarize the results of a class `MMLong` generated using `mm`. Tables with estimated parameters, standard errors and p-value are printed for both the mean model and the dependence model parameters
`anova`	Allows to compare two nested models of class `MMLong` generated using `mm`. fits using mTLV. Currently comparison can be made for two models only

Reference

Diggle, P, PJ Heagerty, KY Liang, and Zeger SL. 2002. Analysis of Longitudinal Data. Oxford University Press.

Heagerty, PJ. 2002. “Marginalized Transition Models and Likelihood Inference for Longitudinal Categorical Data.” Biometrics 58: 342–51.

Heagerty, PJ, and SL Zeger. 1999. “Marginalized Multilevel Models and Likelihood Inference.” Statistical Science 15 (1): 1–19.

Schildcrout, JS, and PJ Heagerty. 2007. “Marginalized Models for Moderate to Long Series of Longitudinal Binary Response Data.” Biometrics 63 (2): 322–31.