Nothing
context("testing utils functions")
# data
# set.seed(123)
p1 <- 10; nsample <- 10; n <- p1 * nsample
d <- data.frame(x1 = rnorm(n = n),
x2 = rnorm(n = n),
u1 = rep(1:p1, each = nsample),
u2 = rep(1:p1, times = nsample))
d$u1 <- as.factor(d$u1); d$u2 <- as.factor(d$u2)
# LMM: y with random intercept
b1 <- 1; b2 <- -1; sd1 <- 1.5
d$y_re_intercept <- b1 * d$x1 + b2 * d$x2 +
rep(rnorm(n = p1, sd = sd1), each = nsample) + # random intercept u1
rep(rnorm(n = p1, sd = sd1), times = nsample) + # random intercept u2
rnorm(n = n)
# LMM: y with random slope
b1 <- 0; sd1 <- 1; sd.x1 <- 2
d$y_re_slope <- b1 * d$x1 +
rep(rnorm(n = p1, sd = sd1), each = nsample) + # random intercept u1
d$x1 * rep(rnorm(n = p1, sd = sd.x1), times = nsample) + # random slope u1
rnorm(n = n)
# GLMM
b1 <- 1; sd1 <- 1.5
prob <- rr2::inv.logit(b1 * d$x1 + rep(rnorm(n = p1, sd = sd1), each = nsample)) # random intercept u1
d$y_binary <- rbinom(n = n, size = 1, prob = prob)
# y for PGLS
b1 <- 1.5; signal <- 0.7
phy <- ape::compute.brlen(ape::rtree(n = n), method = "Grafen", power = 1)
phy.x <- ape::compute.brlen(phy, method = "Grafen", power = .0001)
x_trait <- ape::rTraitCont(phy.x, model = "BM", sigma = 1)
e <- signal^0.5 * ape::rTraitCont(phy, model = "BM", sigma = 1) + (1-signal)^0.5 * rnorm(n=n)
d$x_trait <- x_trait[match(names(e), names(x_trait))]
d$y_pgls <- b1 * x_trait + e
rownames(d) <- phy$tip.label
# y for Phylogenetic logistic regression
b1 <- 1.5; signal <- 2
e <- signal * ape::rTraitCont(phy, model = "BM", sigma = 1)
e <- e[match(phy$tip.label, names(e))]
d$y_phy_binary <- rbinom(n = n, size = 1, prob = rr2::inv.logit(b1 * d$x1 + e))
test_that("inv.logit should return values between 0 and 1", {
expect_lt(inv.logit(10), 1)
expect_gt(inv.logit(-10), 0)
})
test_that("partial R2 should return values between 0 and 1", {
lm.x = lm(y_re_intercept ~ x1 + x2, data = d)
lm.y = lm(y_re_intercept ~ x1, data = d)
expect_lt(partialR2(lm.x, lm.y), 1)
expect_equal(partialR2(lm.x, lm.x), 0)
expect_gt(partialR2(lm.x, lm.y), 0)
})
# there is a typo in ape::binaryPGLMM, which was corrected in rr2::binaryPGLMM, so not equal now...
# test_that("binaryPGLMM within rr2 should have same results as binaryPGLMM from ape", {
# z.f1 <- ape::binaryPGLMM(y_phy_binary ~ x1, data = d, phy = phy)
# z.x1 <- ape::binaryPGLMM(y_phy_binary ~ 1, data = d, phy = phy)
# z.f <- rr2::binaryPGLMM(y_phy_binary ~ x1, data = d, phy = phy)
# z.x <- rr2::binaryPGLMM(y_phy_binary ~ 1, data = d, phy = phy)
# expect_identical(z.f1$s2, z.f$s2)
# expect_identical(z.f1$B, z.f$B)
# })
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.