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

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

Looking at the data

head(Fungi)
str(Fungi)
#View(Fungi) ## Try that too!

Preparing the data

FungiAlive <- subset(Fungi, alive == TRUE)

or

library(dplyr)
Fungi %>%
  filter(alive) -> FungiAlive

Silly attempt

coef(mod_silly <- lm(growth ~ species* (T36 + T38 + PT36 + PT38), data = FungiAlive))

Looking at the data structure to figure out the issue (ignoring species for now)

FungiAlive %>%
  count(T36, T38, PT36, PT38) ## base R function table() is not great with many dimensions

Creating the pre-treatment variable

FungiAlive$PT <- "no"
FungiAlive$PT[FungiAlive$PT36 == 1] <- "36"
FungiAlive$PT[FungiAlive$PT38 == 1] <- "38"
FungiAlive$PT <- factor(FungiAlive$PT)
table(FungiAlive$PT, FungiAlive$species)

Creating the treatment variable

FungiAlive$T <- "no"
FungiAlive$T[FungiAlive$T36 == 1] <- "36"
FungiAlive$T[FungiAlive$T38 == 1] <- "38"
FungiAlive$T <- factor(FungiAlive$T)
table(FungiAlive$T, FungiAlive$species)

The revised dataset

FungiAlive

Fitting the model and predicting the effects: solution 1

mod <- lm(growth ~ species * (PT + T), data = FungiAlive)
newdata <- data.frame(PT = c("36", "38", "no", "no", "36", "38", "no", "no"),
                       T = c("36", "38", "36", "38", "36", "38", "36", "38"),
                       species = c(rep("M_A", 4), rep("M_B", 4)))
newdata$pred <- predict(mod, newdata = newdata)
newdata

plot(newdata[newdata$species == "M_A", "pred"],
     type = "b", axes = FALSE, col = "red",
     ylim = range(newdata$pred),
     ylab = "Growth rate", xlab = "Experimental condition")
points(newdata[newdata$species == "M_B", "pred"],
       type = "b", col = "blue")
axis(side = 2, las = 1)
axis(side = 1, at = 1:4,
     labels = c("36_36", "38_38", "no_36", "no_38"))
box()
legend("top", pch = c(1, 1), col = c("blue", "red"), horiz = TRUE,
       legend = c("M_A", "M_B"), box.lty = 2, title = "Species")

Fitting the model and predicting the effects: solution 2

FungiAlive$condition <-
  factor(paste(FungiAlive$species, FungiAlive$PT, FungiAlive$T, sep = "_"))
  table(FungiAlive$condition)

You can in fact directly get the estimates for each treatment by simply removing the intercept:

mod2 <- lm(growth ~ -1 + condition, data = FungiAlive)
data.frame(coef(mod2))
newdata2 <- data.frame(condition = levels(FungiAlive$condition))
newdata2$pred <- predict(mod2, newdata = newdata2)

plot(c(
  newdata2[newdata2$condition == "M_A_36_36", "pred"],
  newdata2[newdata2$condition == "M_A_38_38", "pred"],
  newdata2[newdata2$condition == "M_A_no_36", "pred"],
  newdata2[newdata2$condition == "M_A_no_38", "pred"]),
  type = "b", axes = FALSE, col = "red",
  ylim = range(newdata2$pred),
  ylab = "Growth rate", xlab = "Experimental condition")
points(c(
  newdata2[newdata2$condition == "M_B_36_36", "pred"],
  newdata2[newdata2$condition == "M_B_38_38", "pred"],
  newdata2[newdata2$condition == "M_B_no_36", "pred"],
  newdata2[newdata2$condition == "M_B_no_38", "pred"]),
  type = "b", col = "blue")
axis(side = 2, las = 1)
axis(side = 1, at = 1:4,
     labels = c("36_36", "38_38", "no_36", "no_38"))
box()
legend("top", pch = c(1, 1), col = c("blue", "red"), horiz = TRUE,
       legend = c("M_A", "M_B"), box.lty = 2, title = "Species")

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