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*.

library(ggplot2) library(GGally) library(reshape2) library(lme4)

library(compiler) library(parallel) library(boot)

**Example 1:**A researcher sampled applications to 40 different colleges to study factor that predict admittance into college. Predictors include student's high school GPA, extracurricular activities, and SAT scores. Some schools are more or less selective, so the baseline probability of admittance into each of the schools is different. School level predictors include whether the school is public or private, the current student-to-teacher ratio, and the school'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 lunger 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.

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.

hdp <- read.csv("http://www.karlin.mff.cuni.cz/~pesta/prednasky/NMFM404/Data/hdp.csv") hdp <- within(hdp, { Married <- factor(Married, levels = 0:1, labels = c("no", "yes")) DID <- factor(DID) HID <- factor(HID) })

ggpairs(hdp[, c("IL6", "CRP", "LengthofStay", "Experience")])

`CancerStage`

. Because `LengthofStay`

is coded discretely in days, we can examine how `CancerStage`

is associated with it using bubble plots. The area of each bubble is proportional to the number of observations with those values. For the continuous predictors, we use violin plots with jittered data values. All of the raw data is presented separated by `CancerStage`

. To alleviate overplotting and see the values better, we add a small amount of random noise (primarily to the x axis) as well as set the alpha transparency. Although the jittered dots are helpful for seeing the raw data, it can be difficult to get a precise sense of the distribution. For that, we add violin plots. Violin plots are just kernel density plots reflected around the plotting axis. We plot the violin plots on top of the jittered points with a transparency so that you can stil see the raw data, but the violin plots are dominant. Because both `IL6`

and `CRP`

tend to have skewed distributions, we use a square root scale on the y axis. The distributions look fairly normal and symmetric, although you can still see the long right tail, even using a square root scale (note that only the scale was shifted, the values themselves are not transformed, which is important because this lets you see and interpret the actual scores, rather than the square root of the scores).
ggplot(hdp, aes(x = CancerStage, y = LengthofStay)) + stat_sum(aes(size = ..n.., group = 1)) + scale_size_area(max_size=10)

tmp <- melt(hdp[, c("CancerStage", "IL6", "CRP")], id.vars="CancerStage") ggplot(tmp, aes(x = CancerStage, y = value)) + geom_jitter(alpha = .1) + geom_violin(alpha = .75) + facet_grid(variable ~ .) + scale_y_sqrt()

tmp <- melt(hdp[, c("remission", "IL6", "CRP", "LengthofStay", "Experience")], id.vars="remission") ggplot(tmp, aes(factor(remission), y = value, fill=factor(remission))) + geom_boxplot() + facet_wrap(~variable, scales="free_y")

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.

Below we use the `glmer`

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.

# estimate the model and store results in m m <- glmer(remission ~ IL6 + CRP + CancerStage + LengthofStay + Experience + (1 | DID), data = hdp, family = binomial, nAGQ = 10)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = ## control$checkConv, : Model failed to converge with max|grad| = 0.00107021 ## (tol = 0.001, component 2)

# print the mod results without correlations among fixed effects print(m, corr = FALSE)

## Generalized linear mixed model fit by maximum likelihood (Adaptive ## Gauss-Hermite Quadrature, nAGQ = 10) [glmerMod] ## Family: binomial ( logit ) ## Formula: ## remission ~ IL6 + CRP + CancerStage + LengthofStay + Experience + ## (1 | DID) ## Data: hdp ## AIC BIC logLik deviance df.resid ## 7397.276 7460.733 -3689.638 7379.276 8516 ## Random effects: ## Groups Name Std.Dev. ## DID (Intercept) 2.015 ## Number of obs: 8525, groups: DID, 407 ## Fixed Effects: ## (Intercept) IL6 CRP CancerStageII ## -2.05164 -0.05677 -0.02148 -0.41398 ## CancerStageIII CancerStageIV LengthofStay Experience ## -1.00340 -2.33682 -0.12119 0.12003

