tests/testthat/test-rptPoisson.R

In mastoffel/rptR: Repeatability Estimation for Gaussian and Non-Gaussian Data

context("rptPoisson")

# Set a seed for reproducibility of the randomization 
suppressWarnings(RNGversion("3.5.0"))
set.seed(23)

# load data
data(BeetlesFemale)

########## checks for one random effect

# run with one random effect, no boot, no permut
R_est_1 <- rptPoisson(Egg ~ Treatment + (1|Container), grname=c("Container"), data = BeetlesFemale,
                      nboot=0, npermut=0)

test_that("rpt estimation works for one random effect, no boot, no permut, no parallelisation, logit link", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        expect_equal(R_est_1$R["R_org", ], 0.5449128, tolerance = 0.001)
        expect_equal(R_est_1$R["R_link", ], 0.5550295, tolerance = 0.001)
})

test_that("LRT works", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        expect_equal(R_est_1$P$LRT_P, 1.98e-46, tolerance = 0.001)
})


# run with one random effect, boot, no permut
R_est_2 <- rptPoisson(Egg ~ Treatment + (1|Container), grname=c("Container"), data = BeetlesFemale,
                      nboot=2, npermut=0)

test_that("rpt estimation works for one random effect, boot, no permut, no parallelisation, logit link", {
        
        expect_that(is.numeric(unlist(R_est_2$R)), is_true()) 
        expect_equal(R_est_2$R["R_org", ], 0.5449128, tolerance = 0.001)
        expect_equal(R_est_2$R["R_link", ], 0.5550295, tolerance = 0.001)
        
        # original scale
        expect_equal(as.numeric(R_est_2$CI_emp$CI_org["2.5%"]),  0.4124189, tolerance = 0.001)
        expect_equal(as.numeric(R_est_2$CI_emp$CI_org["97.5%"]), 0.553597, tolerance = 0.001)
        # link scale
        expect_equal(as.numeric(R_est_2$CI_emp$CI_link["2.5%"]), 0.424651, tolerance = 0.001)
        expect_equal(as.numeric(R_est_2$CI_emp$CI_link["97.5%"]), 0.5656179, tolerance = 0.001)
        
})




# run with one random effect, no boot, permut
R_est_3 <- rptPoisson(Egg ~ Treatment + (1|Container), grname=c("Container"), data = BeetlesFemale,
                      nboot=0, npermut=5)

test_that("rpt estimation works for one random effect, no boot, permut, no parallelisation, logit link", {
        
        expect_that(is.numeric(unlist(R_est_3$R)), is_true()) 
        expect_equal(R_est_3$R["R_org", ], 0.5449128, tolerance = 0.001)
        expect_equal(R_est_3$R["R_link", ], 0.5550295, tolerance = 0.001)
        
        # original scale
        expect_equal(R_est_3$P$P_permut_org, 0.2, tolerance = 0.001)
        # link scale
        expect_equal(R_est_3$P$P_permut_link, 0.2, tolerance = 0.001)
        
})






########## checks for two random effects

# run with two random effect, no boot, no permut
R_est_1 <- rptPoisson(Egg ~ Treatment + (1|Container) + (1|Population), grname=c("Container", "Population"), data = BeetlesFemale,
        nboot=0, npermut=0)

test_that("rpt estimation works for two random effect, no boot, no permut, no parallelisation, logit link", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        # 1st random effect
        expect_equal(R_est_1$R["R_org", 1], 0.03190105, tolerance = 0.001)
        expect_equal(R_est_1$R["R_link", 1],0.03799297, tolerance = 0.001)
        # 2nd random effect
        expect_equal(R_est_1$R["R_org", 2], 0.5041709, tolerance = 0.001)
        expect_equal(R_est_1$R["R_link", 2], 0.5195091, tolerance = 0.001)
})

test_that("LRTs works", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        expect_equal(R_est_1$P[1, "LRT_P"], 8.744869e-03, tolerance = 0.001)
        expect_equal(R_est_1$P[2, "LRT_P"], 1.187500e-16, tolerance = 0.001)
})


