Implications – lsmeans

Understanding a complex model with high-order interactions is tough. The best tool is plotting, but quickly visualizing such models is not easy with flexible plotting libraries.

Fortunately, lsmeans, although fairly inflexible in general, has plotting capabilities designed just for this purpose. Here’s a demo with the full (4-way interacting) model

## look at consequences
if(require(lsmeans)){
    ## here for 4-way interaction model (m2) -- mostly to see how can quickly plot response
    ## at factor levels in complex models
    ## relatively easy to plot mean prediction on response scale account for complex interaction
    lsmip(m2, nutrient_hetero ~  water_add | water_hetero * nutrient_add  , type='response')
    lsmip(m2, water_add ~ nutrient_add  | water_hetero * nutrient_hetero, type='response')
    lsmip(m2, water_add ~ water_hetero  | nutrient_add * nutrient_hetero, type='response')
}

## Loading required package: lsmeans

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'lsmeans'

if(require(lsmeans)){
    ## can also plot quickly mean estimates plus (Wald) confidence intervals
    ## not that it's possible to do finite size correction in lsmeans (but only for linear models)
    ## steps take a min w/ lots of RE for some reason
    m2.lsm <- lsmeans(m2, ~ water_add * water_hetero   | nutrient_hetero  * nutrient_add)
    plot(m2.lsm, layout =c(2, 3), horizontal=TRUE)
}

## Loading required package: lsmeans

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'lsmeans'

    ## for the simpler
    lsmip(m1, nutrient_hetero ~  water_add | water_hetero * nutrient_add  , type='response')
    lsmip(m1, water_add ~ nutrient_add  | water_hetero * nutrient_hetero, type='response')
    lsmip(m1, water_add ~ water_hetero  | nutrient_add * nutrient_hetero, type='response')

    ## overall pattern doesn't differ too much from m1
    m1.lsm <- lsmeans(m1, ~ water_add * water_hetero   | nutrient_hetero  * nutrient_add)
    plot(m1.lsm, layout =c(2, 3), horizontal=TRUE)

quick & dirty

## OK -- back to the full model
## if residuals looking OK -- fastest way to get a quick look at coefficient uncertainty
## I'll make it a function for later use

quick_CI_df <- function(model, drop_intercept = TRUE,...) {
  ## coef table using WALD CIs for a merMod
  ci_dat <- confint(model, method='Wald', ...)
  ci_dat <- cbind(ci_dat, mean=rowMeans(ci_dat))
  ci_df <- data.frame(coef=row.names(ci_dat), ci_dat)
  names(ci_df)[2:3] <- c('lwr', 'upr')
  coef_colon <-   gregexpr("\\:", ci_df$coef)
  coef_order <- sapply(coef_colon, function(x) {
    if(x[1] == -1)
      return(0)
    else
      return(length(x))
  })
  ci_df$coef <- reorder(ci_df$coef, coef_order)
  if (drop_intercept)
    return(subset(ci_df, coef != '(Intercept)'))
  ci_df
}

m1_wald_ci <- quick_CI_df(m2)

lattice::dotplot(coef ~ mean, m1_wald_ci,
                 panel = function(x, y) {
                   panel.xyplot(x, y, pch=19)
                   panel.segments(m1_wald_ci$lwr, y, m1_wald_ci$upr, y)
                   panel.abline(v=0, lty=2)
                 })

Confidence Intervals – rolling your own

For more control than the simple plot above, you could use any of a variety of packages, e.g. coefplot2, arm::coefplot, here we just use builtin lme4::confint to build a dataframe.

I make a dataframe with potentially several different types of confidence intervals: Wald (assume asymptotic multivariate normality of likelihood), profile, and bootstrap. For the random effects, the output of bootstrap and profile methods have different names – as you can see in the plot below. The Wald method only estimates CIs for fixed effects.