The next section gives us basic information that can be used to compare models, followed by the random effect estimates. This represents the estimated variability in the intercept on the logit scale. Had there been other random effects, such as random slopes, they would also appear here. The top section concludes with the total number of observations, and the number of level 2 observations. In our case, this includes the total number of patients (8,525) and doctors (407).

The last 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, 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).

It can be nice to get confidence intervals (CIs). We can get rough estimates using the SEs.

se <- sqrt(diag(vcov(m))) # table of estimates with 95% CI (tab <- cbind(Est = fixef(m), LL = fixef(m) - 1.96 * se, UL = fixef(m) + 1.96 * se))

## Est LL UL ## (Intercept) -2.05164301 -3.09327975 -1.010006262 ## IL6 -0.05677289 -0.07934881 -0.034196973 ## CRP -0.02148351 -0.04151146 -0.001455563 ## CancerStageII -0.41397971 -0.56247659 -0.265482839 ## CancerStageIII -1.00340497 -1.19603829 -0.810771638 ## CancerStageIV -2.33681678 -2.64659735 -2.027036207 ## LengthofStay -0.12118784 -0.18710918 -0.055266504 ## Experience 0.12002708 0.06622201 0.173832158

exp(tab)

## Est LL UL ## (Intercept) 0.12852357 0.04535296 0.3642167 ## IL6 0.94480862 0.92371767 0.9663811 ## CRP 0.97874562 0.95933834 0.9985455 ## CancerStageII 0.66101436 0.56979617 0.7668356 ## CancerStageIII 0.36662895 0.30238982 0.4445149 ## CancerStageIV 0.09663476 0.07089202 0.1317253 ## LengthofStay 0.88586754 0.82935318 0.9462329 ## Experience 1.12752739 1.06846390 1.1898558

Inference from GLMMs is complicated. Except for cases where there are many observations at each level (particularly the highest), assuming that EstimateSE 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. In our case, we first will sample from doctors, and then within each doctor sampled, we will sample from their patients. To do this, we first need to write a function to resample at each level. The Biostatistics Department at Vanderbilt has a nice page describing the idea http://biostat.mc.vanderbilt.edu/wiki/Main/HowToBootstrapCorrelatedData.

sampler <- function(dat, clustervar, replace = TRUE, reps = 1) { cid <- unique(dat[, clustervar[1]]) ncid <- length(cid) recid <- sample(cid, size = ncid * reps, replace = TRUE) if (replace) { rid <- lapply(seq_along(recid), function(i) { cbind(NewID = i, RowID = sample(which(dat[, clustervar] == recid[i]), replace = TRUE)) }) } else { rid <- lapply(seq_along(recid), function(i) { cbind(NewID = i, RowID = which(dat[, clustervar] == recid[i])) }) } dat <- as.data.frame(do.call(rbind, rid)) dat$Replicate <- factor(cut(dat$NewID, breaks = c(1, ncid * 1:reps), include.lowest = TRUE, labels = FALSE)) dat$NewID <- factor(dat$NewID) return(dat) }

set.seed(20) tmp <- sampler(hdp, "DID", reps = 100) bigdata <- cbind(tmp, hdp[tmp$RowID, ])

`lme4`

package on the cluster. Finally, we write a function to fit the model and return the estimates. The call to `glmer()`

is wrapped in `try`

because not all models may converge on the resampled data. This catches the error and returns it, rather than stopping processing.
f <- fixef(m) r <- getME(m, "theta") cl <- makeCluster(4) clusterExport(cl, c("bigdata", "f", "r")) clusterEvalQ(cl, require(lme4))

## [[1]] ## [1] TRUE ## ## [[2]] ## [1] TRUE ## ## [[3]] ## [1] TRUE ## ## [[4]] ## [1] TRUE

myboot <- function(i) { object <- try(glmer(remission ~ IL6 + CRP + CancerStage + LengthofStay + Experience + (1 | NewID), data = bigdata, subset = Replicate == i, family = binomial, nAGQ = 1, start = list(fixef = f, theta = r)), silent = TRUE) if (class(object) == "try-error") return(object) c(fixef(object), getME(object, "theta")) }

`parLapplyLB`

