library("testthat")
library("lme4")
eps <- .Machine$double.eps
oneMeps <- 1 - eps
set.seed(1)
## sample linear predictor values for the unconstrained families
etas <- list(seq.int(-8, 8, by=1), # equal spacing to asymptotic area
runif(17, -8, 8), # random sample from wide uniform dist
rnorm(17, 0, 8), # random sample from wide normal dist
c(-10^30, rnorm(15, 0, 4), 10^30))
## sample linear predictor values for the families in which eta must be positive
etapos <- list(seq.int(1, 20, by=1),
rexp(20),
rgamma(20, 3),
pmax(.Machine$double.eps, rnorm(20, 2, 1)))
## values of mu in the (0,1) interval
mubinom <- list(runif(100, 0, 1),
rbeta(100, 1, 3),
pmin(pmax(eps, rbeta(100, 0.1, 3)), oneMeps),
pmin(pmax(eps, rbeta(100, 3, 0.1)), oneMeps))
context("glmFamily linkInv and muEta")
test_that("inverse link and muEta functions", {
tst.lnki <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
sapply(frm, function(x) expect_that(fam$linkinv(x), equals(ff$linkInv(x))))
}
tst.muEta <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
sapply(frm, function(x) expect_that(fam$mu.eta(x), equals(ff$muEta(x))))
}
tst.lnki(binomial(), etas) # binomial with logit link
tst.muEta(binomial(), etas)
tst.lnki(binomial("probit"), etas) # binomial with probit link
tst.muEta(binomial("probit"), etas)
tst.lnki(binomial("cloglog"), etas) # binomial with cloglog link
tst.muEta(binomial("cloglog"), etas)
tst.lnki(binomial("cauchit"), etas) # binomial with cauchit link
tst.muEta(binomial("cauchit"), etas)
tst.lnki(poisson(), etas) # Poisson with log link
tst.muEta(poisson(), etas)
tst.lnki(gaussian(), etas) # Gaussian with identity link
tst.muEta(gaussian(), etas)
tst.lnki(Gamma(), etapos) # gamma family
tst.muEta(Gamma(), etapos)
tst.lnki(inverse.gaussian(), etapos) # inverse Gaussian
tst.muEta(inverse.gaussian(), etapos)
})
context("glmFamily linkFun and variance")
test_that("link and variance functions", {
tst.link <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
sapply(frm, function(x) expect_that(fam$linkfun(x), equals(ff$link(x))))
}
tst.variance <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
sapply(frm, function(x) expect_that(fam$variance(x), equals(ff$variance(x))))
}
tst.link( binomial(), mubinom)
tst.variance(binomial(), mubinom)
tst.link( binomial("probit"), mubinom)
tst.link( binomial("cauchit"), mubinom)
tst.link( gaussian(), etas)
tst.variance(gaussian(), etas)
tst.link( Gamma(), etapos)
tst.variance(Gamma(), etapos)
tst.link( inverse.gaussian(), etapos)
tst.variance(inverse.gaussian(), etapos)
tst.variance(MASS::negative.binomial(1), etapos)
tst.variance(MASS::negative.binomial(0.5), etapos)
tst.link( poisson(), etapos)
tst.variance(poisson(), etapos)
})
context("glmFamily devResid and aic")
test_that("devResid and aic", {
tst.devres <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
sapply(frm, function(x) {
nn <- length(x)
wt <- rep.int(1, nn)
n <- wt
y <- switch(fam$family,
binomial = rbinom(nn, 1L, x),
gaussian = rnorm(nn, x),
poisson = rpois(nn, x),
error("Unknown family"))
dev <- ff$devResid(y, x, wt)
expect_that(fam$dev.resids(y, x, wt), equals(dev))
dd <- sum(dev)
expect_that(fam$aic(y, n, x, wt, dd), equals(ff$aic(y, n, x, wt, dd)))
})
}
tst.devres(binomial(), mubinom)
tst.devres(gaussian(), etas)
tst.devres(poisson(), etapos)
})
context("negative binomial")
test_that("variance", {
tst.variance <- function(fam, frm) {
ff <- glmFamily$new(family=fam)
sapply(frm, function(x) expect_that(fam$variance(x), equals(ff$variance(x))))
}
tst.variance(MASS::negative.binomial(1.), etapos)
nb1 <- MASS::negative.binomial(1.)
cppnb1 <- glmFamily$new(family=nb1)
expect_that(cppnb1$theta(), equals(1))
nb2 <- MASS::negative.binomial(2.)
cppnb1$setTheta(2)
sapply(etapos, function(x) expect_that(cppnb1$variance(x), equals(nb2$variance(x))))
bfam <- glmFamily$new(family=binomial())
if (FALSE) {
## segfaults on MacOS mavericks 3.1.0 ... ??
expect_error(bfam$theta())#, "theta accessor applies only to negative binomial")
expect_error(bfam$setTheta(2))#, "setTheta applies only to negative binomial")
}
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.