compare_CI_df <- function(modprof, model, BOOT=FALSE, ...) {
  ci <- confint(model, method='Wald')
  cip <- confint(modprof)
  if(BOOT) {
      cib <- confint(model, method="boot", ...)
      cib_dat <- data.frame(cib, parameter=row.names(cib), type='Boot')
      names(cib_dat)[1:2] <- colnames(cib)
  }
  ci_dat <- data.frame(ci, parameter=row.names(ci), type='Wald')
  names(ci_dat)[1:2] <- colnames(ci)
  cip_dat <- data.frame(cip, parameter=row.names(cip), type='Prof')
  names(cip_dat)[1:2] <- colnames(cip)
  if (!BOOT) {
      ci_all_dat <- rbind(ci_dat, cip_dat)
  } else {
      ci_all_dat <- rbind(ci_dat, cip_dat, cib_dat)
  }
  fe_dat <- data.frame(parameter=names(fixef(model)), value=fixef(model))
  list(mean=fe_dat, bounds=ci_all_dat)
}

if(require(ggplot2)) {
    m1_ci <- compare_CI_df(m1.prall, m1, BOOT=TRUE, nsim=100)
    m1_mean <- m1_ci[['mean']]
    m1_bnd <- m1_ci[['bounds']]

    ggplot(m1_bnd[m1_bnd$parameter != '(Intercept)', ]) +
    geom_linerange(aes(parameter, ymax=`97.5 %`, ymin=`2.5 %`, color=type),
                   position=position_dodge(w=.5), size=1.5) +
    geom_hline(yintercept=0, linetype=2) +
    coord_flip() + theme_minimal() +
    geom_point(aes(parameter, value), size=3, shape=3,
               data=m1_mean[m1_mean$parameter != '(Intercept)', ])

} else {
    lattice::dotplot(parameter ~ c(`2.5 %`, `97.5 %`) | type , m1_bnd, abline=c(v=0, lty=2))
    }

## Loading required package: ggplot2

## Computing bootstrap confidence intervals ...

## 
## 38 message(s): boundary (singular) fit: see ?isSingular
## 1 warning(s): Model failed to converge with max|grad| = 0.00445866 (tol = 0.002, component 1)

## Warning: Removed 2 rows containing missing values (geom_linerange).

##   strong effect of nutrient homogeneity
##   effect of adding lots of water
## nutrient has some effect

Pulling data out of models with broom or ggplot2

Method ggplot2::fortify adds fitted and resid to data frame. New-ish CRAN package broom generalizes this (and more), really aiming provide an interface for ‘tidy’ (i.e., easy to plot)representations of model objects from many R functions (including packages and base R). The broom method augment really generalized version of fortify and is recommended instead

if(require(ggplot2)) {
    d_fort <- ggplot2::fortify(m1, d)
    ## adds .fortify, scresid, .resid
    ##head(d_fort)
}

if(require(broom)) {
    glm1 <- glance(m1) # model summary
    tdm1 <- tidy(m1, effects='fixed') #model coefficient summary
    augm1 <- augment(m1, d, se.fit=TRUE, type.predict='response') # model data with fits
    ## head(glm1)
    ## head(tdm1)
    ## head(augm1)
}

## Loading required package: broom

Unfortunately, the predict and confint methods for merMod aren’t yet integrated to broom or ggplot (or my code above for confidence intervals could be much cleaner at the slight cost of an additional package dependency).

GLMM: Lizard mate choice

Mechanisms of speciation – cool question!

• mate compatibility trails
• western North American skinks (Plestiodon skiltonianus spp)
• Richmond, J. Q., E. L. Jockusch, and A. M. Latimer 2010. Mechanical reproductive isolation facilitates parallel speciation in western North American scincid lizards. American Naturalist

Repeated measure of ind crossed with multiple trials, binary outcome

d2 <- read.delim("http://datadryad.org/bitstream/handle/10255/dryad.33579/Richmond%20et%20al%202011%20Datafile.txt?sequence=1")
names(d2) <- gsub('\\.', '', names(d2))