function, which loops through every replicate, giving them out to each node of the cluster to estimate the models. The "LB" stands for load balancing, which means replicates are distributed as a node completes its current job. This is valuable because not all replicates will converge, and if there is an error and it happens early on, one node may be ready for a new job faster than another node. There is some extra communication overhead, but this is small compared to the time it takes to fit each model. The results from all nodes are aggregated back into a single list, stored in the object `res`

. Once that is done, we can shut down the local cluster, which terminates the additional `R`

instances and frees memory.
start <- proc.time() res <- parLapplyLB(cl, X = levels(bigdata$Replicate), fun = myboot) end <- proc.time() # shut down the cluster stopCluster(cl)

# calculate proportion of models that successfully converged success <- sapply(res, is.numeric) mean(success)

## [1] 1

# combine successful results bigres <- do.call(cbind, res[success]) # calculate 2.5th and 97.5th percentiles for 95% CI (ci <- t(apply(bigres, 1, quantile, probs = c(0.025, 0.975))))

## 2.5% 97.5% ## (Intercept) -3.63944293 -1.007201621 ## IL6 -0.08743904 -0.030617950 ## CRP -0.05018041 0.003039692 ## CancerStageII -0.59388626 -0.242818869 ## CancerStageIII -1.30051865 -0.753841874 ## CancerStageIV -2.92950112 -2.029114134 ## LengthofStay -0.22031001 -0.044977250 ## Experience 0.06794416 0.207263575 ## NewID.(Intercept) 2.04852634 2.461902813

# All results finaltable <- cbind(Est = c(f, r), SE = c(se, NA), BootMean = rowMeans(bigres), ci) # round and print round(finaltable, 3)

## Est SE BootMean 2.5% 97.5% ## (Intercept) -2.052 0.531 -2.211 -3.639 -1.007 ## IL6 -0.057 0.012 -0.059 -0.087 -0.031 ## CRP -0.021 0.010 -0.023 -0.050 0.003 ## CancerStageII -0.414 0.076 -0.417 -0.594 -0.243 ## CancerStageIII -1.003 0.098 -1.043 -1.301 -0.754 ## CancerStageIV -2.337 0.158 -2.458 -2.930 -2.029 ## LengthofStay -0.121 0.034 -0.143 -0.220 -0.045 ## Experience 0.120 0.027 0.129 0.068 0.207 ## DID.(Intercept) 2.015 NA 2.260 2.049 2.462

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).

