library("VGAM")
### tests for bracl and its agreement with other methods
## Stem cell research data from Agresti (2017, Analysis of ordinal data, Table 4.1)
freq <- c(34, 41, 58, 21, 30, 64,
67, 83, 63, 52, 52, 50,
30, 23, 15, 24, 18, 16,
25, 14, 12, 15, 11, 11)
fund_research <- factor(rep(c("definitely", "probably", "probably not", "definitely not"), each = 6),
levels = c("definitely", "probably", "probably not", "definitely not"),
ordered = TRUE)
gender <- factor(rep(rep(c("female", "male"), each = 3), 4), levels = c("male", "female"))
## gender <- rep(rep(c(1, 0), each = 3), 4)
religion <- factor(rep(c("fundamentalist", "moderate", "liberal"), 8),
levels = c("fundamentalist", "moderate", "liberal"),
ordered = TRUE)
stemcell <- data.frame(frequency = freq, research = fund_research, gender = gender, religion = religion)
ff <- matrix(freq, nrow = 4, byrow = TRUE)
re <- c(1, 2, 3, 1, 2, 3)
ge <- c(1, 1, 1, 0, 0, 0)
## With proportional odds
suppressWarnings(fit_vgam_p <- vglm(cbind(ff[1, ], ff[2, ], ff[3, ], ff[4, ]) ~ re + ge, family = acat(reverse = TRUE, parallel = TRUE)))
expect_warning(
fit_bracl_p <- bracl(research ~ as.numeric(religion) + gender, weights = frequency, data = stemcell, type = "ML", parallel = TRUE)
)
fit_vgam <- vglm(cbind(ff[1, ], ff[2, ], ff[3, ], ff[4, ]) ~ re + ge, family = acat(reverse = TRUE, parallel = FALSE))
## Without proportional odds
expect_warning(
fit_bracl <- bracl(research ~ as.numeric(religion) + gender, weights = frequency, data = stemcell, type = "ML")
)
tol <- 1e-06
## "VGAM::vglm and bracl return the same coefficients
expect_equal(unname(coef(fit_vgam)), unname(coef(fit_bracl)), tolerance = tol)
expect_equal(unname(coef(fit_vgam_p)), unname(coef(fit_bracl_p)), tolerance = tol)
## Difference in deviance is the same with VGAM::vglm and bracl
expect_equal(logLik(fit_vgam) - logLik(fit_vgam_p),
unclass(logLik(fit_bracl) - logLik(fit_bracl_p)), tolerance = tol,
check.attributes = FALSE)
## logLik returns the correct df for cracl
expect_identical(attr(logLik(fit_bracl), "df"), as.integer(9))
expect_identical(attr(logLik(fit_bracl_p), "df"), as.integer(5))
## bracl returns the correct fitted values
expect_equal(fitted(fit_bracl)[1:6, ], fit_vgam@fitted.values, tolerance = tol,
check.attributes = FALSE)
expect_equal(fitted(fit_bracl_p)[1:6, ], fit_vgam_p@fitted.values, tolerance = tol,
check.attributes = FALSE)
shu <- function(dat) dat[sample(seq.int(nrow(dat)), nrow(dat)), ]
## bracl results is invariance to shuffling of the data
for (j in 1:10) {
expect_warning(
fit_bracl_p_r <- bracl(research ~ as.numeric(religion) + gender, weights = frequency, data = shu(stemcell), type = "ML", parallel = TRUE)
)
expect_warning(
fit_bracl_r <- bracl(research ~ as.numeric(religion) + gender, weights = frequency, data = shu(stemcell), type = "ML", parallel = FALSE)
)
expect_equal(coef(fit_bracl), coef(fit_bracl_r), tolerance = tol)
expect_equal(coef(fit_bracl_p), coef(fit_bracl_p_r), tolerance = tol)
}
tol <- 1e-02
## vcov method for bracl returns the correct vcov matrix
expect_equal(vcov(fit_vgam_p), vcov(fit_bracl_p), tolerance = tol, check.attributes = FALSE)
expect_equal(vcov(fit_vgam), vcov(fit_bracl), tolerance = tol, check.attributes = FALSE)
s1 <- summary(fit_bracl)
s1p <- summary(fit_bracl_p)
s2 <- summary(fit_vgam)
s2p <- summary(fit_vgam_p)
tol <- 1e-06
## summary method for bracl returns the correct coef mat
expect_equal(coef(s1), s2@coef3, tolerance = tol, check.attributes = FALSE)
expect_equal(coef(s1p), s2p@coef3, tolerance = tol, check.attributes = FALSE)
newdata <- expand.grid(gender = c("male", "female"), religion = c("moderate", "fundamentalist"))
## predict.bracl works as expected
pp <- predict(fit_bracl_p, newdata = stemcell, type = "probs")
p <- predict(fit_bracl, newdata = stemcell, type = "probs")
expect_equal(predict(fit_vgam_p, type = "response"),
pp[19:24, ],
tolerance = 1e-06, check.attributes = FALSE)
expect_equal(predict(fit_vgam, type = "response"),
p[19:24, ],
tolerance = 1e-06, check.attributes = FALSE)
## no intercept returns warning
expect_warning(
fit_bracl_p_r <- bracl(research ~ -1 + as.numeric(religion) + gender, weights = frequency, data = shu(stemcell), type = "ML", parallel = FALSE)
)
## prediction with NAs works
newd <- newdata
newd[3, 2] <- NA
expect_true(is.na(predict(fit_bracl_p, newd, "class")[3]))
expect_true(all(is.na(predict(fit_bracl_p, newd, "probs")[3, ])))
## simulate from bracl objects returns a data frame with expected characteristics
simu_df <- simulate(fit_bracl_p)
nam_mf <- names(model.frame(fit_bracl_p))
nam_simu <- names(simu_df)
expect_identical(nrow(simu_df),
nrow(stemcell) * nlevels(stemcell$research))
expect_identical(levels(simu_df$research),
levels(stemcell$research))
expect_identical(is.ordered(simu_df$research),
is.ordered(stemcell$research))
expect_identical(nam_mf[!(nam_mf %in% nam_simu)],
"(weights)")
expect_identical(nam_simu[!(nam_simu %in% nam_mf)],
as.character(fit_bracl_p$call$weights))
expect_identical(nam_simu[(nam_simu %in% nam_mf)],
nam_mf[(nam_mf %in% nam_simu)])
## simulate from bracl objects returns a data frame with expected characteristics
simu_df <- simulate(fit_bracl)
nam_mf <- names(model.frame(fit_bracl))
nam_simu <- names(simu_df)
expect_identical(nrow(simu_df),
nrow(stemcell) * nlevels(stemcell$research))
expect_identical(levels(simu_df$research),
levels(stemcell$research))
expect_identical(is.ordered(simu_df$research),
is.ordered(stemcell$research))
expect_identical(nam_mf[!(nam_mf %in% nam_simu)],
"(weights)")
expect_identical(nam_simu[!(nam_simu %in% nam_mf)],
as.character(fit_bracl$call$weights))
expect_identical(nam_simu[(nam_simu %in% nam_mf)],
nam_mf[(nam_mf %in% nam_simu)])
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.