In stevencarlislewalker/lme4ord: Generalized Linear Mixed Models with Flexible (Co)Variance Structures

library(lme4ord)
library(splines)
library(numDeriv)
library(ape)
library(minqa)
library(ggplot2)

We simulate some data.

set.seed(6)
n <- 100 # 100
m <- 500 # 500
dl <- dims_to_vars(data.list(y = 1 * (matrix(rnorm(n * m), n, m) > 0),
                             x = rnorm(n), z = rnorm(m),
                             dimids = c("sites", "species")))
df <- as.data.frame(dl)
phy <- rtree(n = m)
phy <- compute.brlen(phy, method = "Grafen", power = 0.1)
Vphy <- stanCov(vcv(phy))
dimnames(Vphy) <- rep(list(1:m), 2)
covList <- list(species = Vphy)
form <- y ~ x * z + (1 | species) # + (0 + x | species)
parsedForm <- glmercFormula(form, df, covList = covList)
covarSim <- c(0.5) # c(0.5, -0.2, 0.5)
parsedForm <- within(parsedForm, Lambdat@x[] <- mapToCovFact(covarSim))
X <- model.matrix(nobars(form), df) # fixed effects design matrix
Z <- t(parsedForm$Lambdat %*% parsedForm$Zt) # random effects design
                                             # matrix with
                                             # phylogenetic
                                             # covariances
fixefSim <- rnorm(ncol(X)) # fixed effects
u <- rnorm(ncol(Z)) # whitened random effects
p <- plogis(as.numeric(X %*% fixefSim + Z %*% u)) # probability of observation
dl$y <- rbinom(nrow(df), 1, p) # presence-absence data
df <- as.data.frame(dl) # reconstruct the data frame with new
                        # structured response

Fit the generating model to the simulated data.

mod <- glmerc(form, df, covList = covList,
              optControl = list(iprint = 4L, maxfun = 500))

Look at a slice over the covariance parameter through the optimum.

sfun <- mkSfun(mod, 1)
pp <- seq(0.4, 0.62, length = 100)
fn <- function(par) {
    dv <- sfun(par)
    c(dv, environment(mod$dfun)$pp$ldL2(),
      environment(mod$dfun)$pp$sqrL(1),
      environment(mod$dfun)$resp$aic(),
      2 * log(det(environment(mod$dfun)$pp$L())))
}
ss <- sapply(pp, fn)
par(mfrow = c(5, 1), mar = c(4, 4, 0, 0))
plot(pp, ss[1,], type = "l", ylab = "deviance slice")
plot(pp, apply(ss[-1, ], 2, sum), type = "l", ylab = "deviance slice")
plot(pp, ss[2,], type = "l", ylab = "log determinant")
plot(pp, ss[3,], type = "l", ylab = "squared length of u")
plot(pp, ss[4,], type = "l", ylab = "conditional deviance")

2 * sum(log(environment(mod$dfun)$pp$L()@x))
rho <- environment(mod$dfun)
rho$resp$sqrtrwt * tcrossprod(rho$pp$Lambdat %*% rho$pp$Zt)@x
sum(log( + 1))

That's a tough curve!

Compute the Wald confidence intervals two ways.

centDfun <- function(par) mod$dfun(par) - mod$opt$fval
getWaldCI <- function(Hess, mle) {
    intervalAroundZero <- outer(sqrt(diag(solve(0.5*Hess))), 1.96 * c(-1, 1))
    intervalAroundEst <- sweep(intervalAroundZero, 1, mle, "+")
    colnames(intervalAroundEst) <- c("lower", "upper")
    intervalAroundEst
}
hh <- lme4:::deriv12(centDfun, pars(mod))$Hessian
ci <- getWaldCI(hh, pars(mod))

params <- as.data.frame(cbind(estimated = pars(mod),
                              true = c(covar = covarSim, fixef = fixefSim), ci))
params$parameter <- rownames(params)
ggplot(params) +
    geom_point(aes(estimated, parameter)) +
    geom_point(aes(true, parameter), col = "red") +
    geom_segment(aes(lower, parameter, xend = upper, yend = parameter)) +
    theme_bw()

Compute the deviance profile for a particular parameter.

pfun <- mkPfun(mod, whichPar <- 1)
true <- c(covarSim, fixefSim)[whichPar]
mle <- pars(mod)[whichPar]
expFac <- 100 * apply(ci, 1, diff)
grdMin <- -expFac[whichPar] + ci[whichPar, 1]
grdMax <-  expFac[whichPar] + ci[whichPar, 2]
grd <- seq(min(max(mod$lower[whichPar], grdMin), true),
           max(grdMax, true), length = 10)
pro <- sapply(grd, pfun)

Plot the results

par(mar = c(4, 4, 4, 1))
plot(interpSpline(grd, pro),
     lwd = 3, las = 1,
     xlab = names(pars(mod))[whichPar],
     ylab = "Profiled deviance")
axis(3, c(true, mle), c("true", "mle"), las = 2)
abline(v = true, col = "red")
abline(v = mle, col = "blue")
abline(h = c(0, qchisq(0.95, 1)), lty = 2)
abline(v = ci[whichPar, ], lty = 2)

stevencarlislewalker/lme4ord documentation built on May 30, 2019, 4:43 p.m.