Version info: Code for this page was tested in Stata 12.1
Mixed effects logistic regression is used to model binary outcome variables, in which the log odds of the outcomes are modeled as a linear combination of the predictor variables when data are clustered or there are both fixed and random effects.
Please note: The purpose of this page is to show how to use various data analysis commands. It does not cover all aspects of the research process which researchers are expected to do. In particular, it does not cover data cleaning and checking, verification of assumptions, model diagnostics or potential follow-up analyses.
Examples of mixed effects logistic regression
Example 1: A researcher sampled applications to 40 different colleges to study factors that predict admittance into college. Predictors include student’s high school GPA, extracurricular activities, and SAT scores. Some colleges are more or less selective, so the baseline probability of admittance into each of the colleges is different. College-level predictors include whether the college is public or private, the current student-to-teacher ratio, and the college’s rank.
Example 2: A large HMO wants to know what patient and physician factors are most related to whether a patient’s lung cancer goes into remission after treatment as part of a larger study of treatment outcomes and quality of life in patients with lunge cancer.
Example 3: A television station wants to know how time and advertising campaigns affect whether people view a television show. They sample people from four cities for six months. Each month, they ask whether the people had watched a particular show or not in the past week. After three months, they introduced a new advertising campaign in two of the four cities and continued monitoring whether or not people had watched the show.
Description of the data
In this example, we are going to explore Example 2 about lung cancer using a simulated dataset, which we have posted online. A variety of outcomes were collected on patients, who are nested within doctors, who are in turn nested within hospitals. There are also a few doctor level variables, such as Experience that we will use in our example.
*Grab the most recent version from the internet
insheet using "https://stats.idre.ucla.edu/stat/data/hdp.csv", comma
foreach i of varlist familyhx smokinghx sex cancerstage school {
encode `i', gen(`i'2)
drop `i'
rename `i'2 `i'
}
Here is a general summary of the whole dataset.
summarize
Variable | Obs Mean Std. Dev. Min Max
-------------+--------------------------------------------------------
tumorsize | 8525 70.88067 12.06833 33.96859 116.4579
co2 | 8525 1.605207 .1238528 1.222361 2.128112
pain | 8525 5.473314 1.504302 1 9
wound | 8525 5.731848 1.525207 1 9
mobility | 8525 6.080469 1.484188 1 9
-------------+--------------------------------------------------------
ntumors | 8525 3.066276 2.550696 0 9
nmorphine | 8525 3.62393 2.503595 0 18
remission | 8525 .2957185 .4563918 0 1
lungcapacity | 8525 .7740865 .1756213 .0161208 .9997957
age | 8525 50.97205 6.275041 26.32264 74.48235
-------------+--------------------------------------------------------
married | 8525 .6 .4899267 0 1
lengthofstay | 8525 5.492199 1.04961 1 10
wbc | 8525 5997.58 995.1962 2131.302 9776.412
rbc | 8525 4.995027 .2891493 3.918932 6.06487
bmi | 8525 29.07269 6.648291 18.38268 58
-------------+--------------------------------------------------------
il6 | 8525 4.016984 2.858684 .0352107 23.72776
crp | 8525 4.973017 3.108535 .0451048 28.74212
did | 8525 203.3309 119.4691 1 407
experience | 8525 17.64129 4.075327 7 29
lawsuits | 8525 1.866393 1.486401 0 9
-------------+--------------------------------------------------------
hid | 8525 17.76422 10.21063 1 35
medicaid | 8525 .512513 .2072415 .1415814 .8187299
familyhx | 8525 1.2 .4000235 1 2
smokinghx | 8525 2.4 .8000469 1 3
sex | 8525 1.4 .4899267 1 2
-------------+--------------------------------------------------------
cancerstage | 8525 2.100059 .9436027 1 4
school | 8525 1.24868 .4322735 1 2
We can also get the frequencies for categorical or discrete variables, and the correlations for continuous predictors.
tab1 remission cancerstage lengthofstay
-> tabulation of remission remission | Freq. Percent Cum. ------------+----------------------------------- 0 | 6,004 70.43 70.43 1 | 2,521 29.57 100.00 ------------+----------------------------------- Total | 8,525 100.00 -> tabulation of cancerstage CancerStage | Freq. Percent Cum. ------------+----------------------------------- I | 2,558 30.01 30.01 II | 3,409 39.99 69.99 III | 1,705 20.00 89.99 IV | 853 10.01 100.00 ------------+----------------------------------- Total | 8,525 100.00 -> tabulation of lengthofstay LengthofSta | y | Freq. Percent Cum. ------------+----------------------------------- 1 | 2 0.02 0.02 2 | 14 0.16 0.19 3 | 181 2.12 2.31 4 | 1,196 14.03 16.34 5 | 2,896 33.97 50.31 6 | 2,874 33.71 84.02 7 | 1,168 13.70 97.72 8 | 183 2.15 99.87 9 | 10 0.12 99.99 10 | 1 0.01 100.00 ------------+----------------------------------- Total | 8,525 100.00
cor il6 crp lengthofstay experience
(obs=8,525) | il6 crp length~y experi~e -------------+------------------------------------ il6 | 1.0000 crp | 0.0024 1.0000 lengthofstay | -0.0066 0.0175 1.0000 experience | -0.0039 -0.0052 0.0128 1.0000
Analysis methods you might consider
Below is a list of analysis methods you may have considered.
- Mixed effects logistic regression, the focus of this page.
- Mixed effects probit regression is very similar to mixed effects logistic regression, but it uses the normal CDF instead of the logistic CDF. Both model binary outcomes and can include fixed and random effects.
- Fixed effects logistic regression is limited in this case because it may ignore necessary random effects and/or non independence in the data.
- Fixed effects probit regression is limited in this case because it may ignore necessary random effects and/or non independence in the data.
- Logistic regression with clustered standard errors. These can adjust for non independence but does not allow for random effects.
- Probit regression with clustered standard errors. These can adjust for non independence but does not allow for random effects.
Mixed effects logistic regression
Below we use the xtmelogit command to estimate a mixed effects logistic regression model with il6, crp, and lengthofstay as patient level continuous predictors, cancerstage as a patient level categorical predictor (I, II, III, or IV), experience as a doctor level continuous predictor, and a random intercept by did, doctor ID.
Estimating and interpreting generalized linear mixed models (GLMMs, of which mixed effects logistic regression is one) can be quite challenging. If you are just starting, we highly recommend reading this page first Introduction to GLMMs. It covers some of the background and theory as well as estimation options, inference, and pitfalls in more detail.
xtmelogit remission il6 crp i.cancerstage lengthofstay experience || did: , intpoints(10)
Refining starting values: Iteration 0: log likelihood = -3824.2872 Iteration 1: log likelihood = -3699.8676 Iteration 2: log likelihood = -3690.0759 Performing gradient-based optimization: Iteration 0: log likelihood = -3690.0759 Iteration 1: log likelihood = -3689.6396 Iteration 2: log likelihood = -3689.6382 Iteration 3: log likelihood = -3689.6382 Mixed-effects logistic regression Number of obs = 8,525 Group variable: did Number of groups = 407 Obs per group: min = 2 avg = 20.9 max = 40 Integration points = 10 Wald chi2(7) = 394.80 Log likelihood = -3689.6382 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ remission | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- il6 | -.0567716 .0115181 -4.93 0.000 -.0793466 -.0341966 crp | -.0214829 .0102182 -2.10 0.036 -.0415101 -.0014557 | cancerstage | II | -.4139322 .0757618 -5.46 0.000 -.5624225 -.2654419 III | -1.003461 .0982811 -10.21 0.000 -1.196089 -.8108341 IV | -2.337028 .1580562 -14.79 0.000 -2.646813 -2.027244 | lengthofstay | -.1211816 .0336326 -3.60 0.000 -.1871002 -.055263 experience | .1200848 .0274496 4.37 0.000 .0662847 .173885 _cons | -2.052611 .5313988 -3.86 0.000 -3.094133 -1.011088 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ did: Identity | sd(_cons) | 2.014511 .102283 1.823692 2.225297 ------------------------------------------------------------------------------ LR test vs. logistic model: chibar2(01) = 2435.28 Prob >= chibar2 = 0.0000
The first part gives us the iteration history, tells us the type of model, total number of observations, number of groups, and the grouping variable. Stata also indicates that the estimates are based on 10 integration points and gives us the log likelihood as well as the overall Wald chi square test that all the fixed effects parameters (excluding the intercept) are simultaneously zero. We used 10 integration points (how this works is discussed in more detail here). As we use more integration points, the approximation becomes more accurate converging to the ML estimates; however, more points are more computationally demanding and can be extremely slow or even intractable with today’s technology.
The next section is a table of the fixed effects estimates. For many applications, these are what people are primarily interested in. The estimates represent the regression coefficients. These are unstandardized and are on the logit scale. The estimates are followed by their standard errors (SEs). As is common in GLMs, the SEs are obtained by inverting the observed information matrix (negative second derivative matrix). However, for GLMMs, this is again an approximation. The approximations of the coefficient estimates likely stabilize faster than do those for the SEs. Thus if you are using fewer integration points, the estimates may be reasonable, but the approximation of the SEs may be less accurate. The Wald tests, \(\frac{Estimate}{SE}\), rely on asymptotic theory, here referring to as the highest level unit size converges to infinity, these tests will be normally distributed, and from that, p values (the probability of obtaining the observed estimate or more extreme, given the true estimate is 0). Using the same assumptions, approximate 95% confidence intervals are calculated.
The last section gives us the random effect estimates. This represents the estimated standard deviation in the intercept on the logit scale. Had there been other random effects, such as random slopes, they would also appear here.
If we wanted odds ratios instead of coefficients on the logit scale, we could exponentiate the estimates and CIs. We can do this in Stata by using the OR option. Note that we do not need to refit the model. Note that the random effects parameter estimates do not change. This is not the standard deviation around the exponentiated constant estimate, it is still for the logit scale.
xtmelogit, or
Mixed-effects logistic regression Number of obs = 8,525 Group variable: did Number of groups = 407 Obs per group: min = 2 avg = 20.9 max = 40 Integration points = 10 Wald chi2(7) = 394.80 Log likelihood = -3689.6382 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ remission | Odds Ratio Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- il6 | .9448098 .0108824 -4.93 0.000 .9237197 .9663815 crp | .9787462 .010001 -2.10 0.036 .9593397 .9985454 | cancerstage | II | .6610458 .050082 -5.46 0.000 .569827 .766867 III | .3666082 .0360307 -10.21 0.000 .3023745 .4444872 IV | .0966143 .0152705 -14.79 0.000 .0708768 .131698 | lengthofstay | .8858731 .0297942 -3.60 0.000 .8293606 .9462363 experience | 1.127593 .0309519 4.37 0.000 1.068531 1.189919 _cons | .1283993 .0682312 -3.86 0.000 .0453143 .363823 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ did: Identity | sd(_cons) | 2.014511 .102283 1.823692 2.225297 ------------------------------------------------------------------------------ LR test vs. logistic model: chibar2(01) = 2435.28 Prob >= chibar2 = 0.0000
Multilevel bootstrapping
Inference from GLMMs is complicated. Except for cases where there are many observations at each level (particularly the highest), assuming that \(\frac{Estimate}{SE}\) is normally distributed may not be accurate. A variety of alternatives have been suggested including Monte Carlo simulation, Bayesian estimation, and bootstrapping. Each of these can be complex to implement. We are going to focus on a small bootstrapping example.
Bootstrapping is a resampling method. It is by no means perfect, but it is conceptually straightforward and easy to implement in code. One downside is that it is computationally demanding. For large datasets or complex models where each model takes minutes to run, estimating on thousands of bootstrap samples can easily take hours or days. In the example for this page, we use a very small number of samples, but in practice you would use many more. Perhaps 1,000 is a reasonable starting point.
For single level models, we can implement a simple random sample with replacement for bootstrapping. With multilevel data, we want to resample in the same way as the data generating mechanism. We start by resampling from the highest level, and then stepping down one level at a time. The Biostatistics Department at Vanderbilt has a nice page describing the idea here. Unfortunately, Stata does not have an easy way to do multilevel bootstrapping. However, it can do cluster bootstrapping fairly easily, so we will just do that. The cluster bootstrap is the data generating mechanism if and only if once the cluster variable is selected, all units within it are sampled. In our case, if once a doctor was selected, all of her or his patients were included. If instead, patients were sampled from within doctors, but not necessarily all patients for a particular doctor, then to truly replicate the data generation mechanism, we could write our own program to resample from each level at a time.
Below we use the bootstrap command, clustered by did, and ask for a new, unique ID variable to be generated called newdid. For the purpose of demonstration, we only run 20 replicates. In practice you would probably want to run several hundred or a few thousand. We set the random seed to make the results reproducible. Note for the model, we use the newly generated unique ID variable, newdid and for the sake of speed, only a single integration point. If you take this approach, it is probably best to use the observed estimates from the model with 10 integration points, but use the confidence intervals from the bootstrap, which can be obtained by calling estat bootstrap after the model.
set seed 10
bootstrap _b, cluster(did) idcluster(newdid) rep(2): xtmelogit remission il6 crp i.cancerstage lengthofstay experience || newdid: , intpoints(1)
(running xtmelogit on estimation sample) Bootstrap replications (2) ----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 .. Mixed-effects logistic regression Number of obs = 8,525 Group variable: newdid Number of groups = 407 Obs per group: min = 2 avg = 20.9 max = 40 Integration points = 1 Wald chi2(1) = . Log likelihood = -3695.9925 Prob > chi2 = . (Replications based on 407 clusters in did) ------------------------------------------------------------------------------ | Observed Bootstrap Normal-based remission | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- il6 | -.0567427 .0027501 -20.63 0.000 -.0621329 -.0513525 crp | -.0215193 .0049849 -4.32 0.000 -.0312895 -.0117491 | cancerstage | II | -.4143991 .0176847 -23.43 0.000 -.4490604 -.3797378 III | -1.002911 .0444254 -22.58 0.000 -1.089983 -.9158389 IV | -2.334538 .1448591 -16.12 0.000 -2.618457 -2.05062 | lengthofstay | -.120899 .0870575 -1.39 0.165 -.2915285 .0497306 experience | .1206723 .0245079 4.92 0.000 .0726377 .1687069 _cons | -2.067761 .0991897 -20.85 0.000 -2.262169 -1.873353 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ | Observed Bootstrap Normal-based Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ newdid: Identity | sd(_cons) | 1.981108 .0095011 1.962573 1.999818 ------------------------------------------------------------------------------ LR test vs. logistic model: chibar2(01) = 2422.57 Prob >= chibar2 = 0.0000 Note: Log-likelihood calculations are based on the Laplacian approximation.
estat bootstrap
Mixed-effects logistic regression Number of obs = 8,525 Replications = 2 (Replications based on 407 clusters in did) ------------------------------------------------------------------------------ | Observed Bootstrap remission | Coef. Bias Std. Err. [95% Conf. Interval] -------------+---------------------------------------------------------------- eq1 | il6 | -.05674271 .0002956 .00275013 -.0583917 -.0545025 (BC) crp | -.0215193 .0004014 .00498487 -.0246427 -.017593 (BC) 1b.cancers~e | 0 0 0 . . (BC) 2.cancerst~e | -.41439913 .0212375 .01768467 . . (BC) 3.cancerst~e | -1.0029111 .0298688 .04442539 -1.004456 -.9416288 (BC) 4.cancerst~e | -2.3345383 .0636574 .14485906 -2.373312 -2.16845 (BC) lengthofstay | -.12089897 .0006922 .08705751 -.1817657 -.0586478 (BC) experience | .12067228 .0125137 .0245079 .1158563 .1505157 (BC) _cons | -2.067761 -.2513205 .09918974 . . (BC) -------------+---------------------------------------------------------------- lns1_1_1 | _cons | .68365628 .010208 .00479585 . . (BC) ------------------------------------------------------------------------------ (BC) bias-corrected confidence interval
If you happen to have a multicore version of Stata, that will help with speed. If not, as long as you specify different random seeds, you can run each bootstrap in separate instances of Stata and combine the results. In particular, you can use the saving option to bootstrap to save the estimates from each bootstrap replicate and then combine the results. For example, suppose you ultimately wanted 1000 replicates, you could do 250 replicates on four different cores or machines, save the results, combine the data files, and then get the more stable confidence interval estimates from the greater number of replicates without it taking so long.
Predicted probabilities and graphing
These results are great to put in the table or in the text of a research manuscript; however, the numbers can be tricky to interpret. Visual presentations are helpful to ease interpretation and for posters and presentations. As models become more complex, there are many options. We will discuss some of them briefly and give an example how you could do one.
In a logistic model, the outcome is commonly on one of three scales:
- Log odds (also called logits), which is the linearized scale
- Odds ratios (exponentiated log odds), which are not on a linear scale
- Probabilities, which are also not on a linear scale
For tables, people often present the odds ratios. For visualization, the logit or probability scale is most common. There are some advantages and disadvantages to each. The logit scale is convenient because it is linearized, meaning that a 1 unit increase in a predictor results in a coefficient unit increase in the outcome and this holds regardless of the levels of the other predictors (setting aside interactions for the moment). A downside is the scale is not very interpretable. It is hard for readers to have an intuitive understanding of logits. Conversely, probabilities are a nice scale to intuitively understand the results; however, they are not linear. This means that a one unit increase in the predictor, does not equal a constant increase in the probability—the change in probability depends on the values chosen for the other predictors. In ordinary logistic regression, you could just hold all predictors constant, only varying your predictor of interest. However, in mixed effects logistic models, the random effects also bear on the results. Thus, if you hold everything constant, the change in probability of the outcome over different values of your predictor of interest are only true when all covariates are held constant and you are in the same group, or a group with the same random effect. The effects are conditional on other predictors and group membership, which is quite narrowing. An attractive alternative is to get the average marginal probability. That is, across all the groups in our sample (which is hopefully representative of your population of interest), graph the average change in probability of the outcome across the range of some predictor of interest.
We are going to explore an example with average marginal probabilities. These take more work than conditional probabilities, because you have to calculate separate conditional probabilities for every group and then average them. It is also not easy to get confidence intervals around these average marginal effects in a frequentist framework (although they are trivial to obtain from Bayesian estimation).
First, let’s define the general procedure using the notation from here. We create \(\mathbf{X}_{i}\) by taking \(\mathbf{X}\) and setting a particular predictor of interest, say in column \(j\), to a constant. If we only cared about one value of the predictor, \(i \in \{1\}\). However, more commonly, we want a range of values for the predictor in order to plot how the predicted probability varies across its range. We can do this by taking the observed range of the predictor and taking \(k\) samples evenly spaced within the range. For example, suppose our predictor ranged from 5 to 10, and we wanted 6 samples, \(\frac{10 – 5}{6 – 1} = 1\), so each sample would be 1 apart from the previous and they would be: \(\{5, 6, 7, 8, 9, 10\}\). Then we create \(k\) different \(\mathbf{X}_{i}\)s where \(i \in \{1, \ldots, k\}\) where in each case, the \(j\)th column is set to some constant. Then we calculate: $$ \boldsymbol{\eta}_{i} = \mathbf{X}_{i}\boldsymbol{\beta} + \mathbf{Z}\boldsymbol{\gamma} $$ These are all the different linear predictors. Finally, we take \(h(\boldsymbol{\eta})\), which gives us \(\boldsymbol{\mu}_{i}\), which are the conditional expectations on the original scale, in our case, probabilities. We can then take the expectation of each \(\boldsymbol{\mu}_{i}\) and plot that against the value our predictor of interest was held at. We could also make boxplots to show not only the average marginal predicted probability, but also the distribution of predicted probabilities.
You may have noticed that a lot of variability goes into those estimates. We are using \(\mathbf{X}\) only holding our predictor of interest at a constant, which allows all the other predictors to take on values in the original data. Also, we have left \(\mathbf{Z}\boldsymbol{\gamma}\) as in our sample, which means some groups are more or less represented than others. If we had wanted, we could have re-weighted all the groups to have equal weight. We chose to leave all these things as-is in this example based on the assumption that our sample is truly a good representative of our population of interest. Rather than attempt to pick meaningful values to hold covariates at (even the mean is not necessarily meaningful, particularly if a covariate as a bimodal distribution, it may be that no participant had a value at or near the mean), we used the values from our sample. This also suggests that if our sample was a good representation of the population, then the average marginal predicted probabilities are a good representation of the probability for a new random sample from our population.
Now that we have some background and theory, let’s see how we actually go about calculating these things. First we define a Mata function to do the calculations.
mata
void mypredict(string scalar predictor, string scalar re) {
vars = tokens(st_global("e(datasignaturevars)"))
sol = st_matrix("e(b)")
nf = st_numscalar("e(k_f)")
u = st_data(., re)
b = sol[1::nf]
b = (b[cols(b)], b[1::(cols(b) - 1)])'
dat = st_data(., vars)
for(index = 1; index <= cols(dat); index++) {
if (vars[index] == predictor) {
break
}
}
xvar = dat[., index]
xvals = rangen(min(xvar), max(xvar), 100)
n = rows(xvals)
X = dat[., 2::nf]
r = rows(X)
X = (J(r, 1, 1), X)
avgmarginal = J(n, 1, .)
for (i = 1; i <= n; i++) {
X[., index] = J(r, 1, xvals[i])
linpred = X * b
linpred = linpred + u
linkpred = (1 :/ (1 :+ exp(-linpred)))
avgmarginal[i] = mean(linkpred)
}
results = (xvals, avgmarginal)
st_matrix("e(predictedprobabilities)", results)
}
end
Now we just need to run our model, and then get the average marginal predicted probabilities for lengthofstay. The function mypredict does not work with factor variables, so we will dummy code cancer stage manually.
tabulate cancerstage, gen(cancerstage)
quietly xtmelogit remission il6 crp cancerstage2 cancerstage3 cancerstage4 lengthofstay experience || did: , intpoints(10)
predict re*, reffects
mata mypredict("lengthofstay", "re1")
matrix pr=e(predictedprobabilities)
svmat pr
twoway (line pr2 pr1), ytitle(Average Marginal Predicted Probabilities) xtitle(LengthofStay) ylabel(0(.25)1)
Actually, those predicted probabilities are incorrect. The note from predict indicated that missing values were generated. For this model, Stata seemed unable to provide accurate estimates of the conditional modes. See the R page for a correct example. Nevertheless, in your data, this is the procedure you would use in Stata, and assuming the conditional modes are estimated well, the process works.
Three level mixed effects logistic regression
We have looked at a two level logistic model with a random intercept in depth. This is the simplest mixed effects logistic model possible. Now we are going to briefly look at how you can add a third level and random slope effects as well as random intercepts.
Below we estimate a three level logistic model with a random intercept for doctors and a random intercept for hospitals. In this examples, doctors are nested within hospitals, meaning that each doctor belongs to one and only one hospital. The alternative case is sometimes called “cross classified” meaning that a doctor may belong to multiple hospitals, such as if some of the doctor’s patients are from hospital A and others from hospital B. Note that this model takes several minutes to run on our machines.
xtmelogit remission age lengthofstay i.familyhx il6 crp i.cancerstage experience || hid: || did: , intpoints(10)
Refining starting values: Iteration 0: log likelihood = -3705.8549 Iteration 1: log likelihood = -3627.7993 Iteration 2: log likelihood = -3592.7991 Performing gradient-based optimization: Iteration 0: log likelihood = -3592.7991 Iteration 1: log likelihood = -3583.7752 Iteration 2: log likelihood = -3582.5521 Iteration 3: log likelihood = -3580.7413 Iteration 4: log likelihood = -3580.7376 Iteration 5: log likelihood = -3580.7376 Mixed-effects logistic regression Number of obs = 8,525 ---------------------------------------------------------------------------- | No. of Observations per Group Integration Group Variable | Groups Minimum Average Maximum Points ----------------+----------------------------------------------------------- hid | 35 134 243.6 377 10 did | 407 2 20.9 40 10 ---------------------------------------------------------------------------- Wald chi2(9) = 533.89 Log likelihood = -3580.7376 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ remission | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | -.0160884 .0060669 -2.65 0.008 -.0279793 -.0041975 lengthofstay | -.0419453 .0364578 -1.15 0.250 -.1134013 .0295107 | familyhx | yes | -1.308861 .0955083 -13.70 0.000 -1.496054 -1.121668 il6 | -.0585618 .011737 -4.99 0.000 -.0815658 -.0355577 crp | -.0232041 .0103835 -2.23 0.025 -.0435555 -.0028528 | cancerstage | II | -.321855 .0785607 -4.10 0.000 -.4758313 -.1678788 III | -.8633261 .1027024 -8.41 0.000 -1.064619 -.6620331 IV | -2.162179 .1657534 -13.04 0.000 -2.48705 -1.837308 | experience | .1264035 .0277251 4.56 0.000 .0720633 .1807437 _cons | -1.622896 .5907345 -2.75 0.006 -2.780714 -.4650775 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ hid: Identity | sd(_cons) | .5269101 .1596348 .2909733 .9541572 -----------------------------+------------------------------------------------ did: Identity | sd(_cons) | 1.997495 .1051354 1.801706 2.21456 ------------------------------------------------------------------------------ LR test vs. logistic model: chi2(2) = 2493.70 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference.
We can easily add random slopes to the model as well, and allow them to vary at any level. We are just going to add a random slope for lengthofstay that varies between doctors. We use a single integration point for the sake of time.
xtmelogit remission age lengthofstay i.familyhx il6 crp i.cancerstage experience || hid: || did: lengthofstay, intpoints(1)
Refining starting values: Iteration 0: log likelihood = -3767.3719 (not concave) Iteration 1: log likelihood = -3640.7894 Iteration 2: log likelihood = -3592.6058 Performing gradient-based optimization: Iteration 0: log likelihood = -3592.6058 Iteration 1: log likelihood = -3569.788 (not concave) Iteration 2: log likelihood = -3568.005 Iteration 3: log likelihood = -3560.1972 Iteration 4: log likelihood = -3559.8573 Iteration 5: log likelihood = -3559.8383 Iteration 6: log likelihood = -3559.8372 Iteration 7: log likelihood = -3559.8372 Mixed-effects logistic regression Number of obs = 8,525 ---------------------------------------------------------------------------- | No. of Observations per Group Integration Group Variable | Groups Minimum Average Maximum Points ----------------+----------------------------------------------------------- hid | 35 134 243.6 377 1 did | 407 2 20.9 40 1 ---------------------------------------------------------------------------- Wald chi2(9) = 576.99 Log likelihood = -3559.8372 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ remission | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- age | -.0154194 .0060913 -2.53 0.011 -.0273582 -.0034807 lengthofstay | -.1927663 .0453468 -4.25 0.000 -.2816444 -.1038882 | familyhx | yes | -1.349591 .0970406 -13.91 0.000 -1.539787 -1.159395 il6 | -.0590629 .0117796 -5.01 0.000 -.0821505 -.0359752 crp | -.0214644 .0104164 -2.06 0.039 -.0418803 -.0010486 | cancerstage | II | -.2969842 .078352 -3.79 0.000 -.4505512 -.1434171 III | -.8686914 .1039774 -8.35 0.000 -1.072483 -.6648994 IV | -2.301839 .1721777 -13.37 0.000 -2.639301 -1.964377 | experience | .1060332 .0242676 4.37 0.000 .0584695 .1535969 _cons | -.5547783 .5585291 -0.99 0.321 -1.649475 .5399187 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ hid: Identity | sd(_cons) | .7300311 .1411391 .4997767 1.066367 -----------------------------+------------------------------------------------ did: Independent | sd(length~y) | .3715136 .0287975 .3191499 .4324689 sd(_cons) | .2712679 .8170012 .0007409 99.31988 ------------------------------------------------------------------------------ LR test vs. logistic model: chi2(3) = 2535.50 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference. Note: Log-likelihood calculations are based on the Laplacian approximation.
Things to consider
- Sample size: Often the limiting factor is the sample size at the highest unit of analysis. For example, having 500 patients from each of ten doctors would give you a reasonable total number of observations, but not enough to get stable estimates of doctor effects nor of the doctor-to-doctor variation. 10 patients from each of 500 doctors (leading to the same total number of observations) would be preferable.
- Parameter estimation: Because there are not closed form solutions for GLMMs, you must use some approximation. Three are fairly common.
- Quasi-likelihood approaches use a Taylor series expansion to approximate the likelihood. Thus parameters are estimated to maximize the quasi-likelihood. That is, they are not true maximum likelihood estimates. A Taylor series uses a finite set of differentiations of a function to approximate the function, and power rule integration can be performed with Taylor series. With each additional term used, the approximation error decreases (at the limit, the Taylor series will equal the function), but the complexity of the Taylor polynomial also increases. Early quasi-likelihood methods tended to use a first order expansion, more recently a second order expansion is more common. In general, quasi-likelihood approaches are the fastest (although they can still be quite complex), which makes them useful for exploratory purposes and for large datasets. Because of the bias associated with them, quasi-likelihoods are not preferred for final models or statistical inference.
- The true likelihood can also be approximated using numerical integration. Quadrature methods are common, and perhaps most common among these use the Gaussian quadrature rule, frequently with the Gauss-Hermite weighting function. It is also common to incorporate adaptive algorithms that adaptively vary the step size near points with high error. The accuracy increases as the number of integration points increases. Using a single integration point is equivalent to the so-called Laplace approximation. Each additional integration point will increase the number of computations and thus the speed to convergence, although it increases the accuracy. Adaptive Gauss-Hermite quadrature might sound very appealing and is in many ways. However, the number of function evaluations required grows exponentially as the number of dimensions increases. A random intercept is one dimension, adding a random slope would be two. For three level models with random intercepts and slopes, it is easy to create problems that are intractable with Gaussian quadrature. Consequently, it is a useful method when a high degree of accuracy is desired but performs poorly in high dimensional spaces, for large datasets, or if speed is a concern.
- A final set of methods particularly useful for multidimensional integrals are Monte Carlo methods including the famous Metropolis-Hastings algorithm and Gibbs sampling which are types of Markov chain Monte Carlo (MCMC) algorithms. Although Monte Carlo integration can be used in classical statistics, it is more common to see this approach used in Bayesian statistics.
- Complete or quasi-complete separation: Complete separation means that the outcome variable separate a predictor variable completely, leading perfect prediction by the predictor variable. Particularly if the outcome is skewed, there can also be problems with the random effects. For example, if one doctor only had a few patients and all of them either were in remission or were not, there will be no variability within that doctor.
See also
References
- Agresti, A. (2013). Categorical Data Analysis (3rd Ed.), Hoboken, New Jersey: John Wiley & Sons.