We create \({\bf X}_i\) by taking \({\bf 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, (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 \({\bf X}_i\) where \(i\in\{1,\ldots,k\}\) where in each case, the jth column is set to some constant. Then we calculate: \[ {\boldsymbol\eta}_i={\bf X}_i{\boldsymbol\beta}+{\bf Z}_i{\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 \({\bf 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 \({\bf 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. We get a summary of `LengthofStay`

, our predictor of interest, and then get 100 values across its range to use in prediction. We make a copy of our data so we can fix the values of one of the predictors and then use the `predict`

function to calculate the predicted values. All random effects are included by default, see `?predict.merMod`

for more details. *Note that the predict method for mixed effects models is new and currently is only in the development version of lme4, so make sure that you have that installed.*

# temporary data tmpdat <- hdp[, c("IL6", "CRP", "CancerStage", "LengthofStay", "Experience", "DID")] summary(hdp$LengthofStay)

## Min. 1st Qu. Median Mean 3rd Qu. Max. ## 1.000 5.000 5.000 5.492 6.000 10.000

jvalues <- with(hdp, seq(from = min(LengthofStay), to = max(LengthofStay), length.out = 100)) # calculate predicted probabilities and store in a list pp <- lapply(jvalues, function(j) { tmpdat$LengthofStay <- j predict(m, newdata = tmpdat, type = "response") })

# average marginal predicted probability across a few different Lengths of # Stay sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)

## [1] 0.3652354 0.3366378 0.3075494 0.2796343 0.2530078 0.2277651

# get the means with lower and upper quartiles plotdat <- t(sapply(pp, function(x) { c(M = mean(x), quantile(x, c(0.25, 0.75))) })) # add in LengthofStay values and convert to data frame plotdat <- as.data.frame(cbind(plotdat, jvalues)) # better names and show the first few rows colnames(plotdat) <- c("PredictedProbability", "Lower", "Upper", "LengthofStay") head(plotdat)

## PredictedProbability Lower Upper LengthofStay ## 1 0.3652354 0.08489961 0.6155674 1.000000 ## 2 0.3637090 0.08404758 0.6129569 1.090909 ## 3 0.3621849 0.08320332 0.6103400 1.181818 ## 4 0.3606630 0.08236679 0.6077167 1.272727 ## 5 0.3591434 0.08153791 0.6050872 1.363636 ## 6 0.3576261 0.08071664 0.6024515 1.454545

# plot average marginal predicted probabilities ggplot(plotdat, aes(x = LengthofStay, y = PredictedProbability)) + geom_line() + ylim(c(0, 1))

ggplot(plotdat, aes(x = LengthofStay, y = PredictedProbability)) + geom_linerange(aes(ymin = Lower, ymax = Upper)) + geom_line(size = 2) + ylim(c(0, 1))

`LengthofStay`

we could do it for each level of `CancerStage`

.
# calculate predicted probabilities and store in a list biprobs <- lapply(levels(hdp$CancerStage), function(stage) { tmpdat$CancerStage[] <- stage lapply(jvalues, function(j) { tmpdat$LengthofStay <- j predict(m, newdata = tmpdat, type = "response") }) }) # get means and quartiles for all jvalues for each level of CancerStage plotdat2 <- lapply(biprobs, function(X) { temp <- t(sapply(X, function(x) { c(M=mean(x), quantile(x, c(.25, .75))) })) temp <- as.data.frame(cbind(temp, jvalues)) colnames(temp) <- c("PredictedProbability", "Lower", "Upper", "LengthofStay") return(temp) }) # collapse to one data frame plotdat2 <- do.call(rbind, plotdat2) # add cancer stage plotdat2$CancerStage <- factor(rep(levels(hdp$CancerStage), each = length(jvalues))) # show first few rows head(plotdat2)

## PredictedProbability Lower Upper LengthofStay CancerStage ## 1 0.4474697 0.1547464 0.7328397 1.000000 I ## 2 0.4458035 0.1533108 0.7306772 1.090909 I ## 3 0.4441385 0.1518862 0.7285036 1.181818 I ## 4 0.4424747 0.1504724 0.7263191 1.272727 I ## 5 0.4408123 0.1490695 0.7241237 1.363636 I ## 6 0.4391511 0.1476774 0.7219174 1.454545 I

# graph it ggplot(plotdat2, aes(x = LengthofStay, y = PredictedProbability)) + geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = CancerStage), alpha = .15) + geom_line(aes(colour = CancerStage), size = 2) + ylim(c(0, 1)) + facet_wrap(~ CancerStage)

ggplot(data.frame(Probs = biprobs[[4]][[100]]), aes(Probs)) + geom_histogram() + scale_x_sqrt(breaks = c(0.01, 0.1, 0.25, 0.5, 0.75))

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. In `glmer`

you do not need to specify whether the groups are nested or cross classified, `R`

can figure it out based on the data. We use the same `(1 | ID)`

general syntax to indicate the intercept (1) varying by some ID. For models with more than a single scalar random effect, `glmer`

only supports a single integration point, so we use `nAGQ=1`

.

# estimate the model and store results in m m3a <- glmer(remission ~ Age + LengthofStay + FamilyHx + IL6 + CRP + CancerStage + Experience + (1 | DID) + (1 | HID), data = hdp, family = binomial, nAGQ=1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = ## control$checkConv, : Model failed to converge with max|grad| = 0.0640339 ## (tol = 0.001, component 3)

# print the mod results without correlations among fixed effects print(m3a, corr=FALSE)

## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation) [glmerMod] ## Family: binomial ( logit ) ## Formula: ## remission ~ Age + LengthofStay + FamilyHx + IL6 + CRP + CancerStage + ## Experience + (1 | DID) + (1 | HID) ## Data: hdp ## AIC BIC logLik deviance df.resid ## 7199.094 7283.703 -3587.547 7175.094 8513 ## Random effects: ## Groups Name Std.Dev. ## DID (Intercept) 1.9523 ## HID (Intercept) 0.5487 ## Number of obs: 8525, groups: DID, 407; HID, 35 ## Fixed Effects: ## (Intercept) Age LengthofStay FamilyHxyes ## -1.68828 -0.01493 -0.04467 -1.30662 ## IL6 CRP CancerStageII CancerStageIII ## -0.05686 -0.02216 -0.31846 -0.85692 ## CancerStageIV Experience ## -2.13751 0.12686

It can also be useful to look at the distribution of the conditional modes, which we do with caterpillar polots below. The blue dots are the conditional models with error bars. We do this for both doctors and hospitals. For example for doctors, we can see a bit of a long right tail in that there are more extreme positive than negative values. For the doctors, we suppress their IDs (using the `scales=list(y = list(alternating=0))`

argument) because there are so many, but we leave them in for the hospitals.

dotplot(ranef(m3a, which = "DID", postVar = TRUE), scales = list(y = list(alternating = 0)))

## Warning in ranef.merMod(m3a, which = "DID", postVar = TRUE): 'postVar' is ## deprecated: please use 'condVar' instead

## $DID

dotplot(ranef(m3a, which = "HID", postVar = TRUE))

## Warning in ranef.merMod(m3a, which = "HID", postVar = TRUE): 'postVar' is ## deprecated: please use 'condVar' instead

## $HID

`LengthofStay`

that varies between doctors. As in regular `R`

formulae, we use the + operator to "add" an effect, and we do it in the section for doctor random effects. All terms in one group of parentheses use an unstructured covariance matrix, you can get a diagonal covariance structure by splitting the grouping into separate pieces. Between groupings is assumed indepedent.
# estimate the model and store results in m m3b <- glmer(remission ~ Age + LengthofStay + FamilyHx + IL6 + CRP + CancerStage + Experience + (1 + LengthofStay | DID) + (1 | HID), data = hdp, family = binomial, nAGQ = 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = ## control$checkConv, : Model failed to converge with max|grad| = 0.120807 ## (tol = 0.001, component 1)

# print the mod results without correlations among fixed effects print(m3b, corr = FALSE)

## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation) [glmerMod] ## Family: binomial ( logit ) ## Formula: ## remission ~ Age + LengthofStay + FamilyHx + IL6 + CRP + CancerStage + ## Experience + (1 + LengthofStay | DID) + (1 | HID) ## Data: hdp ## AIC BIC logLik deviance df.resid ## 7147.761 7246.471 -3559.880 7119.761 8511 ## Random effects: ## Groups Name Std.Dev. Corr ## DID (Intercept) 0.5056 ## LengthofStay 0.3715 -0.11 ## HID (Intercept) 0.7301 ## Number of obs: 8525, groups: DID, 407; HID, 35 ## Fixed Effects: ## (Intercept) Age LengthofStay FamilyHxyes ## -0.54702 -0.01527 -0.18974 -1.34010 ## IL6 CRP CancerStageII CancerStageIII ## -0.05864 -0.02103 -0.29421 -0.86536 ## CancerStageIV Experience ## -2.29502 0.10458

dotplot(ranef(m3b, which = "DID", postVar = TRUE), scales = list(y = list(alternating = 0)))

## Warning in ranef.merMod(m3b, which = "DID", postVar = TRUE): 'postVar' is ## deprecated: please use 'condVar' instead

## $DID

dotplot(ranef(m3b, which = "HID", postVar = TRUE), scales = list(y = list(alternating = 0)))

## Warning in ranef.merMod(m3b, which = "HID", postVar = TRUE): 'postVar' is ## deprecated: please use 'condVar' instead

## $HID

*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.

- UCLA: IDRE (Institute for Digital Research and Education). Data Analysis Examples. from http://www.ats.ucla.edu/stat/dae/ (accessed January 31, 2014)
- Agresti, A. (2013). Categorical Data Analysis (3rd Ed.), Hoboken, New Jersey: John Wiley & Sons.