# rescale geo dist, sizediff
d2 <- within(d2, {
    geodist_scl <- (GeoDist - mean(GeoDist)) / var(GeoDist)
    gendist_scl <- (GenDist - mean(GenDist))
    sizediff_scl2 <- ((SizeDiff - mean(SizeDiff)) / var(SizeDiff))^2
    Series <- factor(Series)
    Female <- factor(Female)
    Male <- factor(Male)
  })

## cop - copulation
system.time(gm0 <- glmer(cop ~   sizediff_scl2   + (1 | Series) + (1 | Male)  + (1 | Female) ,
               data=d2, family='binomial'))

## boundary (singular) fit: see ?isSingular

##    user  system elapsed 
##    0.41    0.00    0.41

system.time(gm1 <- glmer(cop ~ gendist_scl + geodist_scl +  sizediff_scl2 + (1 | Series) +
                         (1 | Male)  + (1 | Female) , data=d2, family='binomial'))

## boundary (singular) fit: see ?isSingular

##    user  system elapsed 
##    1.03    0.00    1.03

anova(gm1, gm0) # no need for non-size predictors

## Data: d2
## Models:
## gm0: cop ~ sizediff_scl2 + (1 | Series) + (1 | Male) + (1 | Female)
## gm1: cop ~ gendist_scl + geodist_scl + sizediff_scl2 + (1 | Series) + 
## gm1:     (1 | Male) + (1 | Female)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## gm0  5 159.34 176.47 -74.673   149.34                         
## gm1  7 162.92 186.90 -74.460   148.92 0.4251      2     0.8085

plot(gm0, type=c("p"), id = 0.05 )

## non useful visually but with id = 0.05 gives sense of outlier
## see ?plot.merMod for explanation of test run

## it's more useful to plot residuals in GLM logistic combined across grouping factors
## or binned fitted values
plot(gm0,  factor(Trial) ~ resid(., type='pearson'),  abline=c(v=0), lty=2)

plot(gm0,  factor(Series) ~ resid(., type='pearson'),  abline=c(v=0), lty=2)

## better than default eh

## second plot suggests some issue with series, indicating we may want to
## fit a model accounting for differences between trial series
## (in fact, Richmond et al did fit such a model)
## you can see the outliers here too with
## plot(gm0, resid(., type='pearson') ~ fitted(.) | factor(Series), id = 0.05)

## main thing on pearson or deviance residuals is they're mostly less than 2
## (as assessed in first plot above -- second is the sqrt
## we could later drop these outlier and see if they change

if(require(gridExtra)) {
    do.call(gridExtra::grid.arrange, c(lattice::qqmath(ranef(gm0, condVar=TRUE)), list(nrow=1)))
} else {
    lattice::qqmath(ranef(gm0, condVar=TRUE))
  }

## Loading required package: gridExtra

## Warning: package 'gridExtra' was built under R version 3.6.1

## useful if you have a lot of RE to draw attention to strong ones
## maybe latter == dot = sample size?

## no reason to retain RE on Male, but also drop Female
## no real justification, just very slow profiling for many RE
gm0noindre <- glmer(cop ~   sizediff_scl2   + (1 | Series), data=d2, family='binomial')

system.time(#
    gm0.prall <- profile(gm0noindre, devtol=1e-6) # see https://stat.ethz.ch/pipermail/r-sig-mixed-models/2014q3/022394.html
)

##    user  system elapsed 
##    5.66    0.00    5.67

#system.time(
#    gm0_full.prall <- profile(gm0, devtol=1e-6)
#)
#   user  system elapsed
#344.092   0.539 344.613

lattice::xyplot(logProf(gm0.prall))

lattice::densityplot(gm0.prall)

lattice::splom(gm0.prall)

quick predictions

Again, the quickest way to visualize predictions is with lsmeans. Two things are different here:

