How to plot a binary mixed effect model for visual presentation

Question

I am trying to plot the result of a binary mixed effect model for visual representation in a paper.

I use lme to fit the mixed model:

M2 <- lme(Pass ~ zone.time + length + Fat,
      random =~ 1 | Year)

Pass = binary 1/0 zone.time, length & Fat = continuous

yielding:

Linear mixed-effects model fit by maximum likelihood
Data: DF1 
   AIC      BIC    logLik
39.05604 47.25981 -13.52802

Random effects:
 Formula: ~1 | Year
        (Intercept)  Residual
StdDev: 5.03879e-06 0.3857927

Fixed effects: Pass ~ zone.time + length + Fat 
                Value Std.Error DF   t-value p-value
(Intercept)  4.549716 1.2384118 24  3.673832  0.0012
zone.time    0.299438 0.1239111 24  2.416559  0.0236
length      -0.006718 0.0019492 24 -3.446603  0.0021
Fat         -0.051460 0.0213211 24 -2.413563  0.0238
 Correlation: 
          (Intr) zon.tm length
zone.time  0.045              
length    -0.979 -0.168       
Fat       -0.447 -0.191  0.330

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.9097237 -0.7802111 -0.1410353  0.5683329  2.0908188 

Number of Observations: 29
Number of Groups: 2  

I then go about calculating the predicted values and standard errors:

MyData <- expand.grid(zone.time    = seq(1,3.6, length = 10),
                  length = seq(525, 740, length = 10),
                  Fat = seq(3.7, 17, length = 10))
X <- model.matrix(~zone.time + length + Fat, data = MyData)

Extract parameters and parameter covariance matrix

betas    <- fixef(M2)

for sample data use

betas<- structure(c(4.54971638246632, 0.299438350935228, -0.00671801197327911,-0.0514597408192487), .Names = c("(Intercept)", "zone.time", "length","Fat"))

.

Covbetas <- vcov(M2)

for sample data use:

Covbetas <- structure(c(1.32212400759181, 0.0059001955657893, -0.00203725210229123, 
-0.0101822039057957, 0.0059001955657893, 0.0132361635192455, 
-3.50672281561515e-05, -0.000434188193496185, -0.00203725210229123, 
-3.50672281561515e-05, 3.27522409259271e-06, 1.18250356798504e-05, 
-0.0101822039057957, -0.000434188193496185, 1.18250356798504e-05, 
0.000391886154502855), .Dim = c(4L, 4L), .Dimnames = list(c("(Intercept)", 
"zone.time", "length", "Fat"), c("(Intercept)", "zone.time", 
"length", "Fat")))

Calculate the fitted values in the predictor scale

MyData$eta <- X %*% betas
MyData$Pi  <- exp(MyData$eta) / (1 + exp(MyData$eta))

Calculate the SEs on the scale of the predictor function

MyData$se <- sqrt(diag(X %*% Covbetas %*% t(X)))
MyData$SeUp  <- exp(MyData$eta + 1.96 *MyData$se) / (1 + exp(MyData$eta  + 1.96 *MyData$se))
MyData$SeLo  <- exp(MyData$eta - 1.96 *MyData$se) / (1 + exp(MyData$eta  - 1.96 *MyData$se))

head(MyData)

Is this the correct method for calculating the predicted values?

How do I go about plotting this for visual presentation?

should i be using something like

library(effects)
plot(allEffects(M2, default.levels=50))

or ggplot2

Answer 1

Kind of confusing, but I am gathering you want to extract the fitted results from your mixed-effects model and then plot them. Is that correct?

Create similar data

set.seed(64)

fooDF <- data.frame(Pass = rbinom(n = 100, size = 1, prob = 0.5), zone.time = rnorm(n = 100), length = rnorm(n = 100),
                    Fat = rnorm(n = 100), Year = seq(1913, 2012))

M2 <- lme(Pass ~ zone.time + length + Fat,
              random =~ 1 | Year, data = fooDF)

You can get overall predicted outcome by

    head(fitted(M2, level = 0))

     1913      1914      1915      1916      1917      1918 
    0.4948605 0.7506069 0.5317316 0.5429997 0.6584630 0.7555496 

You could simply plot the fits like this

plot(fitted(M2, level = 0))

You can also use a variable in your dataset, say Fat, on the x-axis with the fitted values on on the y.

plotDF <- data.frame(fat = fooDF$Fat, fitted = fitted(M2, level = 0))
plot(plotDF)

plot(fitted(M2))

As you can see, with this fictitious data the relationship is linear.