### test-previousBug.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: nov 19 2019 (10:17)
## Version:
## Last-Updated: jan 18 2022 (09:53)
## By: Brice Ozenne
## Update #: 35
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * header
## rm(list = ls(all.names = TRUE))
if(TRUE){ ## already called in test-all.R
library(testthat)
library(nlme)
library(lavaSearch2)
}
lava.options(symbols = c("~","~~"))
context("Previous bugs")
## * Brice, 01/23/20 1:41 , sCorrect
## keep.row <- sample.int(139,30)
## dd <- df.longiHCENRER[keep.row,c("neocortex.log","response.w4","age","sex","sert","sb.per.kg","group")]
## dd$sex <- as.numeric(dd$sex)
## dd$group <- as.numeric(dd$group)
## dd$sert <- as.numeric(dd$sert)
## dd$neocortex.log <- dd$neocortex.log+rnorm(NROW(dd), sd = 0.1)
ddW <- data.frame("neocortex.log" = c(-0.5114974, -0.8249681, -0.3910322, -0.2941140, -0.4850710, -0.7354951, -0.3005508, -0.3956807, -0.5279078, -0.4705333, -0.4072993, -0.5494491, -0.9464054, -0.2632149, -0.2781923, -1.0308939, -0.5637023, -0.2689296, -0.4375162, -0.5351115, -0.8410153, -0.3022184, -0.5711485, -0.5155027, -0.2262801, -0.4262566, -0.4353830, -0.7079919, -0.1529665, -0.3954827),
"group" = c(rep(1,10),rep(2,20)),
"id" = 1:30)
ddW$group <- factor(ddW$group, levels = 1:2)
e.GS <- gls(neocortex.log ~ group-1, method = "REML", weight = varIdent(form =~1|group), data = ddW)
## same as
e.lm1 <- lm(neocortex.log ~ 1, ddW[ddW$group==1, ])
e.lm2 <- lm(neocortex.log ~ 1, ddW[ddW$group==2, ])
ddL <- reshape2::dcast(ddW, value.var = "neocortex.log", id~group)
names(ddL) <- c("id","G1","G2")
m <- lvm(G1~1,G2~1)
e.lvm <- estimate(m, data = ddL, missing = TRUE)
test_that("sCorrect in stratified GLS equivalent to separate LM", {
eSSC.res <- estimate2(e.lvm)
## GS1 <- estimate2(lvm(G1~1), data = ddL[!is.na(ddL$G1),])
## GS2 <- estimate2(lvm(G2~1), data = ddL[!is.na(ddL$G2),])
GS <- c(mu1 = mean(ddL$G1, na.rm = TRUE),
mu2 = mean(ddL$G2, na.rm = TRUE),
sigma1 = var(ddL$G1, na.rm = TRUE),
sigma2 = var(ddL$G2, na.rm = TRUE))
expect_equivalent(coef2(eSSC.res), GS, tol = 1e-6)
## sigma(e.GS)^2
expect_equivalent(vcov2(eSSC.res)[1:2,1:2], vcov(e.GS), tol = 1e-6)
})
## * Brice, 01/27/20 6:12, ssc residuals under constrains
## path <- "C:/Users/hpl802/Downloads"
## butils.base:::sourcePackage("lavaSearch2", path = path, c.code = TRUE, trace = TRUE) ## version 1.5.5
dd <- data.frame("Y1" = c(-0.35, -0.87, -2.24, -0.7, 0.04, -1.46, -1.29, 0.6, -1.44, -1.64, -0.33, 1.12, -2, 0.66, 0.09, 1.18, -1.72, -1.02, 1.76, -0.48, -0.63, -1.95, -0.98, -2.8, -0.61),
"eta" = c(-0.37, -0.69, -0.87, -0.1, -0.25, -1.85, -0.08, 0.97, 0.18, -1.38, -1.44, 0.36, -1.76, -0.32, -0.65, 1.09, -0.76, -0.83, 0.83, -0.97, -0.03, 0.23, -0.3, -0.68, 0.66),
"Y2" = c(-0.77, -1.02, 0.5, 2.04, 0.25, -1.07, -0.98, 1.5, -0.46, -1.09, -2.67, -0.09, -2.59, 0.02, 0.41, 2.3, -0.03, -1.31, 1.4, -2.21, 0.35, -1.2, -1.35, -0.9, -0.83))
m <- lvm(Y1[0:sigma1] ~ 1*eta,
Y2[0:sigma2] ~ 1*eta,
eta[mu:1] ~ 1
)
latent(m) <- ~eta
mm <- lvm(Y1[mu:sigma1] ~ 1*eta,
Y2[mu:sigma2] ~ 1*eta,
eta[0:1] ~ 1
)
latent(mm) <- ~eta
e <- estimate(m, dd)
ee <- estimate(mm, dd)
test_that("bug in version 1.5.4 (incorrect handling of the constrain when computing Omega)", {
expect_equal(logLik(e), logLik(ee), tol = 1e-6)
expect_equal(as.double(vcov(e)), as.double(vcov2(e, ssc = FALSE)), tol = 1e-6)
expect_equal(as.double(vcov(ee)), as.double(vcov2(ee, ssc = FALSE)), tol = 1e-6)
test.res1 <- estimate2(e, df = "Satterthwaite", ssc = "residuals", derivative = "analytic")
test.res2 <- estimate2(ee, df = "Satterthwaite", ssc = "residuals", derivative = "analytic")
expect_equal(test.res1$sCorrect$param, c("eta" = -0.583625694846817, "Y1~~Y1" = 0.548634987452149, "Y2~~Y2" = 1.01248330967465),
tol = 1e-6)
## expect_equal(test.cox1$sCorrect$param, c("eta" = -0.58362569, "Y1~~Y1" = 0.50602817, "Y2~~Y2" = 0.95769503),
## tol = 1e-6)
})
## * Brice 04/15/20 9:26 multcomp
mSim <- lvm(Y1~X1,Y2~X1)
n <- 25
set.seed(10)
d <- lava::sim(mSim, n)
dNA <- rbind(d,c(Y1=NA,X1=1,Y2=2))
dNA$id <- 1:26
ls.lmALL <- list("Y1" = estimate(lvm(Y1~X1), data = d),
"Y2" = estimate(lvm(Y2~X1), data = d))
ls.lmRED <- list("Y1" = estimate(lvm(Y1~X1), data = d))
ls.lmNA <- list("Y1" = estimate(lvm(Y1~X1), data = dNA),
"Y2" = estimate(lvm(Y2~X1), data = dNA))
class(ls.lmALL) <- "mmm"
class(ls.lmRED) <- "mmm"
class(ls.lmNA) <- "mmm"
test_that("Stability by subset", {
glht.ALL <- glht2(ls.lmALL, linfct = "X1=0")
glht.RED <- glht2(ls.lmRED, linfct = "X1=0")
glht.NA <- glht2(ls.lmNA, linfct = "X1=0", cluster = "id")
index.model1 <- which(grepl("Y1: ",colnames(glht.ALL$vcov)))
expect_equal(glht.ALL$vcov[index.model1,index.model1], glht.RED$vcov, tol = 1e-6)
expect_equal(glht.RED$vcov[index.model1,index.model1], glht.NA$vcov[index.model1,index.model1], tol = 1e-6)
})
######################################################################
### test-previousBug.R ends here
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.