• to plot GLM predictions on a meaningful scale, you need to pass type = 'response' to the plot function.
• for a quantitative predictor, the default will plot a single point at the mean of the predictor, to see prediction across the range, pass a list to the at argument
• for more meaningful plots, we should label the plot y-axis with unscaled size differences (but I didn’t do this)
• note that you can produce plots at various levels of interacting quantitative predictors using this technique

if(require(lsmeans))

gm0_lsm <-    lsmeans(gm0, ~ sizediff_scl2,
                  at = list(sizediff_scl2 = seq(from=0,
                              to = 1.5 * median(d2$sizediff_scl2), length.out = 10)))

## Loading required package: lsmeans

## Warning: package 'lsmeans' was built under R version 3.6.1

## Loading required package: emmeans

## Warning: package 'emmeans' was built under R version 3.6.1

## The 'lsmeans' package is now basically a front end for 'emmeans'.
## Users are encouraged to switch the rest of the way.
## See help('transition') for more information, including how to
## convert old 'lsmeans' objects and scripts to work with 'emmeans'.

## Warning in vcov.merMod(object, correlation = FALSE): variance-covariance matrix computed from finite-difference Hessian is
## not positive definite or contains NA values: falling back to var-cov estimated from RX

plot(gm0_lsm, type = 'response', ylab = 'scaled size diff', xlab = 'probability of copulation')

## customized predictions

The quick plots using lsmeans are great, but they assume Wald confidence intervals and don’t account for random effects. To deal with the former, we need to bootstrap, which I won’t cover here, but the problem of random effects is easily examined using the re.form argument of predict.merMod in lme4. Here I look at predictions accounting for all,, none, or some of random effects fitted:

## predicted effect
d2$pred_overall_gen <- predict(gm0, type='response', re.form=NA)
d2$pred_overall <- predict(gm0, type='response')
d2$pred_ind <- predict(gm0, type='response', re.form=~ (1|Female))
d2$pred_trial <- predict(gm0, type='response', re.form=~(1|Series))

## NB this changes row order from the fit!!!
d2_pred <- merge(d2, as.data.frame(table(Series=d2$Series)))

if(require(ggplot2)) {

g <- ggplot(d2_pred) +
     geom_point(aes(SizeDiff, cop),
                position=position_jitter(h=0.025, w=0)) +
     theme_minimal() +
     ylab("Probability copulation")


g1 <- g + geom_line(aes(SizeDiff, pred_trial, color=Series))

g2 <- g + geom_line(aes(SizeDiff, pred_overall_gen), color='orange', size = 1.5) +
    geom_point(aes(SizeDiff, pred_ind, alpha=as.numeric(Female)), color = 'blue') +
    theme(legend.position = "none")

gridExtra::grid.arrange(g1, g2, ncol=2)

} else {
    plot(cop ~ SizeDiff, d2_pred, ylab="Probability copulation")
    points(pred_ind ~ SizeDiff, d2_pred, pch=19, cex=.5, col='darkgray')
    points(pred_overall_gen ~ SizeDiff, d2_pred, pch=19, cex=.5)
    }

The plot on the left above shows predictions for various series of experimental trials, and the raw data. The plot on the right shows the overall fit (orange line – the grand mean with no random effects influencing the prediction), the predictions accounting just for random effects at the level of individual females (blue dots), and the raw data.

## could plot them all together -- code below gives some insight ## into the
## ggplot2 fortify method
if (require(ggplot2)) {
  g +  geom_point(aes(SizeDiff, plogis(.fitted)), color='darkgray', size=3,
                  data=fortify(gm0, d2))

## nb -- same as predict with all RE!
  with(fortify(gm0, d2),
       unique(pred_overall - plogis(.fitted)))
}

Checking model adequacy

