Nothing
### test-auto-ttest.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: May 31 2021 (15:20)
## Version:
## Last-Updated: jul 31 2023 (18:07)
## By: Brice Ozenne
## Update #: 56
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
if(FALSE){
library(testthat)
library(lava)
library(multcomp)
library(nlme)
library(LMMstar)
}
context("Check lmm on examples of t-tests")
LMMstar.options(optimizer = "FS", method.numDeriv = "Richardson", precompute.moments = TRUE,
columns.confint = c("estimate","se","df","lower","upper","p.value"))
## * single t-test
## ** simulate data
n <- 15
p <- 3
X.name <- paste0("X",1:p)
link.lvm <- paste0("Y~",X.name)
formula.lvm <- as.formula(paste0("Y~",paste0(X.name,collapse="+")))
m <- lvm(formula.lvm)
distribution(m,~Id) <- Sequence.lvm(a = 1, b = n)
distribution(m,~Gender.num) <- binomial.lvm()
distribution(m,~X3) <- binomial.lvm(size=2)
transform(m,Gene~X3) <- function(x){factor(x,levels=0:2,labels=c("LL","LA","AA"))}
transform(m,Gender~Gender.num) <- function(x){factor(x,levels=0:1,labels=c("M","F"))}
transform(m,id~Gender.num) <- function(x){paste0("id",1:NROW(x))}
latent(m) <- ~Gender.num+X3
set.seed(10)
d <- lava::sim(m,n, latent = FALSE)
d$time <- "t1"
## ** test
test_that("single t-test",{
e.lmm <- lmm(Y ~ Gender-1, repetition = ~Gender|id, data = d, trace = 0, df = TRUE)
e2.lmm <- lmm(Y ~ Gender, repetition = ~Gender|id, data = d, trace = 0, df = TRUE)
## e3.lmm <- lmm(Y ~ Gender, repetition = ~Gender|id, structure = "UN", data = d, trace = 0)
e.gls <- gls(Y ~ Gender, weights = varIdent(form=~1|Gender), data = d)
e.ttest <- t.test(Y~Gender, data = d)
expect_equal(as.double(logLik(e.lmm)),as.double(logLik(e.gls)))
expect_equal(as.double(logLik(e2.lmm)),as.double(logLik(e.gls)))
expect_equal(sort(unname(e.ttest$estimate)), sort(unname(coef(e.lmm, effects = "mean"))))
expect_equal(e.ttest$p.value, anova(e.lmm, effects = c("GenderM-GenderF=0"))$multivariate$p.value, tol = 1e-5)
expect_equal(e.ttest$p.value, confint(e2.lmm)["GenderF","p.value"], tol = 1e-5)
## summary(e.gls)$tTable
## summary(e.lmm)$tTable
})
## summary(anova(e.lmm))
## summary(anova(e2.lmm))
## summary(e.lmm)
## * paired t-test
data(armd.wide, package = "nlmeU")
subjCC <- armd.wide[!is.na(armd.wide$visual0) & !is.na(armd.wide$visual52),"subject"]
armd.wideCC <- armd.wide[armd.wide$subject %in% subjCC,]
armd.longCC <- reshape2::melt(armd.wideCC,
id.var = c("subject","treat.f","lesion","miss.pat"),
measure.vars = c("visual0","visual4","visual12","visual24","visual52"),
variable.name = "week",
value.name = "visual")
armd.longCC$week <- factor(armd.longCC$week,
level = c("visual0","visual4","visual12","visual24","visual52"),
labels = c(0,4,12,24,52))
rownames(armd.longCC) <- NULL
test_that("paired t-test",{
## LMMstar.options(optimizer = "FS")
e.tt <- t.test(visual52-visual0 ~ treat.f, data = armd.wideCC)
e.lmm <- suppressMessages(lmm(visual ~ week*treat.f,
repetition = ~ week | subject, structure = "UN",
data = armd.longCC[armd.longCC$week %in% c("0","52"),]))
e.confintlmm <- confint(e.lmm)
expect_equivalent(abs(diff(e.tt$estimate)), abs(e.confintlmm["week52:treat.fActive","estimate"]), tol = 1e-5)
## expect_equivalent(e.tt$stderr, e.confintlmm["week52:treat.fActive","se"], tol = 1e-5) ## difference 0.01 (2.28 vs. 2.29)
## expect_equivalent(e.tt$parameter, e.confintlmm["week52:treat.fActive","df"], tol = 1e-5) ## difference 1.5 (191.5 vs 193)
})
## * multiple t-test
## ** simulate data
m <- lvm(c(Y1,Y2,Y3,Y4) ~ age + gender)
categorical(m, labels = c("male","female")) <- ~gender
transform(m, id~gender) <- function(x){1:NROW(x)}
distribution(m, ~age) <- gaussian.lvm(mean = 50, sd = 10)
set.seed(10)
dW <- lava::sim(m, 30)
## move to the long format
name.varying <- paste0("Y",1:4)
dL <- reshape(dW, direction = "long",
idvar = c("id","age","gender"),
varying = name.varying,
v.names = "Y",
timevar = "visit")
rownames(dL) <- NULL
dL$visit <- factor(dL$visit,
levels = 1:length(name.varying),
labels = name.varying)
## ** test
test_that("multiple t-test",{
## profvis::profvis(e.lmm <- lmm(Y ~ visit + gender:visit - 1, structure = UN(gender~visit|id), data = dL, trace = 0))
e.lmm <- lmm(Y ~ visit:gender - 1, repetition = gender~visit|id, structure = "UN", data = dL, trace = 0, df = TRUE)
e.lh <- anova(e.lmm, effects = c("gender:Y1"="visitY1:gendermale - visitY1:genderfemale = 0",
"gender:Y2"="visitY2:gendermale - visitY2:genderfemale = 0",
"gender:Y3"="visitY3:gendermale - visitY3:genderfemale = 0",
"gender:Y4"="visitY4:gendermale - visitY4:genderfemale = 0"))
## check likelihood and score
expect_equal(logLik(e.lmm), -240.1934, tol = 1e-5)
e.gls <- gls(Y ~ visit + visit:gender, correlation = corSymm(form=~1|id), weights = varIdent(form=~1|visit), data = dL)
ls.ttest <- list(t.test(Y1~gender, data = dW),
t.test(Y2~gender, data = dW),
t.test(Y3~gender, data = dW),
t.test(Y4~gender, data = dW))
## print(e.lh, method = "none")
expect_equal(confint(e.lh)$se, sapply(ls.ttest,"[[","stderr"), tol = 1e-5)
expect_equal(as.double(unlist(confint(e.lh, method = "none")$df)),
as.double(do.call(rbind,lapply(ls.ttest,"[[","parameter"))), tol = 1e-1)
expect_equal(as.double(unlist(confint(e.lh, method = "none")[,c("lower","upper")])),
as.double(do.call(rbind,lapply(ls.ttest,"[[","conf.int"))), tol = 1e-3)
})
##----------------------------------------------------------------------
### test-auto-ttest.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.