# check that random effect components plus overdispersion variance sum up to one
R_est_1 <- rptPoisson(Egg ~ Treatment + (1|Container) + (1|Population), 
        grname=c("Container", "Population", "Residual"), data = BeetlesFemale,
        nboot=0, npermut=0)

test_that("random effect components sum to up to one", {
        expect_equal(sum(R_est_1$R["R_link", ]), 1)
})


# run with two random effect, boot, no permut
R_est_2 <- rptPoisson(Egg ~ Treatment + (1|Container) + (1|Population), grname=c("Container", "Population"), data = BeetlesFemale,
        nboot=2, npermut=0)

test_that("rpt estimation works for two random effect, boot, no permut, no parallelisation, logit link", {
        
        # original scale
        # 1st random effect
        expect_equal(as.numeric(R_est_2$CI_emp$CI_org[1, "2.5%"]),  0.02431758, tolerance = 0.001)
        expect_equal(as.numeric(R_est_2$CI_emp$CI_org[1, "97.5%"]), 0.0256355, tolerance = 0.001)
        # 2nd random effect
        expect_equal(as.numeric(R_est_2$CI_emp$CI_org[2, "2.5%"]), 0.2346106 , tolerance = 0.001)
        expect_equal(as.numeric(R_est_2$CI_emp$CI_org[2, "97.5%"]),0.5839979, tolerance = 0.001)
        
        
        # link scale
        # 1st random effect
        expect_equal(as.numeric(R_est_2$CI_emp$CI_link[1, "2.5%"]), 0.02642277, tolerance = 0.001)
        expect_equal(as.numeric(R_est_2$CI_emp$CI_link[1, "97.5%"]), 0.03203479, tolerance = 0.001)
        # 2nd random effect
        expect_equal(as.numeric(R_est_2$CI_emp$CI_link[2, "2.5%"]),  0.2449164, tolerance = 0.001)
        expect_equal(as.numeric(R_est_2$CI_emp$CI_link[2, "97.5%"]), 0.5944401, tolerance = 0.001)
        
})



# run with two random effect, no boot, permut
R_est_3 <- rptPoisson(Egg ~ Treatment + (1|Container) + (1|Population), grname=c("Container", "Population"), data = BeetlesFemale,
        nboot=0, npermut=5)

test_that("rpt estimation works for two random effect, no boot, permut, no parallelisation, logit link", {
        # original scale
        expect_equal(R_est_3$P$P_permut_org[1], 0.2, tolerance = 0.001)
        expect_equal(R_est_3$P$P_permut_org[2], 0.2, tolerance = 0.001)
        # link scale
        expect_equal(R_est_3$P$P_permut_link[1], 0.2, tolerance = 0.001)
        expect_equal(R_est_3$P$P_permut_link[2], 0.2, tolerance = 0.001)
        
})




# test that Ratio works

# run with one random effect, no boot, no permut
R_est_1 <- rptPoisson(Egg ~ Treatment + (1|Container), grname=c("Container"), data = BeetlesFemale,
        nboot=5, npermut=5, ratio = FALSE)

test_that("rpt estimation works for one random effect, no boot, no permut, no parallelisation, logit link", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        expect_true(is.na(R_est_1$R["R_org", ]))
        expect_equal(R_est_1$R["R_link", ], 0.3229846, tolerance = 0.001)
})

test_that("LRT works", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        expect_equal(R_est_1$P$LRT_P, 1.98e-46, tolerance = 0.001)
})

# run with two random effect, no boot, no permut
R_est_2 <- rptPoisson(Egg ~ Treatment + (1|Container) + (1|Population), grname=c("Container", "Population"), data = BeetlesFemale,
        nboot=0, npermut=0, ratio = FALSE)

test_that("rpt estimation works for two random effect, no boot, no permut, no parallelisation, logit link", {
        expect_that(is.numeric(unlist(R_est_1$R)), is_true()) 
        # 1st random effect
        expect_true(is.na(R_est_2$R["R_org", 1]))
        expect_equal(R_est_2$R["R_link", 1],  0.02224455, tolerance = 0.001)
        # 2nd random effect
        expect_true(is.na(R_est_2$R["R_org", 2]))
        expect_equal(R_est_2$R["R_link", 2], 0.3041681, tolerance = 0.001)
})