With `lme4, you can also easily produce something akin to posterior-predictive checks, a Bayesian tool for model criticism.

Using simulate on an lme4 produces a set of response data implied by the model. By examining summary statistics of this dataset, and comparing them to the actual data, you can examine how well, or poorly, the model fits. Note that these checks are less rigorous than something like cross-validation. See ?simulate.merMod for additional information.

The upshot: these are relatively quick, and can provide a warning that your model isn’t capturing some aspect of the data. But, be aware there is controversy about their use for anything other than seeing how the model is capturing the data. Often, we’re interested in making more formal, frequentist inference. Though these checks are valuable, the p-like-values implied by them do not satisfy assumptions necessary for this kind of inference (i.e., don’t use them for model comparison or parameter inference).

## without RE
set.seed(1141)
gm0.sim_wo <- simulate(gm0, nsim=1000, re.form=NA)
sims <- sapply(gm0.sim_wo, function(x)  x)
ones <- colSums(sims)
zeros <- -colSums(sims-1)
perc <- ones/nrow(d2)
cvar <- apply(sims, 2, var)

## with RE
gm0.sim <- simulate(gm0, 1000 )
sims <- sapply(gm0.sim, function(x)  x)
re_ones <- colSums(sims)
re_zeros <- -colSums(sims-1)
re_perc <- re_ones/nrow(d2)
re_cvar <- apply(sims, 2, var)

pp_comp <- function(dist, cmp_fn, data, ...) {
  hist(dist, col='orange', xlab='value', ...)
  abline(v=cmp_fn(data), lwd=4, col='blue', lty=2)
  abline(v=mean(dist), lwd=4)
}

op <- par(mfrow=c(2,4))
  pp_comp(ones, function(x) sum(x$cop),
          d2, main='ones - no RE', xlim=c(45, 110))
  pp_comp(zeros, function(x) -sum(x$cop - 1),
          d2, main='zeros -no RE', xlim=c(120, 185))
  pp_comp(perc, function(x) sum(x$cop)/nrow(d2),
          d2, main='% cop - no RE', xlim=c(0, .8))
  pp_comp(cvar, function(x) var(x$cop),
          d2, main='var - no RE')

pp_comp(re_ones, function(x) sum(x$cop),
          d2, main='ones - all RE', xlim=c(45, 110))
  pp_comp(re_zeros, function(x) -sum(x$cop - 1),
          d2, main='zeros - all RE', xlim=c(120, 185))
  pp_comp(re_perc, function(x) sum(x$cop)/nrow(d2),
          d2, main='% cop - RE', xlim=c(0, .8))
  pp_comp(re_cvar, function(x) var(x$cop),
          d2, main='var - all RE')

par(op)

In this case I’ve plotted the simulated distributions of the number of ones, zeros, the proportion of copulations, and the variance in observed copulations – all with and without accounting for all the random effects. For each of these, I show the distribution of the statistic, its mean (black line), and the mean of the same statistic in the observed data.

All test quantities have broader distributions when we account for the random effects.

glmmADMB version (not run)

You can run many of the same graphical procedures (including the quick plots from lsmeans!) on output from other frequentist fitting packages. Try it!

## look at glmmadmb too -- not evaluated
gm1admb <- glmmadmb(cop ~ geodist_scl +  sizediff_scl2, random= ~1 | Series + 1 | Female , data=d2, family='binomial')

d2 <- data.frame(d2, predict(gm1admb, interval='confidence', type='link'))

plot(cop ~ SizeDiff, d2, ylim = c(0, 1))
lines(plogis(fit) ~ SizeDiff, d2, lwd=2)
lines(plogis(lwr) ~ SizeDiff, d2, lty=2, lwd=2)
lines(plogis(upr) ~ SizeDiff, d2, lty=2, lwd=2)
      geom_line(aes(SizeDiff, plogis(lwr)), color='orange', size=1, linetype=2) +
    geom_line(aes(SizeDiff, plogis(upr)), color='orange', size=1, linetype=2) +