In courtiol/LM2GLMM: Data science for biologists: Generalized linear modelling with R

library(LM2GLMM)
library(car)
knitr::opts_chunk$set(cache = TRUE, fig.align = "center", fig.width = 6, fig.height = 6,
cache.path = "./cache_knitr/Exo_GLM_solution/", fig.path = "./fig_knitr/Exo_GLM_solution/")

Data

set.seed(1L)
Aliens <- simulate_Aliens_GLM()
fit_binar <- glm(happy ~ humans_eaten, family = binomial(), data = Aliens)
fit_binom <- glm(cbind(blue_eyes, pink_eyes) ~ humans_eaten, family = binomial(),  data = Aliens)
fit_poiss <- glm(eggs  ~ humans_eaten, family = poisson(),  data = Aliens)

Goal

• plot the predictions and confidence intervals for fit_binar, fit_binom and fit_poiss for a number of humans eaten varying between 0 and 15.

Solution for fit_binar

Without packages

data.for.pred <- data.frame(humans_eaten = 0:15)
p <- predict(fit_binar, newdata = data.for.pred, se.fit = TRUE)
p$upr <- p$fit + qnorm(0.975) * p$se.fit
p$lwr <- p$fit + qnorm(0.025) * p$se.fit
p.mu <- data.frame(fit = binomial()$linkinv(p$fit),
                   upr = binomial()$linkinv(p$upr),
                   lwr = binomial()$linkinv(p$lwr))

plot(p.mu$fit ~ data.for.pred$humans_eaten, type = "l", lwd = 2, las = 1, ylim = c(0, 1),
     ylab = "Probability of being happy", xlab = "Number of humans eaten")
points(p.mu$upr ~  data.for.pred$humans_eaten, type = "l", lwd = 2, lty = 2)
points(p.mu$lwr ~  data.for.pred$humans_eaten, type = "l", lwd = 2, lty = 2)
rug(x = jitter(Aliens$humans_eaten[Aliens$happy == 1]), col = "red", side = 3)
rug(x = jitter(Aliens$humans_eaten[Aliens$happy == 0]), col = "red", side = 1)

Using {spaMM}

library(spaMM)
fit_binar_spaMM <- fitme(happy ~ humans_eaten, family = binomial(), data = Aliens)
p <- predict(fit_binar_spaMM, newdata = data.for.pred, binding = "pred")
p <- cbind(p, get_intervals(fit_binar_spaMM, newdata = data.for.pred, intervals = "fixefVar"))
plot(pred ~ humans_eaten, type = "l", lwd = 2, las = 1, ylim = c(0, 1),
     ylab = "Probability of being happy", xlab = "Number of humans eaten", data = p)
points(fixefVar_0.025 ~  humans_eaten, type = "l", lwd = 2, lty = 2, data = p)
points(fixefVar_0.975 ~  humans_eaten, type = "l", lwd = 2, lty = 2, data = p)
rug(x = jitter(Aliens$humans_eaten[Aliens$happy == 1]), col = "red", side = 3)
rug(x = jitter(Aliens$humans_eaten[Aliens$happy == 0]), col = "red", side = 1)

Solution for fit_binom

Without packages

data.for.pred <- data.frame(humans_eaten = 0:15)
p <- predict(fit_binom, newdata = data.for.pred, se.fit = TRUE)
p$upr <- p$fit + qnorm(0.975) * p$se.fit
p$lwr <- p$fit + qnorm(0.025) * p$se.fit
p.mu <- data.frame(fit = binomial()$linkinv(p$fit),
                   upr = binomial()$linkinv(p$upr),
                   lwr = binomial()$linkinv(p$lwr))

plot(p.mu$fit ~ data.for.pred$humans_eaten, type = "l", lwd = 2, las = 1, ylim = c(0, 1),
     ylab = "Proportion of blue eyes", xlab = "Number of humans eaten")
points(p.mu$upr ~  data.for.pred$humans_eaten, type = "l", lwd = 2, lty = 2)
points(p.mu$lwr ~  data.for.pred$humans_eaten, type = "l", lwd = 2, lty = 2)
points(blue_eyes/(blue_eyes + pink_eyes) ~ jitter(humans_eaten), data = Aliens)

Using {spaMM}

fit_binom_spaMM <- fitme(cbind(blue_eyes, pink_eyes) ~ humans_eaten,
                         family = binomial(), data = Aliens)
p <- predict(fit_binom_spaMM, newdata = data.for.pred, binding = "pred")
p <- cbind(p, get_intervals(fit_binom_spaMM, newdata = data.for.pred, intervals = "fixefVar"))
plot(pred ~ humans_eaten, type = "l", lwd = 2, las = 1, ylim = c(0, 1),
     ylab = "Proportion of blue eyes", xlab = "Number of humans eaten", data = p)
points(fixefVar_0.025 ~  humans_eaten, type = "l", lwd = 2, lty = 2, data = p)
points(fixefVar_0.975 ~  humans_eaten, type = "l", lwd = 2, lty = 2, data = p)
points(blue_eyes/(blue_eyes + pink_eyes) ~ jitter(humans_eaten), data = Aliens)

Solution for fit_poiss

Without packages

p <- predict(fit_poiss, newdata = data.for.pred, se.fit = TRUE)
p$upr <- p$fit + qnorm(0.975) * p$se.fit
p$lwr <- p$fit + qnorm(0.025) * p$se.fit
p.mu <- data.frame(fit = poisson()$linkinv(p$fit),
                   upr = poisson()$linkinv(p$upr),
                   lwr = poisson()$linkinv(p$lwr))
plot(p.mu$fit ~ data.for.pred$humans_eaten, type = "l", lwd = 2, las = 1,
     ylim = range(Aliens$eggs),
     ylab = "Number of eggs", xlab = "Number of humans eaten")
points(p.mu$upr ~  data.for.pred$humans_eaten, type = "l", lwd = 2, lty = 2)
points(p.mu$lwr ~  data.for.pred$humans_eaten, type = "l", lwd = 2, lty = 2)
points(eggs ~ jitter(humans_eaten), data = Aliens)

Using {spaMM}

fit_poiss_spaMM <- fitme(eggs  ~ humans_eaten, family = poisson(), data = Aliens)
p <- predict(fit_poiss_spaMM, newdata = data.for.pred, binding = "pred")
p <- cbind(p, get_intervals(fit_poiss_spaMM, newdata = data.for.pred, intervals = "fixefVar"))
plot(pred ~ humans_eaten, type = "l", lwd = 2, las = 1, ylim = range(Aliens$eggs),
     ylab = "Number of eggs", xlab = "Number of humans eaten", data = p)
points(fixefVar_0.025 ~  humans_eaten, type = "l", lwd = 2, lty = 2, data = p)
points(fixefVar_0.975 ~  humans_eaten, type = "l", lwd = 2, lty = 2, data = p)
points(eggs ~ jitter(humans_eaten), data = Aliens)

Using {effects} or {visreg}

library(effects)
plot(allEffects(fit_poiss))

library(visreg)
visreg(fit_poiss, scale = "response")

courtiol/LM2GLMM documentation built on July 3, 2022, 7:42 a.m.