# test that Residual and Overdispersion works

R_est_1 <-  rptPoisson(formula = Egg ~ Treatment + (1|Container) + (1|Habitat) ,
grname=c("Container", "Habitat", "Residual", "Overdispersion"), data = BeetlesFemale,
nboot=5, npermut=5, ratio = FALSE)

test_that("rpt estimation works for two random effects, estimation of Variance and Residual / Overdispersion", {
        expect_true(is.na(R_est_1$R["R_org", "Overdispersion"]))
        expect_true(is.na(R_est_1$R["R_org", "Residual"]))
        expect_equal(R_est_1$R["R_link", "Overdispersion"], 0.1004924, tolerance = 0.001)
        expect_equal(R_est_1$R["R_link", "Residual"],0.2570103, tolerance = 0.001)
})


R_est_1 <-  rptPoisson(formula = Egg ~ Treatment + (1|Container) + (1|Habitat) ,
        grname=c("Container", "Habitat", "Residual", "Overdispersion"), data = BeetlesFemale,
        nboot=0, npermut=0, ratio = TRUE)
R_est_2 <-  rptPoisson(formula = Egg ~ Treatment + (1|Container) + (1|Habitat) ,
        grname=c("Container", "Habitat"), data = BeetlesFemale,
        nboot=0, npermut=0, ratio = TRUE)


# test whether repeatabilities are equal for grouping factors independent of residual and overdispersion specification

test_that("repeatabilities are equal for grouping factors independent of residual and overdispersion specification", {
        expect_equal(R_est_1$R$Habitat, R_est_2$R$Habitat)
        expect_equal(R_est_1$R$Container, R_est_2$R$Container)
})


# check that grname sequence doesnt play a role
R_est_3 <-  rptPoisson(formula = Egg ~ Treatment + (1|Container) + (1|Habitat) ,
        grname=c("Habitat", "Container"), data = BeetlesFemale,
        nboot=0, npermut=0, ratio = TRUE)

test_that("Repeatabilities are equal for different order in grname argument", {
        expect_false(any(R_est_2$R$Container == R_est_3$R$Container) == FALSE)
        expect_false(any(R_est_2$R$Habitat == R_est_3$R$Habitat) == FALSE)
})

test_that("LRTs are equal for different different sequence in grname argument", {
        expect_equal(R_est_2$P["Container", "LRT_P"], R_est_3$P["Container", "LRT_P"])
        expect_equal(R_est_2$P["Habitat", "LRT_P"], R_est_3$P["Habitat", "LRT_P"])
})


R_est_4 <-  rptPoisson(formula = Egg ~ Treatment + (1|Habitat) + (1|Container),
        grname=c("Container", "Habitat"), data = BeetlesFemale,
        nboot=0, npermut=0, ratio = TRUE)

test_that("Repeatabilities are equal for different order in formula argument", {
        expect_false(any(R_est_2$R$Container == R_est_4$R$Container) == FALSE)
        expect_false(any(R_est_2$R$Habitat == R_est_4$R$Habitat) == FALSE)
})

test_that("LRTs are equal for different for different order in formula argument", {
        expect_equal(R_est_2$P["Container", "LRT_P"], R_est_4$P["Container", "LRT_P"])
        expect_equal(R_est_2$P["Habitat", "LRT_P"], R_est_4$P["Habitat", "LRT_P"])
})

# random slopes
R_est_5 <- rptPoisson(formula = Egg ~ Treatment + (1 + Treatment|Container),
        grname=c("Container"), data = BeetlesFemale,
        nboot=0, npermut=0, ratio = TRUE)

test_that("Random Slope repeatability point estimate is correct", {
        expect_equal(R_est_5$R["R_org", "Container"], 0.5539859, tolerance = 0.001)
})

test_that("Random slopes fitted without correlation (i.e. with grname occuring in multiple random effect terms) throw an error", {
        expect_error(rptPoisson(formula = Egg ~ (1|Treatment) + (0 + Treatment|Container),
                grname=c("Container"), data = BeetlesFemale,
                nboot=0, npermut=0, ratio = TRUE))
})