Nothing
### test-mmm2.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: nov 29 2017 (15:22)
## Version:
## Last-Updated: Jan 12 2022 (14:46)
## By: Brice Ozenne
## Update #: 130
##----------------------------------------------------------------------
##
### Commentary:
## Test battery:
## - list of linear regressions:
## Compare multcomp:::glht to lavaSearch2:::glht2 (no small sample adjustment)
## * standard error derived from the information matrix
## * robust standard error derived from the iid decomposition
##
## - latent variable models:
## Compare multcomp:::glht to lavaSearch2:::glht2 (no small sample adjustment)
## * standard error derived from the information matrix, with or without second member
##
## - list of latent variable model (for linear regressions):
## Compare multcomp:::glht to lavaSearch2:::glht2 (no small sample adjustment)
## * standard error derived from the information matrix
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * header
## rm(list = ls())
if(FALSE){ ## already called in test-all.R
library(testthat)
library(lavaSearch2)
}
library(multcomp)
library(sandwich)
library(emmeans)
lava.options(symbols = c("~","~~"))
context("multcomp - mmm")
## * simulation
mSim <- lvm(c(Y1,Y2,Y3,Y4)~ beta * eta,
E ~ 1, Y1 ~ 0.25*T1 + 0.5*T2 + 0.05*T3)
latent(mSim) <- "eta"
set.seed(10)
n <- 1e2
df.data <- lava::sim(mSim, n, latent = FALSE, p = c(beta = 1))
df.data$eY1 <- exp(df.data$Y1)
## * linear regressions with logical constrains
e.lm <- lm(Y1 ~ T1 + T2 + T3, data = df.data)
e.lvm <- estimate(lvm(Y1 ~ T1 + T2 + T3), data = df.data)
## summary(e.lm)
test_that("glht vs. glht2 (logical constrains)", {
e.glht <- glht(e.lm, linfct = c("T2-T1=0",
"T2-T3=0",
"T1-T3=0"))
## summary(e.glht, test = adjusted("none"))
## summary(e.glht, test = adjusted("bonferroni"))
e.glht2 <- glht2(e.lvm, linfct = c("Y1~T2-Y1~T1=0",
"Y1~T2-Y1~T3=0",
"Y1~T1-Y1~T3=0"))
expect_equal(unname(e.glht$vcov),unname(e.glht2$vcov[1:4,1:4]), tol = 1e-6)
expect_equal(unname(e.glht$coef),unname(e.glht2$coef[1:4]), tol = 1e-6)
eS.glht <- summary(e.glht, test = adjusted("Shaffer"))
eS.glht2 <- summary(e.glht2, test = adjusted("Shaffer"))
expect_equivalent(eS.glht$test[c("coefficients","sigma","tstat","pvalues")],
eS.glht2$test[c("coefficients","sigma","tstat","pvalues")], tol = 1e-6)
})
test_that("glht2 (back-transformation)", {
e.log.lvm <- estimate(lvm(log(eY1) ~ T1 + T2 + T3), data = df.data)
e.glht2 <- glht2(e.log.lvm, linfct = c("eY1~T1","eY1~T2","eY1~T3"))
df.glht2 <- summary(e.glht2, transform = "exp", test = adjusted("none"))$table2
e.glht2.bis <- glht2(e.log.lvm, linfct = "eY1~T3")
df.glht2.bis <- summary(e.glht2.bis, transform = exp, test = adjusted("none"))$table2
expect_equal(as.double(df.glht2[3,]) , as.double(df.glht2.bis[1,]))
})
## * list of linear regressions
name.Y <- setdiff(endogenous(mSim),"E")[1:2]
n.Y <- length(name.Y)
ls.lm <- lapply(name.Y, function(iY){
eval(parse( text = paste("lm(",iY,"~E, data = df.data)")))
})
names(ls.lm) <- name.Y
class(ls.lm) <- "mmm"
ls.lvm <- lapply(name.Y, function(iY){
eval(parse( text = paste("estimate(lvm(",iY,"~E), data = df.data)")))
})
names(ls.lvm) <- name.Y
class(ls.lvm) <- "mmm"
test_that("glht vs. glht2 (list lm): information std", {
e.glht <- glht(ls.lm, mlf("E = 0"))
e.glht2 <- glht2(ls.lvm, linfct = "E")
e.glht2C <- glht2(ls.lvm, linfct = createContrast(ls.lvm, linfct = "E")$contrast)
eS.glht <- summary(e.glht)
eS.glht2 <- summary(e.glht2)
eS.glht2C <- summary(e.glht2C)
expect_equivalent(eS.glht2$test, eS.glht2C$test, tol = 1e-6)
expect_equal(unname(eS.glht$test$tstat), unname(eS.glht2$test$tstat), tol = 1e-6)
})
test_that("glht vs. glht2 (list ml): robust std", {
e.glht <- summary(glht(ls.lm, mlf("E = 0"), vcov = sandwich))
e.lava <- rbind(estimate(ls.lm[[1]])$coefmat[2,,drop=FALSE],
estimate(ls.lm[[2]])$coefmat[2,,drop=FALSE])
## no correction for the score
e.glht0 <- summary(glht2(ls.lvm, linfct = "E", robust = TRUE, ssc = "residuals0"))
## correction for the score by inflating the residuals such that they have correct variance
e.glht2 <- summary(glht2(ls.lvm, linfct = "E", robust = TRUE))
e.glht2C <- summary(glht2(ls.lvm, linfct = createContrast(ls.lvm, linfct = "E")$contrast, robust = TRUE))
expect_equivalent(e.glht0$test$tstat, e.glht$test$tstat, tol = 1e-6)
## cannot compare p.values
## because some are based on a student law and others on a gaussian law
expect_equivalent(e.glht2$test, e.glht2C$test, tol = 1e-6)
expect_equivalent(e.glht2$test$tstat, e.glht$test$tstat*sqrt(coef(estimate2(ls.lvm[[1]], ssc = "none"))["Y1~~Y1"])/sigma(ls.lm[[1]]), tol = 1e-6)
})
test_that("glht vs. calcDistMaxIntegral", {
e.glht <- glht(ls.lm, mlf("E = 0"), vcov = sandwich)
res.GS <- summary(e.glht)
iid.tempo <- do.call(cbind,lapply(ls.lm, iid)) %*% t(e.glht$linfct)
beta <- unlist(lapply(ls.lm, coef)) %*% t(e.glht$linfct)
beta.var <- diag(crossprod(iid.tempo))
z.value <- beta/sqrt(beta.var)
res.Search <- calcDistMaxIntegral(as.vector(z.value),
iid = iid.tempo, quantile.compute = FALSE,
df = NULL, trace = FALSE, alpha = 0.05)
expect_equal(as.double(res.Search$p.adjust),
as.double(res.GS$test$pvalues),
tolerance = attr(res.GS$test$pvalues, "error")
)
## cannot compare p.values
## because some are based on a student law and others on a gaussian law
})
## * lvm
## names(df.data)
m.lvm <- lvm(c(Y1,Y2,Y3)~ eta, eta ~ E,
Y1~Y4)
e.lvm <- estimate(m.lvm, df.data)
test_that("glht vs. glht2 (lvm): information std", {
resC <- createContrast(e.lvm, linfct = c("eta~E","Y2=1","Y3=1"))
e.glht.null <- glht(e.lvm, linfct = resC$contrast)
e.glht.H1 <- glht(e.lvm, linfct = resC$contrast, rhs = resC$null)
e.glht2.null <- glht2(e.lvm, linfct = resC$contrast, rhs = rep(0,3),
ssc = "none")
e.glht2.H1 <- glht2(e.lvm, linfct = resC$contrast, rhs = resC$null,
ssc = "none")
eS.glht.null <- summary(e.glht.null)
eS.glht.H1 <- summary(e.glht.H1)
eS.glht2.null <- summary(e.glht2.null)
eS.glht2.H1 <- summary(e.glht2.H1)
expect_equal(unname(eS.glht.null$test$tstat),
unname(eS.glht2.null$test$tstat))
expect_equal(unname(eS.glht.H1$test$tstat),
unname(eS.glht2.H1$test$tstat))
## cannot compare p.values
## because some are based on a student law and others on a gaussian law
})
## * list of lvm
mmm.lvm <- mmm(Y1 = estimate(lvm(Y1~E), data = df.data),
Y2 = estimate(lvm(Y2~E), data = df.data),
Y3 = estimate(lvm(Y3~E), data = df.data),
Y4 = estimate(lvm(Y4~E), data = df.data)
)
test_that("glht2 (list lvm): information std", {
##
resC <- createContrast(mmm.lvm, linfct = paste0("Y",1:4,": Y",1:4,"~E"))
lvm2.glht <- glht2(mmm.lvm, linfct = resC$contrast,
ssc = NA, robust = FALSE)
lvm2.sglht <- summary(lvm2.glht)
expect_equal(lvm2.sglht$df,100)
lvm3.glht <- glht2(mmm.lvm, linfct = resC$contrast,
rhs = rnorm(4),
ssc = NA, robust = FALSE)
lvm3.sglht <- summary(lvm3.glht)
expect_equal(lvm3.sglht$df,100)
})
##----------------------------------------------------------------------
### test-mmm2.R ends here
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.