diff --git a/R/compute_variances.R b/R/compute_variances.R
index 35631d6bb..215b00477 100644
--- a/R/compute_variances.R
+++ b/R/compute_variances.R
@@ -20,7 +20,7 @@
   faminfo <- model_info(x, verbose = FALSE)
 
   # check argument
-  approx_method <- match.arg(approximation, c("lognormal", "delta", "trigamma"))
+  approx_method <- match.arg(approximation, c("lognormal", "delta", "trigamma", "observation_level"))
 
   if (any(faminfo$family == "truncated_nbinom1")) {
     if (verbose) {
@@ -601,11 +601,24 @@
       y_factor <- 1
     }
 
+    # for observation level approximation
+    fe_null <- as.numeric(lme4::fixef(model_null))
+    pmean <- as.numeric(stats::plogis(fe_null - 0.5 * sum(revar_null) * tanh(fe_null * (1 + 2 * exp(-0.5 * sum(revar_null))) / 6))) # nolint
+
     resid.variance <- switch(faminfo$link_function,
-      logit = pi^2 / (3 * y_factor),
-      probit = 1 / y_factor,
+      logit = switch(approx_method,
+        observation_level = 1 / (y_factor * pmean * (1 - pmean)),
+        pi^2 / (3 * y_factor)
+      ),
+      probit = switch(approx_method,
+        observation_level = 2 * pi / y_factor * pmean * (1 - pmean) * exp((stats::qnorm(pmean) / sqrt(2))^2)^2, # nolint
+        1 / y_factor
+      ),
       cloglog = ,
-      clogloglink = pi^2 / (6 * y_factor),
+      clogloglink = switch(approx_method,
+        observation_level = pmean / y_factor / log(1 - pmean)^2 / (1 - pmean),
+        pi^2 / (6 * y_factor)
+      ),
       .badlink(faminfo$link_function, faminfo$family, verbose = verbose)
     )
   } else if (faminfo$is_count) {
diff --git a/tests/testthat/test-r2_nakagawa_bernoulli.R b/tests/testthat/test-r2_nakagawa_bernoulli.R
index 8d9dd6622..acbb0620a 100644
--- a/tests/testthat/test-r2_nakagawa_bernoulli.R
+++ b/tests/testthat/test-r2_nakagawa_bernoulli.R
@@ -195,3 +195,72 @@ test_that("lme4, bernoulli", {
   expect_equal(out1[1, "R2m"], out2$R2_marginal, ignore_attr = TRUE, tolerance = 1e-4)
   expect_equal(out1[1, "R2c"], out2$R2_conditional, ignore_attr = TRUE, tolerance = 1e-4)
 })
+
+
+# ==============================================================================
+# Validate against Nakagawa et al. 2017 paper!
+test_that("glmer, Bernoulli", {
+  # example data from Nakagawa et al. 2017
+  Population <- gl(12, 80, 960)
+  Container <- gl(120, 8, 960)
+  Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))
+  Habitat <- factor(rep(rep(c("Dry", "Wet"), each = 4), 120))
+  Treatment <- factor(rep(c("Cont", "Exp"), 480))
+  Data <- data.frame(
+    Population = Population, Container = Container, Sex = Sex,
+    Habitat = Habitat, Treatment = Treatment
+  )
+  DataFemale <- Data[Data$Sex == "Female", ]
+  set.seed(777)
+  PopulationE <- rnorm(12, 0, sqrt(0.4))
+  ContainerE <- rnorm(120, 0, sqrt(0.05))
+  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)))
+  DataFemale$Egg <- rpois(length(EggL), exp(EggL))
+  DataAll <- Data
+  PopulationE <- rnorm(12, 0, sqrt(0.5))
+  ContainerE <- rnorm(120, 0, sqrt(0.8))
+  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])
+  DataAll$Parasite <- rnbinom(length(ParasiteL), size = 5, mu = exp(ParasiteL))
+  PopulationE <- rnorm(12, 0, sqrt(1.3))
+  ContainerE <- rnorm(120, 0, sqrt(0.3))
+  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))
+  PopulationE <- rnorm(12, 0, sqrt(0.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  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])
+  DataAll$Exploration <- rgamma(length(ExplorationL), shape = exp(ExplorationL) * 0.3, rate = 0.3)
+  DataMale <- subset(Data, Sex == "Male")
+  PopulationE <- rnorm(12, 0, sqrt(1.2))
+  ContainerE <- rnorm(120, 0, sqrt(0.2))
+  ColourL <- with(DataMale, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 * (as.numeric(Habitat) - 1) + PopulationE[Population] + ContainerE[Container])
+  DataMale$Colour <- rbinom(length(ColourL), 1, plogis(ColourL))
+
+  morphmodGLMERr <- lme4::glmer(
+    Colour ~ 1 + (1 | Population) + (1 | Container),
+    family = binomial(),
+    data = DataMale
+  )
+  # Fit alternative model including fixed and all random effects
+  morphmodGLMERf <- lme4::glmer(
+    Colour ~ Treatment + Habitat + (1 | Population) + (1 | Container),
+    family = binomial(),
+    data = DataMale
+  )
+
+  VarF <- var(as.vector(model.matrix(morphmodGLMERf) %*% lme4::fixef(morphmodGLMERf)))
+  VarDS <- pi^2 / 3
+  Vt <- lme4::VarCorr(morphmodGLMERr)$Population + lme4::VarCorr(morphmodGLMERr)$Container
+  pmean <- as.numeric(plogis(as.vector(lme4::fixef(morphmodGLMERr)) - 0.5 * Vt * tanh(as.vector(lme4::fixef(morphmodGLMERr)) * (1 + 2 * exp(-0.5 * Vt)) / 6)))
+  VarOL <- 1 / (pmean * (1 - pmean))
+
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarDS)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarDS)
+  out <- performance::r2_nakagawa(morphmodGLMERf)
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+
+  R2glmmM <- VarF / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarOL)
+  R2glmmC <- (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf)))) / (VarF + sum(as.numeric(lme4::VarCorr(morphmodGLMERf))) + VarOL)
+  out <- performance::r2_nakagawa(morphmodGLMERf, approximation = "observation_level")
+  expect_equal(out$R2_conditional, R2glmmC, tolerance = 1e-4, ignore_attr = TRUE)
+  expect_equal(out$R2_marginal, R2glmmM, tolerance = 1e-4, ignore_attr = TRUE)
+})