works with glmmadmb nbinom1

R/compute_variances.R

@@ -697,6 +697,7 @@
       beta = ,
       ordbeta = stats::plogis(mu),
       poisson = ,
+      nbinom = ,
       nbinom1 = ,
       nbinom2 = ,
       `negative binomial` = exp(mu + 0.5 * as.vector(revar_null)),
@@ -726,6 +727,7 @@
 
         # (zero-inflated) negative binomial ----
         # --------------------------------------
+        nbinom = ,
         nbinom1 = ,
         nbinom2 = ,
         `negative binomial` = sig,
@@ -756,6 +758,7 @@
       switch(faminfo$family,
         nbinom2 = ,
         `negative binomial` = (1 / mu) + (1 / sig),
+        nbinom = ,
         poisson = ,
         nbinom1 = vv / mu,
         vv / mu^2

R/get_sigma.R

@@ -87,7 +87,11 @@ get_sigma <- function(x, ci = NULL, verbose = TRUE) {
 .get_sigma.glmmadmb <- function(x, ...) {
   check_if_installed("lme4")
   vc <- lme4::VarCorr(x)
-  attr(vc, "sc")
+  out <- attr(vc, "sc")
+  # sanity check
+  if (is.null(out)) {
+    out <- .safe(x$alpha)
+  }
 }
 
 

WIP/nakagawa_suppl.R

@@ -0,0 +1,141 @@
+# installing glmmADMB
+library(R2admb)
+library(glmmADMB)
+library(lme4)
+library(MASS)
+library(performance)
+
+
+# 12 different populations n = 960
+Population <- gl(12, 80, 960)
+# 120 containers (8 individuals in each container)
+Container <- gl(120, 8, 960)
+# Sex of the individuals. Uni-sex within each container (individuals are
+# sorted at the pupa stage)
+Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
+# Habitat at the collection site: dry or wet soil (four indiviudal from each
+# Habitat in each container)
+Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
+# Food treatment at the larval stage: special food ('Exp') or standard food
+# ('Cont')
+Treatment <- factor(rep(c("Cont", "Exp"), 480))
+# Data combined in a data frame
+Data <- data.frame(
+  Population = Population, Container = Container, Sex = Sex,
+  Habitat = Habitat, Treatment = Treatment
+)
+
+# Subset the design matrix (only females lay eggs)
+DataFemale <- Data[Data$Sex == "Female", ]
+# set seed for reproduciblity (this will enable one to get the same data
+# every time)
+set.seed(777)
+# simulation of the underlying random effects (Population and Container with
+# variance of 0.4 and 0.05, respectively)
+PopulationE <- rnorm(12, 0, sqrt(0.4))
+ContainerE <- rnorm(120, 0, sqrt(0.05))
+# generation of response values on latent scale (!) based on fixed effects,
+# random effects and residual errors
+EggL <- with(DataFemale, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 * (as.numeric(Habitat) -
+  1) + PopulationE[Population] + ContainerE[Container] + rnorm(480, 0, sqrt(0.1)))
+# data generation (on data scale!) based on Poisson distribution
+DataFemale$Egg <- rpois(length(EggL), exp(EggL))
+
+# Data frame for both sex
+DataAll <- Data
+# simulation of the underlying random effects (Population and Container with
+# variance of 0.5 and 0.8, respectively)
+PopulationE <- rnorm(12, 0, sqrt(0.5))
+ContainerE <- rnorm(120, 0, sqrt(0.8))
+# generation of response values on latent scale (!) based on fixed effects
+# and random effects
+ParasiteL <- with(DataAll, 1.8 + 2 * (-1) * (as.numeric(Sex) - 1) + 0.8 * (-1) *
+  (as.numeric(Treatment) - 1) + 0.7 * (as.numeric(Habitat) - 1) + PopulationE[Population] +
+  ContainerE[Container])
+# data generation (on data scale!) based on negative binomial distributions;
+# size = theta
+DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
+
+# simulation of the underlying random effects (Population and Container with
+# variance of 1.3 and 0.3, respectively)
+PopulationE <- rnorm(12, 0, sqrt(1.3))
+ContainerE <- rnorm(120, 0, sqrt(0.3))
+# data generation based on fixed effects, random effects and random
+# residuals errors
+DataAll$BodyL <- 15 + 3 * (-1) * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) -
+  1) + 0.15 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container] +
+  rnorm(960, 0, sqrt(1.2))
+
+# simulation of the underlying random effects (Population and Container with
+# variance of 0.2 and 0.2, respectively)
+PopulationE <- rnorm(12, 0, sqrt(0.2))
+ContainerE <- rnorm(120, 0, sqrt(0.2))
+# generation of response values on latent scale (!) based on fixed effects
+# and random effects
+ExplorationL <- with(DataAll, 4 + 1 * (-1) * (as.numeric(Sex) - 1) + 2 * (as.numeric(Treatment) -
+  1) + 0.5 * (-1) * (as.numeric(Habitat) - 1) + PopulationE[Population] +
+  ContainerE[Container])
+# data generation (on data scale!) based on gamma distribution; size = theta
+DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) *
+  0.3, rate = 0.3)
+
+# Subset the design matrix (only males express colour morphs)
+DataMale <- subset(Data, Sex == "Male")
+# simulation of the underlying random effects (Population and Container with
+# variance of 1.2 and 0.2, respectively)
+PopulationE <- rnorm(12, 0, sqrt(1.2))
+ContainerE <- rnorm(120, 0, sqrt(0.2))
+# generation of response values on latent scale (!) based on fixed effects
+# and random effects
+ColourL <- with(DataMale, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 *
+  (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+# data generation (on data scale!) based on binomial distribution
+DataMale$Colour <- rbinom(length(ColourL), 1, plogis(ColourL))
+
+
+# Quasi-Poisson with log link (page 5) -------------------------
+
+# Fit null model without fixed effects (but including all random effects)
+fecmodADMBr <- glmmadmb(Egg ~ 1 + (1 | Population) + (1 | Container), family = "nbinom1",
+data = DataFemale)
+# Fit alternative model including fixed and all random effects
+fecmodADMBf <- glmmadmb(Egg ~ Treatment + Habitat + (1 | Population) + (1 |
+Container), family = "nbinom1",data = DataFemale)
+
+# Calculation of the variance in fitted values
+VarF <- var(as.vector(model.matrix(fecmodADMBf) %*% fixef(fecmodADMBf)))
+# getting the observation-level variance Null model
+omegaN <- fecmodADMBr$alpha # overdispersion omega is alpha in glmmadmb
+lambda <- as.numeric(exp(fixef(fecmodADMBr) + 0.5 * (as.numeric(VarCorr(fecmodADMBr)[1]) + as.numeric(VarCorr(fecmodADMBr)[2]))))
+# lambda2 <- mean(DataFemale$Egg) # for lambda we use the mean of all
+# observations
+VarOdN <- omegaN / lambda # the delta method
+VarOlN <- log(1 + omegaN / lambda) # log-normal approximation
+VarOtN <- trigamma(lambda / omegaN) # trigamma function
+# comparing the three
+c(VarOdN = VarOdN, VarOlN = VarOlN, VarOlN = VarOtN)
+
+# Full model
+omegaF <- fecmodADMBf$alpha # overdispersion omega is alpha in glmmadmb
+VarOdF <- omegaF / lambda # the delta method
+VarOlF <- log(1 + omegaF / lambda) # log-normal approximation
+VarOtF <- trigamma(lambda / omegaF) # trigamma function
+# comparing the three
+c(VarOdF = VarOdF, VarOlF = VarOlF, VarOlF = VarOtF)
+
+# R2[GLMM(m)] - marginal R2[GLMM]
+R2glmmM <- VarF/(VarF + sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# R2[GLMM(c)] - conditional R2[GLMM] for full model
+R2glmmC <- (VarF + sum(as.numeric(VarCorr(fecmodADMBf)))) / (VarF + sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# Raw unadjusted ICC[Population]
+ICCrawPop <- VarCorr(fecmodADMBr)$Population[1]/(sum(as.numeric(VarCorr(fecmodADMBr))) + VarOlN)
+# adjusted ICC[Population]
+ICCadjPop <- VarCorr(fecmodADMBf)$Population[1]/(sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# Raw unadjusted ICC[Container]
+ICCrawCont <- VarCorr(fecmodADMBr)$Container[1]/(sum(as.numeric(VarCorr(fecmodADMBr))) + VarOlN)
+# adjusted ICC[Container]
+ICCadjCont <- VarCorr(fecmodADMBf)$Container[1]/(sum(as.numeric(VarCorr(fecmodADMBf))) + VarOlF)
+# comparing the results
+c(R2glmmM = R2glmmM, R2glmmC = R2glmmC, ICCrawPop = ICCrawPop, ICCadjPop = ICCadjPop, ICCrawCont = ICCrawCont, ICCadjCont = ICCadjCont)
+
+performance::r2_nakagawa(fecmodADMBf)