Nothing
tests/testthat/extralong/test-mv-extra.R
In spaMM: Mixed-Effect Models, with or without Spatial Random Effects
cat(crayon::yellow("\ntest-mv-extra:")) # not part of the testthat.R tests (neither test-composite-extra.R)
library(spaMM)
options(error=recover)
spaMM_tol <- local_tol <- spaMM.getOption("spaMM_tol")
local_tol$logL_tol <- 5e-07
spaMM.options(spaMM_tol=local_tol) # to control strictness of checks in independent-fit tests
{ # test missing data
set.seed(1)
d <- data.frame(y1 = rnorm(20), y2 = rnorm(20), x = runif(20))
d_y_missing <- d
d_y_missing[1, "y1"] <- NA
blaNA <- fitmv(submodels = list(one = list(y1 ~ x), two = list(y2 ~ x)),
# X2X=matrix(c(1,0,1,0,0,1,0,0,0,0,0,1),ncol=3,nrow=4,dimnames=list(NULL,c("(Intercept)","x_1","x_2"))),
data = d_y_missing)
testthat::expect_true(diff(range(predict(blaNA)[,1] - predict(blaNA, newdata=blaNA$data)[-1,1]))< 2e-16) # must work with X2X too
# predict(blaNA) does not predict for the line with NA response
# predict(blaNA, newdata=blaNA$data) predict for this line.
# It is important that the locvars provided to ..get_locdata() do not contain the response variable in the second case =>
# Cf .strip_cF_args(locformS[[mv_it]][-2]) in .calc_new_X_ZAC_mv()
testthat::expect_true(diff(range(predict(blaNA, blockSize=3)[,1] - predict(blaNA, newdata=blaNA$data, blockSize=3)[-1,1]))< 2e-16)
# there is slicing neither in first call (nrX = 0) nor in second (validrownames not NULL): cf conditions in predict.HLfit()
# With 'real' newdata without the validrownames attribute, slicing can occur:
testthat::expect_true(diff(range(predict(blaNA, newdata=d)[,1] - predict(blaNA, newdata=d, blockSize=3)[,1]))< 2e-16)
}
{
data("wafers")
me <- fitme(y ~ 1+(1|batch), family=Gamma(log), data=wafers)
set.seed(123)
y2 <- simulate(me, type="residual")
wafmv <- wafers
wafmv$batch2 <- wafmv$batch
wafmv$y2 <- y2
wafmv$ly <- log(wafmv$y)
wafmv$y3 <- log(y2)
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly ~ 1+(1|batch), family=gaussian()),
mod2=list(formula=y3 ~ 1+(1|batch2), family=gaussian())),
data=wafmv))
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14) ## there a resid.model so nothing is done with phi
simulate(zut1)
simulate(zut1, newdata=wafmv[1:3,])
cat(crayon::yellow("[ upper, get_ranPars() (-> VarCorr()) ] ; "))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly ~ 1+(1|batch), family=gaussian()),
mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian())),
data=wafmv, init=list(lambda=1), upper=list(lambda=0.02)))
testthat::expect_true(diff(range(get_ranPars(zut1, which="lambda"),0.02))<1e-12)
cat(crayon::yellow("(mv()|.) and (0+mv()|.); "))
(zut0 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(0+mv(1,2)|batch)),
mod2=list(formula=y3~X1+(0+mv(1,2)|batch), family=gaussian())),
data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(mv(1,2)|batch)),
mod2=list(formula=y3~X1+(mv(1,2)|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(mv(1,2)|batch)),
mod1=list(formula=ly~X1+(mv(1,2)|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut0),logLik(zut1),logLik(zut2)))<1e-05)
# # test equivalence of the two paremtrization, and that they are indeed distinguished by the code.
testthat::expect_true(diff(range(ranef(zut0)[[1]][,2]-rowSums(ranef(zut1)[[1]]))/(max(abs(ranef(zut0)[[1]]))))<1e-4)
(numinfo <- numInfo(zut0))
(numSEs <- sqrt(diag(solve(numinfo))))
condSEs <- summary(zut0,verbose=FALSE)$beta_table[,"Cond. SE"]
crit <- max(abs(numSEs[5:8]-condSEs))
testthat::test_that("numInfo() consistent with cond.SEs", testthat::expect_true(crit<1e-10))
{ # check of anova
zutLM <- fitmv(submodels=list(mod1=list(formula=ly~X1+batch),
mod2=list(formula=y3~1+batch, family=gaussian())),
data=wafmv)
crit <- max(abs(c(anova(fitme(ly~X1+batch, data=wafmv), method = "t.Chisq")[,3],
anova( fitme(y3~1+batch, data=wafmv), method = "t.Chisq")[,3]) -
anova(zutLM, type = "2")[,3]))
testthat::test_that("anova() OK on mv-LM",
testthat::expect_true(crit<1e-8))
}
#
anova(zut0, type = "2")
{ # missing data
afers <- wafmv
afers$y3[71:140] <- NA
afers$ly[1:70] <- NA
(zut0 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(0+mv(1,2)|batch)),
mod2=list(formula=y3~X1+(0+mv(1,2)|batch))),
data=afers))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(mv(1,2)|batch)),
mod2=list(formula=y3~X1+(mv(1,2)|batch))),
data=afers))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(mv(1,2)|batch)),
mod1=list(formula=ly~X1+(mv(1,2)|batch))),
data=afers))
logliks <- c(l0=logLik(zut0),l1=logLik(zut1),l2=logLik(zut2))
crit <- diff(range(logliks))
testthat::test_that(paste0("Hey, zut2 is still poor (inner estimation of phi?) logliks are", # so I made var_ranCoefs a reason for outer estim
paste(signif(logliks,6),collapse=",")),
testthat::expect_true(crit<1e-05)) # previous reasons for testing this were use_ZA_L or .calc_r22()
}
cat(crayon::yellow("ranCoefs; "))
(mod1 <- fitme(ly~X1, data=wafmv))
(mod2 <- fitme(y3~X1+(X2|batch), data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(X2|batch), family=gaussian()),
mod1=list(formula=ly~X1)),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-08)
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch)),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(X2|batch), family=gaussian()),
mod1=list(formula=ly~X1+(X2|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-07)
cat(crayon::yellow("fixing ranCoefs in two ways: "))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch), fixed=list(ranCoefs=list("1"=c(0.02, -0.1, 0.005)))),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch)),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv, fixed=list(ranCoefs=list("1"=c(0.02, -0.1, 0.005)))))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-08)
cat(crayon::yellow("isDiagFamily in two ways: "))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch), fixed=list(ranCoefs=list("1"=c(NA, 0, NA)))),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch)),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv, fixed=list(ranCoefs=list("1"=c(NA, 0, NA)))))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-08)
testthat::expect_true(unique(length(attr(zut1, "optimInfo")$optim.pars$trRanCoefs[["1"]]),
length(attr(zut2, "optimInfo")$optim.pars$trRanCoefs[["1"]]))==2L)
cat(crayon::yellow("with resid.model; "))
# independent-fit test
(mod1 <- fitme(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2), data=wafmv))
(mod2 <- fitme(formula=y2 ~ 1+(1|batch2), family=Gamma(log), data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2)),
mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(mod1)+logLik(mod2)))<2e-7)
{ # check of anova
mod1 <- fitme(formula=y ~ X1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2), data=wafmv)
mod2 <- fitme(formula=y2 ~ X1+(1|batch2), family=Gamma(log), data=wafmv)
zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ X1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2)),
mod2=list(formula=y2 ~ X1+(1|batch2), family=Gamma(log))),
data=wafmv)
crit <- max(abs(c(anova(mod1, method = "t.Chisq")[,3],
anova(mod2, method = "t.Chisq")[,3]) -
anova(zut1, type = "2")[,3]))
testthat::test_that("anova() OK on mv-GLMM",
testthat::expect_true(crit<1e-6))
}
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
cat(crayon::yellow("full dhglm (the 5 'phi' message are from different fits); "))
(mod1 <- fitme(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch), data=wafmv))
(mod2 <- fitme(formula=y2 ~ 1+(1|batch2), family=Gamma(log), data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch)),
mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log))),
data=wafmv))
get_residVar(zut1, newdata=zut1$data) # with a phiHGLM: residVar extractor has ignored this case for a long time.
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-4)
## To 'cover' the init procedure of refits in fitmv_body when there is a mixed-effect phi-resid.model:
# Here one must force some outer optimization for the mean-response model hence the lambda init:
cover_residM_reinit <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch))),
data=wafmv, init=list(lambda=c(0.01,0.01)))
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-4)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14) ## there a resid.model so nothin is done with phi
spaMM_boot(zut1, function(v) var(v), nsim=3L, type ="marginal")$bootreps
confint(zut1,"(Intercept)_1")
cat(crayon::yellow("fixing phi: four different ways; "))#
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.001)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.002))),
data=wafmv, fixed=list(lambda=0.1)))
# Check that <processed>$phi.Fix was correctly set up otherwise phi might still be fixed in fit but not considered as such in post-fit)
testthat::expect_true(all(sapply(zut1$phi.object, function(v) attr(v$phi_outer,"type"))==rep("fix",2)))
(zut2 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("1"=0.001,"2"=0.002), lambda=0.1)))
(zut3 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list(0.001,0.002), lambda=0.1)))
(zut4 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.001)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("2"=0.002), lambda=0.1)))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2),logLik(zut3),logLik(zut4)))<1e-14)
cat(crayon::yellow("visual checks; "))# and low fixed phis are useful for the following visual checks:
set.seed(123)
ressim1 <- simulate(zut1, type="residual")
ressim2 <- simulate(zut1, newdata=zut1$data, type="residual")
margsim1 <- simulate(zut1)
margsim2 <- simulate(zut1, newdata=zut1$data)
plot(ressim1,ressim2);abline(0,1) # clusters on the diagonal
plot(margsim1,margsim2);abline(0,1) # TWO indep draws of the 11 levels => clusters are not on the diagonal
# the two submodels share the ranef values so the clusters are always on the diagonal, even for margsim:
plot(ressim2[1:198], ressim2[1L+(1:198)]);abline(0,1) # on diagonal
plot(margsim2[1:198], margsim2[1L+(1:198)]);abline(0,1) # on diagonal
cat(crayon::yellow("confint with 'fixed'; "))# Checking that confint() obeys fixed values fixed in different ways
(zute <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(lambda=0.1)))
confint(zute,"(Intercept)_1")$lowerfit # estimated phis
(zutf <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.2224)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.2103))),
data=wafmv, fixed=list(lambda=0.1)))
confint(zutf,"(Intercept)_1")$lowerfit # fixed phis
(zutfm <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.2224)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("2"=0.2103), lambda=0.1)))
confint(zutfm,"(Intercept)_1")$lowerfit # fixed phis
(zutf2 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("2"=0.2103), lambda=0.1)))
confint(zutf2,"(Intercept)_1")$lowerfit # fixed 2nd phi
cat(crayon::yellow("some of the early tests; "))# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-8)
# permutation test. The init from hazards of development
(zut1 <- fitmv(submodels=list(mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian()),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log))),
control.HLfit=list(LevenbergM=TRUE),
data=wafmv,init=list(phi=list("1"=1 ))))
(zut2 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian())),
control.HLfit=list(LevenbergM=TRUE),
data=wafmv,init=list(phi=list("2"=1 ))))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
(zut <- fitmv(submodels=list(mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian()),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log))),
data=wafmv,init=list(phi=list("1"=1,"2"=0.2 ))))
attr(zut,"optimInfo")$LUarglist$canon.init
}
{
data("clinics")
#
(fitClinics <- HLfit(cbind(npos,nneg)~treatment+(1|clinic),family=binomial(),data=clinics))
set.seed(123)
y2 <- simulate(fitClinics, type="residual")
climv <- clinics
climv$np2 <- y2
climv$nn2 <- climv$npos + climv$nneg - y2
cat(crayon::yellow("binomial-poisson; "))# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(1|clinic),family=poisson())),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(1|clinic),family=poisson()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
cat(crayon::yellow("predVar cov gaussian-poisson; "))
# independent-fits test
(fg <- fitme(formula=npos~treatment+(1|clinic),family=gaussian(),data=climv))
(fp <- fitme(formula=np2~treatment+(1|clinic),family=poisson(),data=climv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=npos~treatment+(1|clinic),family=gaussian()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(+1|clinic),family=poisson()),
mod1=list(formula=npos~treatment+(1|clinic),family=gaussian())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2),logLik(fg)+logLik(fp)))<3e-7)
testthat::expect_true(diff(range(get_predVar(zut1)-c(get_predVar(fg),get_predVar(fp))))<1e-5)
testthat::expect_true(diff(range(get_predVar(zut2)-c(get_predVar(fp),get_predVar(fg))))<1e-5)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data,variances=list(cov=TRUE)) -
get_predVar(zut1,variances=list(cov=TRUE))))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data,variances=list(cov=TRUE)) -
Matrix::bdiag(get_predVar(fg,variances=list(cov=TRUE)),
get_predVar(fp,variances=list(cov=TRUE)))))<1e-5)
get_intervals(zut1)
get_intervals(zut1, intervals = "respVar")
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=npos~treatment+(1|clinic),family=gaussian()),
mod2=list(formula=np2~treatment+(1|clinic),family=poisson())),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(1|clinic),family=poisson()),
mod1=list(formula=npos~treatment+(1|clinic),family=gaussian())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
testthat::expect_true(diff(range(get_predVar(zut1)-get_predVar(zut2)[outer((1:16),c(16,0),"+")]))<1e-10)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data,variances=list(cov=TRUE)) -
get_predVar(zut1,variances=list(cov=TRUE))))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut2, newdata=zut2$data)-get_predVar(zut2)))<1e-14)
}
{
meth <- "ML(1,1)" # REML seems to work but the tests are OK only with lower accuracy
# independent-fit test
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv, method=meth))
(fp <- fitme(formula=np2~treatment+(1|clinic),family=poisson(), data=climv, method=meth))
(zut1 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
data=climv, method=meth))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(+1|clinic),family=poisson()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv, method=meth))
# test simple init
(zut3 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
init=list(lambda=c("1"=1.1,"2"=2.2)), data=climv, method=meth))
testthat::expect_true(identical(attr(zut3,"optimInfo")$LUarglist$canon.init, list(lambda=c("1"=1.1,"2"=2.2)))) # inits heeded
# test fancy names in init
zut4 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
init=list(lambda=c("clinic.1"=2.2,"clinic"=1.1)), data=climv, method=meth)
testthat::expect_true(identical(attr(zut4,"optimInfo")$LUarglist$canon.init, list(lambda=c("1"=1.1,"2"=2.2)))) # inits heeded
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(zut3), logLik(zut4), logLik(fb)+logLik(fp)))<2e-6) # 1e-5 for REML # stricter conv check decreased the precision of the comparison.
pfb <- get_predVar(fb)
pfp <- get_predVar(fp)
pzut1 <- get_predVar(zut1)
pzut1n <- get_predVar(zut1, newdata = zut1$data)
testthat::expect_true(diff(range(c(pfb,pfp)-pzut1, pzut1n-pzut1))<1e-5)
pzut2 <- get_predVar(zut2)
pzut2n <- get_predVar(zut2, newdata = zut2$data)
testthat::expect_true(diff(range(c(pfp,pfb)-pzut2, pzut2n-pzut2))<1e-5)
cat(crayon::yellow("deliberate warnings; ")) # (Deliberately generating warnings:)
oldopt <- options(warn=0L)
(zut5 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), init=list(lambda=1.1)),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson(), init=list(lambda=2.2))),
data=climv, method=meth))
options(oldopt)
testthat::expect_true(identical(attr(zut5,"optimInfo")$LUarglist$canon.init, NULL)) # inits NOT heeded as init is not preprocessed
cat(crayon::yellow("rand family and many post-fit fns; "))# rand.family independent-fit test
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(+1|clinic),family=poisson(), rand.family=Gamma(log)),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson(), rand.family=Gamma(log))),
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv))
(fg <- fitme(formula=np2~treatment+(1|clinic),family=poisson(), data=climv, rand.family=Gamma(log)))
testthat::expect_true(diff(range(logLik(zut1), logLik(fb)+logLik(fg)))<5e-6) # stricter conv check decreased the precision of the comparison.
cat(crayon::yellow("some extractors; "))
testthat::expect_true(diff(range(residuals(zut1)-c(residuals(fb),residuals(fg))))<1e-5)
formula(zut1) # list
terms(zut1) # list
nobs(zut1) # single value # decision made
# rand.family permutation test (rand.family need to be specified in each submodel)
(zut1 <- fitmv(submodels=list(mod1=list(cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), rand.family=Gamma(log)),
mod2=list(np2~treatment+(1|clinic),family=poisson(), rand.family=Gamma(log))),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(1|clinic),family=poisson(), rand.family=Gamma(log)),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), rand.family=Gamma(log))),
data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-9)
update_formulas(zut1,formula(zut1))
update(zut1, formula.=formula(zut1))
update(zut1, formula.=list(cbind(npos,nneg)~1+(1|clinic),
np2~1+(1|clinic)))
# step(zut1) dos not stop() but does not steps...
simulate(zut1,nsim=3) # checks that mv simulate Tpoisson works (possibly depending on same attributes as Tnegbin)
update_resp(zut1, newresp=simulate(zut1))
# CI for the variance of the random effect:
( ci <- confint(zut1,parm=function(fit){VarCorr(fit)[1,"Variance"]},
boot_args=list(nb_cores=9, nsim=9, seed=123)) )
# The distribution of bootstrap replicates:
plot(ecdf(ci$call$t))
zutnull <- update(zut1, submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), rand.family=Gamma(log)),
mod2=list(formula=np2~1+(1|clinic),family=poisson(), rand.family=Gamma(log))))
anova(zut1,zutnull)
anova(zut1,zutnull, boot.repl=99, nb_cores=9)
anova(zut1,zutnull, boot.repl=99, nb_cores=9)
# MSFDR(zutnull,zut1) # documented problem
cat(crayon::yellow("Tpoisson; "))# Tpoisson independent-fit test
(zut <- fitmv(submodels=list(mod2=list(formula=I(1L+np2)~treatment+(+1|clinic),family=Tpoisson()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=I(1L+np2)~treatment+(+1|clinic),family=Tpoisson())),
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv))
(fTp <- fitme(formula=I(1L+np2)~treatment+(1|clinic),family=Tpoisson(), data=climv))
testthat::expect_true(diff(range(logLik(zut), logLik(fb)+logLik(fTp)))<2e-6) # stricter conv check decreased the precision of the comparison.
## negbin(): outer-optimized dispersion parameters.
meth <- "ML(1,1)"
{ # quite slow 2*40s
# independent-fit test without fixed effects
(zut1 <- fitmv(submodels=list(mod2=list(formula=I(20*np2)~0+(+1|clinic),family=negbin()),
mod1=list(formula=cbind(npos,nneg)~0+(1|clinic),family=binomial())),
method=meth,
data=climv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~0+(1|clinic),family=binomial()),
mod2=list(formula=I(20*np2)~0+(+1|clinic),family=negbin())),
method=meth,
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~0+(1|clinic),family=binomial(),
method=meth,
data=climv))
(fn <- fitme(formula=I(20*np2)~0+(+1|clinic),family=negbin(),
method=meth,
data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(fb)+logLik(fn)))<1e-06)
}
cat(crayon::yellow("Tnegbin; "))## independent-fit test Tnegbin; outer-optimized dispersion parameters
(zut1 <- fitmv(submodels=list(mod2=list(formula=I(1L+20*np2)~treatment+(+1|clinic),family=Tnegbin()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=I(1L+20*np2)~treatment+(+1|clinic),family=Tnegbin())),
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv))
(fTn <- fitme(formula=I(1L+20*np2)~treatment+(1|clinic),family=Tnegbin(), data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(fb)+logLik(fTn)))<1e-06) # This did not work at first by obsInfo bc
# there were convergence problems for IRLS but LevM crit was not true => change in crit using pot4improv.
## permutation test Tnegbin; outer-optimized dispersion parameters
(zut1 <- fitmv(submodels=list(mod2=list(formula=I(1L+20*np2)~treatment+(1|clinic),family=Tnegbin()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=I(1L+20*np2)~treatment+(1|clinic),family=Tnegbin())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-06)
simulate(zut1,nsim=3) # checks that mv simulate Tnegbin keeps required attributes
{
cat(crayon::yellow("COMPoisson; "))## independent-fit test; outer-optimized dispersion parameters
data("freight") ## example from Sellers & Shmueli, Ann. Appl. Stat. 4: 943–961 (2010)
(mod1 <- fitme(broken ~ transfers+(1|id), data=freight, family = COMPoisson()))
freimv <- freight
set.seed(123)
freimv$brok2 <- simulate(mod1)
(mod2 <- fitme(brok2 ~ transfers,family=poisson(), data=freimv))
(zut1 <- fitmv(submodels=list(mod1=list(broken ~ transfers+(1|id),family=COMPoisson()),
mod2=list(brok2 ~ transfers,family=poisson())),
data=freimv))
(zut2 <- fitmv(submodels=list(mod2=list(brok2 ~ transfers,family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-06)
## permutation test COMPoisson; outer-optimized dispersion parameters
(zut1 <- fitmv(submodels=list(mod1=list(broken ~ transfers+(1|id),family=COMPoisson()),
mod2=list(brok2 ~ transfers+(1|id),family=poisson())),
data=freimv))
(zut2 <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-6)
simulate(zut1,nsim=3) # checks that mv simulate COMPoisson
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=COMPoisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("1"=0.5,"2"=1.1)))) # names needed !
attr(zutxx,"optimInfo")$LUarglist$init.optim$COMP_nu
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=COMPoisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
attr(zutxx,"optimInfo")$LUarglist$init.optim$COMP_nu
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
attr(zutxx,"optimInfo")$LUarglist$init.optim$COMP_nu
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson(1.1))),
data=freimv))
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson(1.1))),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=COMPoisson(1.1)),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
## assign("last.warning", NULL, envir = baseenv()) # flush the .COMP_maxn() warnings()
}
}
data("Loaloa")
lll <- Loaloa
lll$ID <- lll$ID2 <- 1+(seq(nrow(lll)) %% 2)
lll$resp <- 1+lll$npos + 10*lll$ID
(tnb <- fitme(resp~1+(1|ID), data=lll,family=Tnegbin()))
set.seed(124)
lll$r1 <- simulate(tnb)
lll$r2 <- simulate(tnb)
lll$long2 <- lll$longitude
## independent-fit test 2 Tnegbin
(tnb1 <- fitme(r1~1+(1|ID), data=lll,family=Tnegbin()))
(tnb2 <- fitme(r2~1+(1|ID), data=lll,family=Tnegbin()))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+(1|ID),family=Tnegbin()),
mod2=list(r2~1+(1|ID2),family=Tnegbin())),
data=lll))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+(1|ID2),family=Tnegbin()),
mod1=list(r1~1+(1|ID),family=Tnegbin())),
data=lll))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tnb1)+logLik(tnb2)))<1e-06)
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+(1|ID),family=Tnegbin()),
mod2=list(r2~1+(1|ID),family=Tnegbin())),
data=lll))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+(1|ID),family=Tnegbin()),
mod1=list(r1~1+(1|ID),family=Tnegbin())),
data=lll))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-8)
if (FALSE) { # slowish ~4s
cat(crayon::yellow("Matern; "))## independent-fit test with a Matern
(tp1 <- fitme(r1~1+Matern(1|longitude+latitude), data=lll,family=poisson()))
(tp2 <- fitme(r2~1+(1|ID), data=lll,family=poisson()))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson()),
mod2=list(r2~1+(1|ID2),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+(1|ID2),family=poisson()),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tp1)+logLik(tp2)))<1e-08)
}
if (FALSE) { # slow! > 4mn total
## independent-fit test with two Matern
(tp1 <- fitme(r1~1+Matern(1|longitude+latitude), data=lll,family=poisson()))
(tp2 <- fitme(r2~1+Matern(1|long2+latitude), data=lll,family=poisson()))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson()),
mod2=list(r2~1+Matern(1|long2+latitude),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+Matern(1|long2+latitude),family=poisson()),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tp1)+logLik(tp2)))<1e-08)
}
if (FALSE) { # slow! > 150s total
## independent-fit test with two Matern and fixed nu's
(tp1 <- fitme(r1~1+Matern(1|longitude+latitude), data=lll,family=poisson(), fixed=list(nu=0.5)))
(tp2 <- fitme(r2~1+Matern(1|long2+latitude), data=lll,family=poisson(), fixed=list(nu=1)))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=0.5)),
mod2=list(r2~1+Matern(1|long2+latitude),family=poisson(), fixed=list(nu=1))),
data=lll, verbose=c(TRACE=TRUE)))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+Matern(1|long2+latitude),family=poisson(), fixed=list(nu=1)),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=0.5))),
data=lll, verbose=c(TRACE=TRUE)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tp1)+logLik(tp2)))<1e-08)
}
if (FALSE) { # slow 22s
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=1)),
mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, verbose=c(TRACE=1)))
(zut2 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson()),
mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=1))),
data=lll, verbose=c(TRACE=1)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
predict(zut1)
simulate(zut1)
}
# default-name test
cat(crayon::yellow("'fixed' both in submodel and global call; "))
(zut1 <- fitmv(submodels=list(mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson()),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(rho=1,nu=1))),
data=lll, fixed=list(lambda=c("1"=666)), verbose=c(TRACE=1)))
(zut2 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(rho=1,nu=1)),
mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, fixed=list(lambda=c(666)), verbose=c(TRACE=1)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
{
data("blackcap")
MLdistMat2 <- as.matrix(proxy::dist(blackcap[,c("latitude","longitude")]))
MLcorMat2 <- MaternCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
nu=0.6285603,rho=0.0544659)
cap_mv <- blackcap
cap_mv$name <- as.factor(rownames(blackcap))
cap_mv$grp <- 1L+(blackcap$migStatus>1)
set.seed(123)
cap_mv$status2 <- blackcap$migStatus+ rnorm(14,sd=0.001)
cat(crayon::yellow("corrMatrix vs Matern; (map_ranef too)"))# corrMatrix
## independent-fit test with compar to equivalent Matern:
# need to fix this phi to avoid logLik uncertainty for low phi:
(mod1 <- fitme(migStatus ~ 1+ corrMatrix(1|name),data=cap_mv, corrMatrix=MLcorMat2, fixed=list(phi=0.1)))
(mod2 <- fitme(status2 ~ 1+ (1|grp),data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut1b <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1), corrMatrix=MLcorMat2),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
(zut1c <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, covStruct=list(corrMatrix=MLcorMat2,NULL,corrMatrix=NULL))) # a bit of abuse...
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ (1|grp)),
mod1=list(migStatus ~ 1+ corrMatrix(1|name),fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, fixed=list(rho=0.0544659,nu=0.6285603)))
(zut4 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
map_ranef(zut4)
simulate(zut1, newdata=cap_mv[1:3,])
simulate(zut1b, newdata=cap_mv[1:3,])
simulate(zut1c, newdata=cap_mv[1:3,])
simulate(zut2, newdata=cap_mv[1:3,])
simulate(zut3, newdata=cap_mv[1:3,])
simulate(zut4, newdata=cap_mv[1:3,])
get_predVar(zut4, variances=list(cov=TRUE))
if (FALSE) {
get_predVar(zut4, variances=list(cov=TRUE)) -> bla
bli <- bla
dim(bli) <- c(14,2,14,2)
bli[,1,,2]
bli[,1,,1]
bli[,2,,2]
}
testthat::expect_true(diff(range(logLik(zut1), logLik(zut1b), logLik(zut1c),
logLik(zut2), logLik(zut3), logLik(zut4), logLik(mod1)+logLik(mod2)))<1e-10)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
## permutation test with compar to equivalent Matern:
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ corrMatrix(1|name))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ corrMatrix(1|name)),
mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603)),
mod2=list(status2 ~ 1+ Matern(1|longitude+latitude))),
data=cap_mv))
(zut4 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ Matern(1|longitude+latitude)),
mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude),fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(zut3), logLik(zut4)))<1e-9)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
get_predVar(zut4, variances=list(cov=TRUE))
if (FALSE) { # may be used to check the mapping, but to have asym off-diag blocks, use migStatus ~ 1+ Matern(1|longitude+latitude)+(1|grp), the latter with fixed large variance
get_predVar(zut4, variances=list(cov=TRUE)) -> bli
nresp <- length(attr(bli, "respnames")) # but I commented out code providing this odd attribute, see get_predVar() and .predict_body()
nobs <- nrow(bli)/nresp
dim(bli) <- c(nobs,nresp,nobs,nresp)
dimnames(bli)[[2]] <- dimnames(bli)[[4]] <- attr(bli, "respnames")
bli[,1,,2]
bli[,1,,1]
bli[,2,,2]
bli[1,,1,]
}
cat(crayon::yellow("distMatrix; "))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Matern(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Matern(1|name))),
distMatrix=MLdistMat2, fixed=list(rho=0.0544659,nu=0.6285603),
data=cap_mv))
(zut2 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Matern(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Matern(1|name))),
distMatrix=2*MLdistMat2, fixed=list(rho=0.0544659/2,nu=0.6285603),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-9)
if (FALSE) { # Amusingly one can fit the same data by two sub-models (even tow diferent ones).
# we have twice the conditional likelihood, and once the ranef lik, so (as in any mv model with shared ranefs) marginal likelihoods do not add up.
# The phi's are low and they lower value must be controlled in order to make detailed numerical comparisons
(mod1 <- fitme(migStatus ~ 1+ corrMatrix(1|name),data=cap_mv, corrMatrix=MLcorMat))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name)),
mod2=list(migStatus ~ 1+ corrMatrix(1|name))),
data=cap_mv, corrMatrix=MLcorMat2))
vcov(zut1) # with phi->0 and a common ranef, this may be logical
}
cat(crayon::yellow("corrMatrix vs Cauchy; "))
MLcorMat3 <- CauchyCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
shape=1,longdep=0.5) # and default rho=1!
## independent-fit test with compar to equivalent Matern:
# need to fix this phi to avoid logLik uncertainty for low phi:
(mod1 <- fitme(migStatus ~ 1+ corrMatrix(1|name),data=cap_mv, corrMatrix=MLcorMat3, fixed=list(phi=0.1)))
(mod2 <- fitme(status2 ~ 1+ (1|grp),data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut1b <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1), corrMatrix=MLcorMat3),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
(zut1c <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, covStruct=list(corrMatrix=MLcorMat3,NULL,corrMatrix=NULL))) # a bit of abuse...
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ (1|grp)),
mod1=list(migStatus ~ 1+ corrMatrix(1|name),fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, fixed=list(shape=1,rho=1,longdep=0.5)))
(zut4 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude), fixed=list(phi=0.1,rho=1,shape=1,longdep=0.5)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut1b), logLik(zut1c),
logLik(zut2), logLik(zut3), logLik(zut4), logLik(mod1)+logLik(mod2)))<1e-10)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
## permutation test with compar to equivalent Cauchy:
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ corrMatrix(1|name))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ corrMatrix(1|name)),
mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude), fixed=list(phi=0.1,rho=1,shape=1,longdep=0.5)),
mod2=list(status2 ~ 1+ Cauchy(1|longitude+latitude))),
data=cap_mv))
(zut4 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ Cauchy(1|longitude+latitude)),
mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude),fixed=list(phi=0.1,rho=1,shape=1,longdep=0.5))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(zut3), logLik(zut4)))<1e-8)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
cat(crayon::yellow("distMatrix with Cauchy; "))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Cauchy(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Cauchy(1|name))),
distMatrix=MLdistMat2, fixed=list(rho=1,shape=1,longdep=0.5),
data=cap_mv))
(zut2 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Cauchy(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Cauchy(1|name))),
distMatrix=2*MLdistMat2, fixed=list(rho=1/2,shape=1,longdep=0.5),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-9)
}
{cat(crayon::yellow("IMRF; "))# fit IMRF
# There is a check of simulate of a large IMRF + Matern hurdle model
{ # create IMRF model
## Creating the mesh
oldMDCopt <- options(Matrix.warnDeprecatedCoerce = 0) # INLA issue
mesh <- INLA::inla.mesh.2d(loc = blackcap[, c("longitude", "latitude")],
cutoff=30,
max.edge = c(3, 20))
mesh$n ## 40
matern <- INLA::inla.spde2.matern(mesh)
options(oldMDCopt)
}
{ # Aim was first to check mv IMRF predictions, but putting a (G)LM as first model caught many problems
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
fixed=list(phi=c(0.02,0.02)),
data=cap_mv))
p1 <- predict(zut1)
p2 <- predict(zut1, newdata=zut1$data) # note difference in frame attribute. What do we want? _FIXME_
testthat::expect_true(diff(range(p2-p1))<1e-14)
}
{ # independent-fit test
# a bit tricky bc nloptr does not find the univariate optimum without some help (... tiny blackcap data)
# and in mv case, it finds it by default or not depending on order...
(mod1 <- fitme(migStatus ~ means + IMRF(1|longitude+latitude, model=matern), init=list(lambda=0.02),
data=cap_mv, verbose=c(TRACE=FALSE)))
(mod2 <- fitme(status2 ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603), data=cap_mv,
verbose=c(TRACE=FALSE)))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern)),
mod2=list(status2 ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1))),
data=cap_mv, fixed=list(rho=0.0544659,nu=0.6285603), init=list(lambda=c(0.02,NA)))) # Matern in spprec...
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1)),
mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern))),
data=cap_mv, fixed=list(rho=0.0544659,nu=0.6285603), init=list(lambda=c(NA,0.02)))) # Matern in spprec...
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-08)
}
{ # permutation test
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern)),
mod2=list(status2 ~ means+ IMRF(1|longitude+latitude, model=matern))),
data=cap_mv))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ means+ IMRF(1|longitude+latitude, model=matern)),
mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern))),
data=cap_mv)) # sensitive to .solve_crossr22( ., use_crossr22=TRUE)
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-07)
# The following with status2 ~ 1+ ... rather than status2 ~ means+ ... makes lambda diverge (jointly with kappa: does not occur for fixed kappa).
# Then all comparisons are sensitive to floating point precision (and no stringent permutation test is passed).
# but let us keep the zut2_testr22 computation as it has been useful in the past to 'torture-test' .calc_r22() algorithms
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern)),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
init=list(lambda=c(0.02),phi=c(0.02,0.02)),
data=cap_mv))
crit <- diff(range(logLik(zut1), -6.48825776602))
try(testthat::test_that(paste0("criterion was ",signif(crit,6)," from -6.48825776602"), # affected by use_ZA_L or .calc_r22() ... and minKappa...
testthat::expect_true(crit<1e-08)))
(zut2_testr22 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern)),
mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern))),
init=list(lambda=c(0.02),phi=c(0.02,0.02)),
data=cap_mv)) # has been sensitive to .solve_crossr22( ., use_crossr22=TRUE)
crit <- diff(range(-6.48825776601 , logLik(zut2_testr22)))
try(testthat::test_that(paste0("criterion was ",signif(crit,6)," from -6.48830317593"), # affected by use_ZA_L or .calc_r22() ... and minKappa...
testthat::expect_true(crit<1e-08)))
}
if (spaMM.getOption("example_maxtime")>87) { cat(crayon::yellow("multIMRF indep-fit tests; "))
(mrf1fixx <- fitme(migStatus ~ 1 + (1|pos) +
multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(phi=1,lambda=c("1"=0.5),
hyper=list("1"=list(hy_kap=0.1,hy_lam=1)))) )
(mrf2 <- fitme(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern),
data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + (1|pos) + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
fixed=list(phi=1,lambda=c("1"=0.5),
hyper=list("1"=list(hy_kap=0.1,hy_lam=1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fixx)+logLik(mrf2)))<1e-7) # precision modified (effect of changing bobyqa's rhoend?)
#
(mrf1fix <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(hyper=list("1"=list(hy_kap=0.1,hy_lam=1)))) )
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
fixed=list(hyper=list("1"=list(hy_kap=0.1,hy_lam=1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fix)+logLik(mrf2)))<1e-10)
(mrf1fix <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(hyper=list("1"=list(hy_kap=0.1)))) )
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
fixed=list(hyper=list("1"=list(hy_kap=0.1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fix)+logLik(mrf2)))<1e-10)
(mrf1fix <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(hyper=list("1"=list(hy_lam=1)))) )
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
fixed=list(hyper=list("1"=list(hy_lam=1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fix)+logLik(mrf2)))<1e-10)
{ # particularly slow
{
# now there is a global maximum at phi=0 which may be missed when phi outer optim is used,
# But not missed when default augZXy is used. augZXy is used by default for mrf1, but not in mv => discrepancy
spaMM.options(allow_augZXy=FALSE)
logLik(mrf1F <- fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv)) # bad
logLik(fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv, init=list(phi=(1e-4)))) # good
spaMM.options(allow_augZXy=NULL) #
# This may lose accuracy when .calc_r22() is modified! :
logLik(mrf1T <- fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv)) # good; use y-augmented matrix
logLik(mrf1I <- fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv, init.HLfit=list(phi=1))) # good; use y-augmented matrix !
# so to compare with mv fit we must provide a low init phi to the submodel
}
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut0 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
data=cap_mv, init=list(phi=list("1"=1e-4))))
testthat::expect_true(diff(range(logLik(zut0), logLik(mrf1T)+logLik(mrf2)))<1e-5)
(mrf3 <- fitme(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern), data=cap_mv))
# inits to try to speed the fit () (and lambda=c("3"=.) works as it should)
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
data=cap_mv, init=list(lambda=c("3"=5),phi=list("1"=1e-4,"2"=0.05)),verbose=c(TRACE=0.5)))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern)),
mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(lambda=c("1"=5),phi=list("1"=0.05,"2"=1e-4))))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mrf1T)+logLik(mrf3)))<1e-5)
}
} else cat(crayon::bgGreen("\n multIMRF indep-fit tests are slow (~87s). Run them once in a while; "))
{ cat(crayon::yellow("multIMRF permutation tests; ")) # ~ 13s
# reason for init phi as above: to avoid a local maximum
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1+ multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(phi=list("1"=1e-4,"2"=1e-4))))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(phi=list("1"=1e-4,"2"=1e-4))))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
if (FALSE) {
# The following code fits a single (repeated) kappa value, and only one lambda hyperparam.
# So the single IMRF from multIMRF(. ,levels=1) is possibly
# recognized as first level of multIMRF(. ,levels=2). Which may be nice, or may not be the intent. (__FIXME__).
# (What does .calc_normalized_ZAlist() do ?)
# Adding a fictitious argument bla=666 has no effect as it is not retained in the expanded formula.
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ multIMRF(1|longitude+latitude,margin=5,levels=1)),
mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(phi=list("1"=1e-4,"2"=1e-4))))
}
}
}
{ cat(crayon::yellow("adjacency; "))
data("scotlip")
(mod1 <- fitme(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)),
adjMatrix=Nmatrix, family=poisson(), data=scotlip) )
scotmv <- scotlip
set.seed(123)
scotmv$cases2 <- simulate(mod1)
scotmv$code2 <- scotmv$gridcode
(mod2 <- fitme(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)),
adjMatrix=Nmatrix, family=poisson(), data=scotmv) )
{ # independent-fit test
# a bit tricky bc nloptr does not find the univariate optimum without some help (... tiny blackcap data)
# and in mv case, it finds it by default or not depending on order...
(zut1 <- fitmv(submodels=list(mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson()),
mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|code2)+offset(log(expec)), family=poisson())),
data=scotmv,covStruct=list(adjMatrix=Nmatrix,adjMatrix=Nmatrix))) # Matern in spprec...
(zut2 <- fitmv(submodels=list(mod1=list(cases2 ~ I(prop.ag/10)+adjacency(1|code2)+offset(log(expec)), family=poisson()),
mod2=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson())),
data=scotmv,covStruct=list(adjMatrix=Nmatrix,adjMatrix=Nmatrix))) # Matern in spprec...
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-04) # __F I X M E__ was more accurate by outer lambda
}
{ # permutation test
# here to the order affects nloptr... tiny blackcap data again. Note lambda divergence.
# we can avoid that by either providing fixed or init phi (the latter slow)
(zut1 <- fitmv(submodels=list(mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson()),
mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson())),
data=scotmv,adjMatrix=Nmatrix))
(zut2 <- fitmv(submodels=list(mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson()),
mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson())),
data=scotmv,adjMatrix=Nmatrix))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
#
(zut1 <- fitmv(submodels=list(mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec))),
mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)))),
data=scotmv,adjMatrix=Nmatrix))
confint(zut1, parm="(Intercept)_1")
}
}
{
data("sleepstudy",package = "lme4")
(mod0 <- fitme(Reaction ~ Days + AR1(1|Days), data = sleepstudy, fixed=list(ARphi=0.5,lambda=1000)))
# bc the unconstrained fit has essentially lambda=0, which cause later ambiguities
# Don't give to much importance to the result: the fixed Days capture the trend and the 'correlated' rnaefs are not really correlated
# plot(ranef(mod0)[[1]])
sleepmv <- sleepstudy
set.seed(124)
sleepmv$Reaction <- simulate(mod0)
sleepmv$reac2 <- simulate(mod0)
sleepmv$days2 <- sleepmv$Days
cat(crayon::yellow("non-standard REML (both ways); "))
(re1 <- fitme(Reaction ~ 1 + (1|Days), REMLformula=~Days, data = sleepmv, method="REML"))
(re2 <- fitme(reac2 ~ Days + (1|days2), REMLformula=~1, data = sleepmv, method="REML"))
(rezut1 <- fitmv(submodels=list(mod1=list(Reaction ~ 1 + (1|Days), REMLformula=~Days),
mod2=list(reac2 ~ Days + (1|days2), REMLformula=~1)),
data=sleepmv, method="REML"))
testthat::expect_true(diff(range(logLik(rezut1), logLik(re1)+logLik(re2)))<1e-08)
cat(crayon::yellow("AR1; "))
(mod1 <- fitme(Reaction ~ Days + AR1(1|Days), data = sleepmv))
(mod2 <- fitme(reac2 ~ Days + AR1(1|Days), data = sleepmv))
(zut1 <- fitmv(submodels=list(mod1=list(Reaction ~ Days + AR1(1|Days)),
mod2=list(reac2 ~ Days + AR1(1|days2))),
data=sleepmv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mod1)+logLik(mod2)))<1e-08)
(zut1 <- fitmv(submodels=list(mod1=list(Reaction ~ Days + AR1(1|Days)),
mod2=list(reac2 ~ Days + AR1(1|Days))),
data=sleepmv))
(zut2 <- fitmv(submodels=list(mod2=list(reac2 ~ Days + AR1(1|Days)),
mod1=list(Reaction ~ Days + AR1(1|Days))),
data=sleepmv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
}
if(FALSE) { cat(crayon::yellow("simulation study; "))
{
set.seed(123)
replic <- function(nind=1000L, cor=-0.5, lambda=c(1,0.2), binsize=1L, method=if (binsize<4L) {"PQL/L"} else {"ML"}, return.fit=FALSE, ...) {
cat(".")
u1 <- rnorm(nind)
u2 <- rnorm(nind)
lam1 <- lambda[1]
lam2 <- lambda[2]
L <- chol(matrix(c(lam1,sqrt(lam1*lam2)*cor,sqrt(lam1*lam2)*cor,lam2),ncol=2))
v <- t(L) %*% rbind(u1,u2)
eta1 <- v[1,]
eta2 <- 1+v[2,]
surv <- rbinom(length(eta1),binsize,inv.logit(eta1)) # was it boot::inv.logit ?? Also called once more below.
feco <- rpois(length(eta2), lambda = exp(eta2))
lfh <- data.frame(id=seq_along(eta1), surv=surv, feco=feco, binsize=binsize)
(fitlfh <- fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id), family=poisson())),
data=lfh, method=method, ...))
if (return.fit) {
return(fitlfh)
} else {
corr <- VarCorr(fitlfh)[2,"Corr."]
cat(corr)
corr
}
}
replicate(3L, replic(nind=1000L))
replicate(10L, replic(nind=3000L, binsize=10)) # much better...
replicate(10L, replic(nind=300L, binsize=10)) # much better...
PQLLdist <- replicate(50L, replic(nind=300L, binsize=2)) # passable
PQLdist <- replicate(50L, replic(nind=300L, binsize=2, method="PQL")) # passable
# one can guess a *small* advantage of REML:
plot(ecdf(PQLdist))
plot(ecdf(PQLLdist), add=TRUE, col="red")
#
PQLdist <- replicate(10L, replic(nind=1000L, binsize=1, method="PQL")) # binary...
PQLdist <- replicate(10L, replic(nind=1000L, cor=0, binsize=1, method="PQL")) # binary...
PQLdist <- replicate(10L, replic(nind=300L, cor=-0.5, binsize=1, method="PQL")) # binary & small sample size -> extreme negative
PQLdist <- replicate(10L, replic(nind=10000L, cor=-0.5, binsize=1, method="PQL")) # large sample size does not improve this
(zut <- replic(nind=300L, cor=-0.5, binsize=1, method="PQL", return.fit=TRUE, verbose=c(TRACE=TRUE)))
replic(nind=1000L, binsize=1, return.fit=TRUE)
plot(eta1,eta2)
cov(t(v))
plot(eta1,surv)
plot(eta2,feco)
plot(eta1,feco)
cov(cbind(surv,feco))
}
{
set.seed(123)
replic2 <- function(nind=1000L, lambda=0.2, binsize=1L, method=if (binsize<4L) {"PQL/L"} else {"ML"}, return.fit=FALSE, ...) {
cat(".")
u1 <- sqrt(lambda)*rnorm(nind)
eta_f <- 1+u1
feco <- rpois(length(eta_f), lambda = exp(eta_f))
surv <- rbinom(length(eta_f),binsize,inv.logit(-(log(feco)-1)))
lfh <- data.frame(id=seq_along(eta_f), id2=seq_along(eta_f), surv=surv, feco=feco, binsize=binsize)
(fitlfh <- fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id), family=poisson())),
data=lfh, method=method, ...))
if (return.fit) {
return(fitlfh)
} else {
corr <- VarCorr(fitlfh)[2,"Corr."]
cat(corr)
corr
}
}
replicate(10L, replic2(nind=1000L))
replicate(10L, replic2(nind=1000L, binsize=10))
replicate(10L, replic2(nind=3000L))
replicate(10L, replic2(nind=3000L, binsize=10)) #
(zut <- replic2(nind=3000L, binsize=10, return.fit=TRUE))
with(zut$data, plot(surv,feco))
fitme(cbind(surv,binsize-surv) ~ 1+(1|id), family=binomial(),
data=zut$data)
fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(1|id), family=binomial()),
list(feco ~ 1+(1|id2), family=poisson())),
data=zut$data)
fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id), family=poisson())),
data=zut$data)
(debu <- fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id2), family=poisson())),
data=zut$data[1:20,]))
plot(eta1,eta2)
cov(t(v))
plot(eta1,surv)
plot(eta2,feco)
plot(eta1,feco)
cov(cbind(surv,feco))
}
{
simfun <- function(nind=570L, cor=-0.5, lambda=c(0.2,0.1), return.fit=FALSE, ...) {
cat(".")
u <- rnorm(2*nind)
lam1 <- lambda[1]
lam2 <- lambda[2]
L <- t(chol(matrix(c(lam1,sqrt(lam1*lam2)*cor,sqrt(lam1*lam2)*cor,lam2),ncol=2)))
v <- L %*% matrix(u,nrow=2)
if (return.fit) plot(t(v))
lfh <- data.frame(id=seq_len(nind), id2=seq_len(nind), feco= rpois(nind, lambda = exp(1+v[1,])),
growth=rgamma(nind,shape=1/0.2, scale=0.2*exp(1+v[2,]))) # mean=exp(1+v[2,]), var= 0.2*mean^2
(fitlfh <- fitmv(submodels=list(list(feco ~ 1+(0+mv(1,2)|id), family=poisson()),
list(growth ~ 1+(0+mv(1,2)|id), family=Gamma(log))),
data=lfh, method=method, ...))
if (return.fit) {
return(fitlfh)
} else {
corr <- VarCorr(fitlfh)[2,"Corr."]
cat(corr)
corr
}
}
set.seed(123)
(zut <- simfun(cor=-0.5,return.fit=TRUE))
plot(zut$data)
replicate(3L, simfun(3000))
replicate(3L, simfun(10000))
}
}
if (FALSE) {
library(aster)
data(echinacea)
asNA.acyclic <- function(data, pred) {
order_pred <- order(pred)
predvars <- names(pred)
ord_desc_vars <- names(pred[order_pred])
for (ordered_it in order_pred) {
descvar <- ord_desc_vars[ordered_it]
predecessor <- pred[descvar]
if (predecessor>0L) {
predvar <- predvars[predecessor]
# Predecessor variable being 0 (no survival, or no flowers) means
# that descendant variable cannot be observed:
data[is.na(data[[predvar]]) | data[[predvar]]==0L, descvar] <- NA
}
}
return(data)
}
# Note the names, important here:
varpred <- c(ld02=0, ld03=1, ld04=2, fl02=1, fl03=2, fl04=3,
hdct02=4, hdct03=5, hdct04=6)
NAechin <- asNA.acyclic(echinacea, pred=varpred)
yearvars <- c("ld02", "ld03", "ld04")
byyear <- reshape(NAechin, varying = list(yearvars), direction = "long",
timevar = "varld",times = as.factor(yearvars), v.names = "respld")
yearvars <- c("fl02", "fl03", "fl04")
flbyyear <- reshape(NAechin, varying = list(yearvars), direction = "long",
timevar = "varfl",times = as.factor(yearvars), v.names = "respfl")
yearvars <- c("hdct02", "hdct03", "hdct04")
hdctbyyear <- reshape(NAechin, varying = list(yearvars), direction = "long",
timevar = "varhdct",times = as.factor(yearvars), v.names = "resphdct")
byyear$varfl <- flbyyear$varfl
byyear$varhdct <- hdctbyyear$varhdct
byyear$respfl <- flbyyear$respfl
byyear$resphdct <- hdctbyyear$resphdct
byyear <- byyear[ ,c(12,4:6,10,13,14,11,15,16)]
head(byyear)
## ----asMM, size="small"---------------------------------------------------------------------------------------------------------------------------------------
asMM <- fitmv(submodels=list(
list(respld ~ varld + nsloc + ewloc+(1|id), family=binomial()),
list(respfl ~ varfl + nsloc + ewloc+(1|id), family=binomial()),
list(resphdct ~ varhdct + nsloc + ewloc+pop, family=Tpoisson())),
data=byyear)
if (requireNamespace("multcomp", quietly = TRUE)) {
library(multcomp)
#summary(glht(asMM,mcp("varld" = "Tukey"), coef.=fixef.HLfit)) # documented limitation
}
}
spaMM.options(spaMM_tol=spaMM_tol)
summary(warnings())
Try the spaMM package in your browser
Any scripts or data that you put into this service are public.
spaMM documentation built on Aug. 30, 2023, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
cat(crayon::yellow("\ntest-mv-extra:")) # not part of the testthat.R tests (neither test-composite-extra.R)
library(spaMM)
options(error=recover)
spaMM_tol <- local_tol <- spaMM.getOption("spaMM_tol")
local_tol$logL_tol <- 5e-07
spaMM.options(spaMM_tol=local_tol) # to control strictness of checks in independent-fit tests
{ # test missing data
set.seed(1)
d <- data.frame(y1 = rnorm(20), y2 = rnorm(20), x = runif(20))
d_y_missing <- d
d_y_missing[1, "y1"] <- NA
blaNA <- fitmv(submodels = list(one = list(y1 ~ x), two = list(y2 ~ x)),
# X2X=matrix(c(1,0,1,0,0,1,0,0,0,0,0,1),ncol=3,nrow=4,dimnames=list(NULL,c("(Intercept)","x_1","x_2"))),
data = d_y_missing)
testthat::expect_true(diff(range(predict(blaNA)[,1] - predict(blaNA, newdata=blaNA$data)[-1,1]))< 2e-16) # must work with X2X too
# predict(blaNA) does not predict for the line with NA response
# predict(blaNA, newdata=blaNA$data) predict for this line.
# It is important that the locvars provided to ..get_locdata() do not contain the response variable in the second case =>
# Cf .strip_cF_args(locformS[[mv_it]][-2]) in .calc_new_X_ZAC_mv()
testthat::expect_true(diff(range(predict(blaNA, blockSize=3)[,1] - predict(blaNA, newdata=blaNA$data, blockSize=3)[-1,1]))< 2e-16)
# there is slicing neither in first call (nrX = 0) nor in second (validrownames not NULL): cf conditions in predict.HLfit()
# With 'real' newdata without the validrownames attribute, slicing can occur:
testthat::expect_true(diff(range(predict(blaNA, newdata=d)[,1] - predict(blaNA, newdata=d, blockSize=3)[,1]))< 2e-16)
}
{
data("wafers")
me <- fitme(y ~ 1+(1|batch), family=Gamma(log), data=wafers)
set.seed(123)
y2 <- simulate(me, type="residual")
wafmv <- wafers
wafmv$batch2 <- wafmv$batch
wafmv$y2 <- y2
wafmv$ly <- log(wafmv$y)
wafmv$y3 <- log(y2)
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly ~ 1+(1|batch), family=gaussian()),
mod2=list(formula=y3 ~ 1+(1|batch2), family=gaussian())),
data=wafmv))
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14) ## there a resid.model so nothing is done with phi
simulate(zut1)
simulate(zut1, newdata=wafmv[1:3,])
cat(crayon::yellow("[ upper, get_ranPars() (-> VarCorr()) ] ; "))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly ~ 1+(1|batch), family=gaussian()),
mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian())),
data=wafmv, init=list(lambda=1), upper=list(lambda=0.02)))
testthat::expect_true(diff(range(get_ranPars(zut1, which="lambda"),0.02))<1e-12)
cat(crayon::yellow("(mv()|.) and (0+mv()|.); "))
(zut0 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(0+mv(1,2)|batch)),
mod2=list(formula=y3~X1+(0+mv(1,2)|batch), family=gaussian())),
data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(mv(1,2)|batch)),
mod2=list(formula=y3~X1+(mv(1,2)|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(mv(1,2)|batch)),
mod1=list(formula=ly~X1+(mv(1,2)|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut0),logLik(zut1),logLik(zut2)))<1e-05)
# # test equivalence of the two paremtrization, and that they are indeed distinguished by the code.
testthat::expect_true(diff(range(ranef(zut0)[[1]][,2]-rowSums(ranef(zut1)[[1]]))/(max(abs(ranef(zut0)[[1]]))))<1e-4)
(numinfo <- numInfo(zut0))
(numSEs <- sqrt(diag(solve(numinfo))))
condSEs <- summary(zut0,verbose=FALSE)$beta_table[,"Cond. SE"]
crit <- max(abs(numSEs[5:8]-condSEs))
testthat::test_that("numInfo() consistent with cond.SEs", testthat::expect_true(crit<1e-10))
{ # check of anova
zutLM <- fitmv(submodels=list(mod1=list(formula=ly~X1+batch),
mod2=list(formula=y3~1+batch, family=gaussian())),
data=wafmv)
crit <- max(abs(c(anova(fitme(ly~X1+batch, data=wafmv), method = "t.Chisq")[,3],
anova( fitme(y3~1+batch, data=wafmv), method = "t.Chisq")[,3]) -
anova(zutLM, type = "2")[,3]))
testthat::test_that("anova() OK on mv-LM",
testthat::expect_true(crit<1e-8))
}
#
anova(zut0, type = "2")
{ # missing data
afers <- wafmv
afers$y3[71:140] <- NA
afers$ly[1:70] <- NA
(zut0 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(0+mv(1,2)|batch)),
mod2=list(formula=y3~X1+(0+mv(1,2)|batch))),
data=afers))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(mv(1,2)|batch)),
mod2=list(formula=y3~X1+(mv(1,2)|batch))),
data=afers))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(mv(1,2)|batch)),
mod1=list(formula=ly~X1+(mv(1,2)|batch))),
data=afers))
logliks <- c(l0=logLik(zut0),l1=logLik(zut1),l2=logLik(zut2))
crit <- diff(range(logliks))
testthat::test_that(paste0("Hey, zut2 is still poor (inner estimation of phi?) logliks are", # so I made var_ranCoefs a reason for outer estim
paste(signif(logliks,6),collapse=",")),
testthat::expect_true(crit<1e-05)) # previous reasons for testing this were use_ZA_L or .calc_r22()
}
cat(crayon::yellow("ranCoefs; "))
(mod1 <- fitme(ly~X1, data=wafmv))
(mod2 <- fitme(y3~X1+(X2|batch), data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(X2|batch), family=gaussian()),
mod1=list(formula=ly~X1)),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-08)
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch)),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3~X1+(X2|batch), family=gaussian()),
mod1=list(formula=ly~X1+(X2|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-07)
cat(crayon::yellow("fixing ranCoefs in two ways: "))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch), fixed=list(ranCoefs=list("1"=c(0.02, -0.1, 0.005)))),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch)),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv, fixed=list(ranCoefs=list("1"=c(0.02, -0.1, 0.005)))))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-08)
cat(crayon::yellow("isDiagFamily in two ways: "))
(zut1 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch), fixed=list(ranCoefs=list("1"=c(NA, 0, NA)))),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=ly~X1+(X2|batch)),
mod2=list(formula=y3~X1+(X2|batch), family=gaussian())),
data=wafmv, fixed=list(ranCoefs=list("1"=c(NA, 0, NA)))))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-08)
testthat::expect_true(unique(length(attr(zut1, "optimInfo")$optim.pars$trRanCoefs[["1"]]),
length(attr(zut2, "optimInfo")$optim.pars$trRanCoefs[["1"]]))==2L)
cat(crayon::yellow("with resid.model; "))
# independent-fit test
(mod1 <- fitme(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2), data=wafmv))
(mod2 <- fitme(formula=y2 ~ 1+(1|batch2), family=Gamma(log), data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2)),
mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(mod1)+logLik(mod2)))<2e-7)
{ # check of anova
mod1 <- fitme(formula=y ~ X1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2), data=wafmv)
mod2 <- fitme(formula=y2 ~ X1+(1|batch2), family=Gamma(log), data=wafmv)
zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ X1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2)),
mod2=list(formula=y2 ~ X1+(1|batch2), family=Gamma(log))),
data=wafmv)
crit <- max(abs(c(anova(mod1, method = "t.Chisq")[,3],
anova(mod2, method = "t.Chisq")[,3]) -
anova(zut1, type = "2")[,3]))
testthat::test_that("anova() OK on mv-GLMM",
testthat::expect_true(crit<1e-6))
}
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ X3+I(X3^2))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
cat(crayon::yellow("full dhglm (the 5 'phi' message are from different fits); "))
(mod1 <- fitme(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch), data=wafmv))
(mod2 <- fitme(formula=y2 ~ 1+(1|batch2), family=Gamma(log), data=wafmv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch)),
mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log))),
data=wafmv))
get_residVar(zut1, newdata=zut1$data) # with a phiHGLM: residVar extractor has ignored this case for a long time.
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-4)
## To 'cover' the init procedure of refits in fitmv_body when there is a mixed-effect phi-resid.model:
# Here one must force some outer optimization for the mean-response model hence the lambda init:
cover_residM_reinit <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch2), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch))),
data=wafmv, init=list(lambda=c(0.01,0.01)))
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log),resid.model= ~ 1+(1|batch))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-4)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14) ## there a resid.model so nothin is done with phi
spaMM_boot(zut1, function(v) var(v), nsim=3L, type ="marginal")$bootreps
confint(zut1,"(Intercept)_1")
cat(crayon::yellow("fixing phi: four different ways; "))#
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.001)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.002))),
data=wafmv, fixed=list(lambda=0.1)))
# Check that <processed>$phi.Fix was correctly set up otherwise phi might still be fixed in fit but not considered as such in post-fit)
testthat::expect_true(all(sapply(zut1$phi.object, function(v) attr(v$phi_outer,"type"))==rep("fix",2)))
(zut2 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("1"=0.001,"2"=0.002), lambda=0.1)))
(zut3 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list(0.001,0.002), lambda=0.1)))
(zut4 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.001)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("2"=0.002), lambda=0.1)))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2),logLik(zut3),logLik(zut4)))<1e-14)
cat(crayon::yellow("visual checks; "))# and low fixed phis are useful for the following visual checks:
set.seed(123)
ressim1 <- simulate(zut1, type="residual")
ressim2 <- simulate(zut1, newdata=zut1$data, type="residual")
margsim1 <- simulate(zut1)
margsim2 <- simulate(zut1, newdata=zut1$data)
plot(ressim1,ressim2);abline(0,1) # clusters on the diagonal
plot(margsim1,margsim2);abline(0,1) # TWO indep draws of the 11 levels => clusters are not on the diagonal
# the two submodels share the ranef values so the clusters are always on the diagonal, even for margsim:
plot(ressim2[1:198], ressim2[1L+(1:198)]);abline(0,1) # on diagonal
plot(margsim2[1:198], margsim2[1L+(1:198)]);abline(0,1) # on diagonal
cat(crayon::yellow("confint with 'fixed'; "))# Checking that confint() obeys fixed values fixed in different ways
(zute <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(lambda=0.1)))
confint(zute,"(Intercept)_1")$lowerfit # estimated phis
(zutf <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.2224)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.2103))),
data=wafmv, fixed=list(lambda=0.1)))
confint(zutf,"(Intercept)_1")$lowerfit # fixed phis
(zutfm <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log), fixed=list(phi=0.2224)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("2"=0.2103), lambda=0.1)))
confint(zutfm,"(Intercept)_1")$lowerfit # fixed phis
(zutf2 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y2 ~ 1+(1|batch), family=Gamma(log))),
data=wafmv, fixed=list(phi=list("2"=0.2103), lambda=0.1)))
confint(zutf2,"(Intercept)_1")$lowerfit # fixed 2nd phi
cat(crayon::yellow("some of the early tests; "))# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian(log))),
data=wafmv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian(log)),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log))),
data=wafmv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-8)
# permutation test. The init from hazards of development
(zut1 <- fitmv(submodels=list(mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian()),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log))),
control.HLfit=list(LevenbergM=TRUE),
data=wafmv,init=list(phi=list("1"=1 ))))
(zut2 <- fitmv(submodels=list(mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log)),
mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian())),
control.HLfit=list(LevenbergM=TRUE),
data=wafmv,init=list(phi=list("2"=1 ))))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
(zut <- fitmv(submodels=list(mod2=list(formula=y3 ~ 1+(1|batch), family=gaussian()),
mod1=list(formula=y ~ 1+(1|batch), family=Gamma(log))),
data=wafmv,init=list(phi=list("1"=1,"2"=0.2 ))))
attr(zut,"optimInfo")$LUarglist$canon.init
}
{
data("clinics")
#
(fitClinics <- HLfit(cbind(npos,nneg)~treatment+(1|clinic),family=binomial(),data=clinics))
set.seed(123)
y2 <- simulate(fitClinics, type="residual")
climv <- clinics
climv$np2 <- y2
climv$nn2 <- climv$npos + climv$nneg - y2
cat(crayon::yellow("binomial-poisson; "))# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(1|clinic),family=poisson())),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(1|clinic),family=poisson()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
cat(crayon::yellow("predVar cov gaussian-poisson; "))
# independent-fits test
(fg <- fitme(formula=npos~treatment+(1|clinic),family=gaussian(),data=climv))
(fp <- fitme(formula=np2~treatment+(1|clinic),family=poisson(),data=climv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=npos~treatment+(1|clinic),family=gaussian()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(+1|clinic),family=poisson()),
mod1=list(formula=npos~treatment+(1|clinic),family=gaussian())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2),logLik(fg)+logLik(fp)))<3e-7)
testthat::expect_true(diff(range(get_predVar(zut1)-c(get_predVar(fg),get_predVar(fp))))<1e-5)
testthat::expect_true(diff(range(get_predVar(zut2)-c(get_predVar(fp),get_predVar(fg))))<1e-5)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data,variances=list(cov=TRUE)) -
get_predVar(zut1,variances=list(cov=TRUE))))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data,variances=list(cov=TRUE)) -
Matrix::bdiag(get_predVar(fg,variances=list(cov=TRUE)),
get_predVar(fp,variances=list(cov=TRUE)))))<1e-5)
get_intervals(zut1)
get_intervals(zut1, intervals = "respVar")
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(formula=npos~treatment+(1|clinic),family=gaussian()),
mod2=list(formula=np2~treatment+(1|clinic),family=poisson())),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(1|clinic),family=poisson()),
mod1=list(formula=npos~treatment+(1|clinic),family=gaussian())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2)))<1e-10)
testthat::expect_true(diff(range(get_predVar(zut1)-get_predVar(zut2)[outer((1:16),c(16,0),"+")]))<1e-10)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data)-get_predVar(zut1)))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut1, newdata=zut1$data,variances=list(cov=TRUE)) -
get_predVar(zut1,variances=list(cov=TRUE))))<1e-14)
testthat::expect_true(diff(range( get_predVar(zut2, newdata=zut2$data)-get_predVar(zut2)))<1e-14)
}
{
meth <- "ML(1,1)" # REML seems to work but the tests are OK only with lower accuracy
# independent-fit test
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv, method=meth))
(fp <- fitme(formula=np2~treatment+(1|clinic),family=poisson(), data=climv, method=meth))
(zut1 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
data=climv, method=meth))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(+1|clinic),family=poisson()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv, method=meth))
# test simple init
(zut3 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
init=list(lambda=c("1"=1.1,"2"=2.2)), data=climv, method=meth))
testthat::expect_true(identical(attr(zut3,"optimInfo")$LUarglist$canon.init, list(lambda=c("1"=1.1,"2"=2.2)))) # inits heeded
# test fancy names in init
zut4 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson())),
init=list(lambda=c("clinic.1"=2.2,"clinic"=1.1)), data=climv, method=meth)
testthat::expect_true(identical(attr(zut4,"optimInfo")$LUarglist$canon.init, list(lambda=c("1"=1.1,"2"=2.2)))) # inits heeded
testthat::expect_true(diff(range(logLik(zut1),logLik(zut2), logLik(zut3), logLik(zut4), logLik(fb)+logLik(fp)))<2e-6) # 1e-5 for REML # stricter conv check decreased the precision of the comparison.
pfb <- get_predVar(fb)
pfp <- get_predVar(fp)
pzut1 <- get_predVar(zut1)
pzut1n <- get_predVar(zut1, newdata = zut1$data)
testthat::expect_true(diff(range(c(pfb,pfp)-pzut1, pzut1n-pzut1))<1e-5)
pzut2 <- get_predVar(zut2)
pzut2n <- get_predVar(zut2, newdata = zut2$data)
testthat::expect_true(diff(range(c(pfp,pfb)-pzut2, pzut2n-pzut2))<1e-5)
cat(crayon::yellow("deliberate warnings; ")) # (Deliberately generating warnings:)
oldopt <- options(warn=0L)
(zut5 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), init=list(lambda=1.1)),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson(), init=list(lambda=2.2))),
data=climv, method=meth))
options(oldopt)
testthat::expect_true(identical(attr(zut5,"optimInfo")$LUarglist$canon.init, NULL)) # inits NOT heeded as init is not preprocessed
cat(crayon::yellow("rand family and many post-fit fns; "))# rand.family independent-fit test
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(+1|clinic),family=poisson(), rand.family=Gamma(log)),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut1 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=np2~treatment+(+1|clinic),family=poisson(), rand.family=Gamma(log))),
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv))
(fg <- fitme(formula=np2~treatment+(1|clinic),family=poisson(), data=climv, rand.family=Gamma(log)))
testthat::expect_true(diff(range(logLik(zut1), logLik(fb)+logLik(fg)))<5e-6) # stricter conv check decreased the precision of the comparison.
cat(crayon::yellow("some extractors; "))
testthat::expect_true(diff(range(residuals(zut1)-c(residuals(fb),residuals(fg))))<1e-5)
formula(zut1) # list
terms(zut1) # list
nobs(zut1) # single value # decision made
# rand.family permutation test (rand.family need to be specified in each submodel)
(zut1 <- fitmv(submodels=list(mod1=list(cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), rand.family=Gamma(log)),
mod2=list(np2~treatment+(1|clinic),family=poisson(), rand.family=Gamma(log))),
data=climv))
(zut2 <- fitmv(submodels=list(mod2=list(formula=np2~treatment+(1|clinic),family=poisson(), rand.family=Gamma(log)),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), rand.family=Gamma(log))),
data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-9)
update_formulas(zut1,formula(zut1))
update(zut1, formula.=formula(zut1))
update(zut1, formula.=list(cbind(npos,nneg)~1+(1|clinic),
np2~1+(1|clinic)))
# step(zut1) dos not stop() but does not steps...
simulate(zut1,nsim=3) # checks that mv simulate Tpoisson works (possibly depending on same attributes as Tnegbin)
update_resp(zut1, newresp=simulate(zut1))
# CI for the variance of the random effect:
( ci <- confint(zut1,parm=function(fit){VarCorr(fit)[1,"Variance"]},
boot_args=list(nb_cores=9, nsim=9, seed=123)) )
# The distribution of bootstrap replicates:
plot(ecdf(ci$call$t))
zutnull <- update(zut1, submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), rand.family=Gamma(log)),
mod2=list(formula=np2~1+(1|clinic),family=poisson(), rand.family=Gamma(log))))
anova(zut1,zutnull)
anova(zut1,zutnull, boot.repl=99, nb_cores=9)
anova(zut1,zutnull, boot.repl=99, nb_cores=9)
# MSFDR(zutnull,zut1) # documented problem
cat(crayon::yellow("Tpoisson; "))# Tpoisson independent-fit test
(zut <- fitmv(submodels=list(mod2=list(formula=I(1L+np2)~treatment+(+1|clinic),family=Tpoisson()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=I(1L+np2)~treatment+(+1|clinic),family=Tpoisson())),
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv))
(fTp <- fitme(formula=I(1L+np2)~treatment+(1|clinic),family=Tpoisson(), data=climv))
testthat::expect_true(diff(range(logLik(zut), logLik(fb)+logLik(fTp)))<2e-6) # stricter conv check decreased the precision of the comparison.
## negbin(): outer-optimized dispersion parameters.
meth <- "ML(1,1)"
{ # quite slow 2*40s
# independent-fit test without fixed effects
(zut1 <- fitmv(submodels=list(mod2=list(formula=I(20*np2)~0+(+1|clinic),family=negbin()),
mod1=list(formula=cbind(npos,nneg)~0+(1|clinic),family=binomial())),
method=meth,
data=climv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~0+(1|clinic),family=binomial()),
mod2=list(formula=I(20*np2)~0+(+1|clinic),family=negbin())),
method=meth,
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~0+(1|clinic),family=binomial(),
method=meth,
data=climv))
(fn <- fitme(formula=I(20*np2)~0+(+1|clinic),family=negbin(),
method=meth,
data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(fb)+logLik(fn)))<1e-06)
}
cat(crayon::yellow("Tnegbin; "))## independent-fit test Tnegbin; outer-optimized dispersion parameters
(zut1 <- fitmv(submodels=list(mod2=list(formula=I(1L+20*np2)~treatment+(+1|clinic),family=Tnegbin()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=I(1L+20*np2)~treatment+(+1|clinic),family=Tnegbin())),
data=climv))
(fb <- fitme(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial(), data=climv))
(fTn <- fitme(formula=I(1L+20*np2)~treatment+(1|clinic),family=Tnegbin(), data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(fb)+logLik(fTn)))<1e-06) # This did not work at first by obsInfo bc
# there were convergence problems for IRLS but LevM crit was not true => change in crit using pot4improv.
## permutation test Tnegbin; outer-optimized dispersion parameters
(zut1 <- fitmv(submodels=list(mod2=list(formula=I(1L+20*np2)~treatment+(1|clinic),family=Tnegbin()),
mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial())),
data=climv))
(zut2 <- fitmv(submodels=list(mod1=list(formula=cbind(npos,nneg)~treatment+(1|clinic),family=binomial()),
mod2=list(formula=I(1L+20*np2)~treatment+(1|clinic),family=Tnegbin())),
data=climv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-06)
simulate(zut1,nsim=3) # checks that mv simulate Tnegbin keeps required attributes
{
cat(crayon::yellow("COMPoisson; "))## independent-fit test; outer-optimized dispersion parameters
data("freight") ## example from Sellers & Shmueli, Ann. Appl. Stat. 4: 943–961 (2010)
(mod1 <- fitme(broken ~ transfers+(1|id), data=freight, family = COMPoisson()))
freimv <- freight
set.seed(123)
freimv$brok2 <- simulate(mod1)
(mod2 <- fitme(brok2 ~ transfers,family=poisson(), data=freimv))
(zut1 <- fitmv(submodels=list(mod1=list(broken ~ transfers+(1|id),family=COMPoisson()),
mod2=list(brok2 ~ transfers,family=poisson())),
data=freimv))
(zut2 <- fitmv(submodels=list(mod2=list(brok2 ~ transfers,family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-06)
## permutation test COMPoisson; outer-optimized dispersion parameters
(zut1 <- fitmv(submodels=list(mod1=list(broken ~ transfers+(1|id),family=COMPoisson()),
mod2=list(brok2 ~ transfers+(1|id),family=poisson())),
data=freimv))
(zut2 <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-6)
simulate(zut1,nsim=3) # checks that mv simulate COMPoisson
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=COMPoisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("1"=0.5,"2"=1.1)))) # names needed !
attr(zutxx,"optimInfo")$LUarglist$init.optim$COMP_nu
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=COMPoisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
attr(zutxx,"optimInfo")$LUarglist$init.optim$COMP_nu
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
attr(zutxx,"optimInfo")$LUarglist$init.optim$COMP_nu
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson(1.1))),
data=freimv))
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=poisson()),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson(1.1))),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
(zutxx <- fitmv(submodels=list(mod2=list(brok2 ~ transfers+(1|id),family=COMPoisson(1.1)),
mod1=list(broken ~ transfers+(1|id),family=COMPoisson())),
data=freimv, lower=list(COMP_nu=c("2"=1.1))))
## assign("last.warning", NULL, envir = baseenv()) # flush the .COMP_maxn() warnings()
}
}
data("Loaloa")
lll <- Loaloa
lll$ID <- lll$ID2 <- 1+(seq(nrow(lll)) %% 2)
lll$resp <- 1+lll$npos + 10*lll$ID
(tnb <- fitme(resp~1+(1|ID), data=lll,family=Tnegbin()))
set.seed(124)
lll$r1 <- simulate(tnb)
lll$r2 <- simulate(tnb)
lll$long2 <- lll$longitude
## independent-fit test 2 Tnegbin
(tnb1 <- fitme(r1~1+(1|ID), data=lll,family=Tnegbin()))
(tnb2 <- fitme(r2~1+(1|ID), data=lll,family=Tnegbin()))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+(1|ID),family=Tnegbin()),
mod2=list(r2~1+(1|ID2),family=Tnegbin())),
data=lll))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+(1|ID2),family=Tnegbin()),
mod1=list(r1~1+(1|ID),family=Tnegbin())),
data=lll))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tnb1)+logLik(tnb2)))<1e-06)
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+(1|ID),family=Tnegbin()),
mod2=list(r2~1+(1|ID),family=Tnegbin())),
data=lll))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+(1|ID),family=Tnegbin()),
mod1=list(r1~1+(1|ID),family=Tnegbin())),
data=lll))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-8)
if (FALSE) { # slowish ~4s
cat(crayon::yellow("Matern; "))## independent-fit test with a Matern
(tp1 <- fitme(r1~1+Matern(1|longitude+latitude), data=lll,family=poisson()))
(tp2 <- fitme(r2~1+(1|ID), data=lll,family=poisson()))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson()),
mod2=list(r2~1+(1|ID2),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+(1|ID2),family=poisson()),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tp1)+logLik(tp2)))<1e-08)
}
if (FALSE) { # slow! > 4mn total
## independent-fit test with two Matern
(tp1 <- fitme(r1~1+Matern(1|longitude+latitude), data=lll,family=poisson()))
(tp2 <- fitme(r2~1+Matern(1|long2+latitude), data=lll,family=poisson()))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson()),
mod2=list(r2~1+Matern(1|long2+latitude),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+Matern(1|long2+latitude),family=poisson()),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, verbose=c(TRACE=TRUE)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tp1)+logLik(tp2)))<1e-08)
}
if (FALSE) { # slow! > 150s total
## independent-fit test with two Matern and fixed nu's
(tp1 <- fitme(r1~1+Matern(1|longitude+latitude), data=lll,family=poisson(), fixed=list(nu=0.5)))
(tp2 <- fitme(r2~1+Matern(1|long2+latitude), data=lll,family=poisson(), fixed=list(nu=1)))
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=0.5)),
mod2=list(r2~1+Matern(1|long2+latitude),family=poisson(), fixed=list(nu=1))),
data=lll, verbose=c(TRACE=TRUE)))
(zut2 <- fitmv(submodels=list(mod2=list(r2~1+Matern(1|long2+latitude),family=poisson(), fixed=list(nu=1)),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=0.5))),
data=lll, verbose=c(TRACE=TRUE)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(tp1)+logLik(tp2)))<1e-08)
}
if (FALSE) { # slow 22s
# permutation test
(zut1 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=1)),
mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, verbose=c(TRACE=1)))
(zut2 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson()),
mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(nu=1))),
data=lll, verbose=c(TRACE=1)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
predict(zut1)
simulate(zut1)
}
# default-name test
cat(crayon::yellow("'fixed' both in submodel and global call; "))
(zut1 <- fitmv(submodels=list(mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson()),
mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(rho=1,nu=1))),
data=lll, fixed=list(lambda=c("1"=666)), verbose=c(TRACE=1)))
(zut2 <- fitmv(submodels=list(mod1=list(r1~1+Matern(1|longitude+latitude),family=poisson(), fixed=list(rho=1,nu=1)),
mod2=list(r2~1+Matern(1|longitude+latitude),family=poisson())),
data=lll, fixed=list(lambda=c(666)), verbose=c(TRACE=1)))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
{
data("blackcap")
MLdistMat2 <- as.matrix(proxy::dist(blackcap[,c("latitude","longitude")]))
MLcorMat2 <- MaternCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
nu=0.6285603,rho=0.0544659)
cap_mv <- blackcap
cap_mv$name <- as.factor(rownames(blackcap))
cap_mv$grp <- 1L+(blackcap$migStatus>1)
set.seed(123)
cap_mv$status2 <- blackcap$migStatus+ rnorm(14,sd=0.001)
cat(crayon::yellow("corrMatrix vs Matern; (map_ranef too)"))# corrMatrix
## independent-fit test with compar to equivalent Matern:
# need to fix this phi to avoid logLik uncertainty for low phi:
(mod1 <- fitme(migStatus ~ 1+ corrMatrix(1|name),data=cap_mv, corrMatrix=MLcorMat2, fixed=list(phi=0.1)))
(mod2 <- fitme(status2 ~ 1+ (1|grp),data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut1b <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1), corrMatrix=MLcorMat2),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
(zut1c <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, covStruct=list(corrMatrix=MLcorMat2,NULL,corrMatrix=NULL))) # a bit of abuse...
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ (1|grp)),
mod1=list(migStatus ~ 1+ corrMatrix(1|name),fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, fixed=list(rho=0.0544659,nu=0.6285603)))
(zut4 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
map_ranef(zut4)
simulate(zut1, newdata=cap_mv[1:3,])
simulate(zut1b, newdata=cap_mv[1:3,])
simulate(zut1c, newdata=cap_mv[1:3,])
simulate(zut2, newdata=cap_mv[1:3,])
simulate(zut3, newdata=cap_mv[1:3,])
simulate(zut4, newdata=cap_mv[1:3,])
get_predVar(zut4, variances=list(cov=TRUE))
if (FALSE) {
get_predVar(zut4, variances=list(cov=TRUE)) -> bla
bli <- bla
dim(bli) <- c(14,2,14,2)
bli[,1,,2]
bli[,1,,1]
bli[,2,,2]
}
testthat::expect_true(diff(range(logLik(zut1), logLik(zut1b), logLik(zut1c),
logLik(zut2), logLik(zut3), logLik(zut4), logLik(mod1)+logLik(mod2)))<1e-10)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
## permutation test with compar to equivalent Matern:
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ corrMatrix(1|name))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ corrMatrix(1|name)),
mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat2))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603)),
mod2=list(status2 ~ 1+ Matern(1|longitude+latitude))),
data=cap_mv))
(zut4 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ Matern(1|longitude+latitude)),
mod1=list(migStatus ~ 1+ Matern(1|longitude+latitude),fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(zut3), logLik(zut4)))<1e-9)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
get_predVar(zut4, variances=list(cov=TRUE))
if (FALSE) { # may be used to check the mapping, but to have asym off-diag blocks, use migStatus ~ 1+ Matern(1|longitude+latitude)+(1|grp), the latter with fixed large variance
get_predVar(zut4, variances=list(cov=TRUE)) -> bli
nresp <- length(attr(bli, "respnames")) # but I commented out code providing this odd attribute, see get_predVar() and .predict_body()
nobs <- nrow(bli)/nresp
dim(bli) <- c(nobs,nresp,nobs,nresp)
dimnames(bli)[[2]] <- dimnames(bli)[[4]] <- attr(bli, "respnames")
bli[,1,,2]
bli[,1,,1]
bli[,2,,2]
bli[1,,1,]
}
cat(crayon::yellow("distMatrix; "))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Matern(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Matern(1|name))),
distMatrix=MLdistMat2, fixed=list(rho=0.0544659,nu=0.6285603),
data=cap_mv))
(zut2 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Matern(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Matern(1|name))),
distMatrix=2*MLdistMat2, fixed=list(rho=0.0544659/2,nu=0.6285603),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-9)
if (FALSE) { # Amusingly one can fit the same data by two sub-models (even tow diferent ones).
# we have twice the conditional likelihood, and once the ranef lik, so (as in any mv model with shared ranefs) marginal likelihoods do not add up.
# The phi's are low and they lower value must be controlled in order to make detailed numerical comparisons
(mod1 <- fitme(migStatus ~ 1+ corrMatrix(1|name),data=cap_mv, corrMatrix=MLcorMat))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name)),
mod2=list(migStatus ~ 1+ corrMatrix(1|name))),
data=cap_mv, corrMatrix=MLcorMat2))
vcov(zut1) # with phi->0 and a common ranef, this may be logical
}
cat(crayon::yellow("corrMatrix vs Cauchy; "))
MLcorMat3 <- CauchyCorr(proxy::dist(blackcap[,c("latitude","longitude")]),
shape=1,longdep=0.5) # and default rho=1!
## independent-fit test with compar to equivalent Matern:
# need to fix this phi to avoid logLik uncertainty for low phi:
(mod1 <- fitme(migStatus ~ 1+ corrMatrix(1|name),data=cap_mv, corrMatrix=MLcorMat3, fixed=list(phi=0.1)))
(mod2 <- fitme(status2 ~ 1+ (1|grp),data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut1b <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1), corrMatrix=MLcorMat3),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
(zut1c <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, covStruct=list(corrMatrix=MLcorMat3,NULL,corrMatrix=NULL))) # a bit of abuse...
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ (1|grp)),
mod1=list(migStatus ~ 1+ corrMatrix(1|name),fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv, fixed=list(shape=1,rho=1,longdep=0.5)))
(zut4 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude), fixed=list(phi=0.1,rho=1,shape=1,longdep=0.5)),
mod2=list(status2 ~ 1+ (1|grp))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut1b), logLik(zut1c),
logLik(zut2), logLik(zut3), logLik(zut4), logLik(mod1)+logLik(mod2)))<1e-10)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
## permutation test with compar to equivalent Cauchy:
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ corrMatrix(1|name))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ corrMatrix(1|name)),
mod1=list(migStatus ~ 1+corrMatrix(1|name), fixed=list(phi=0.1))),
data=cap_mv, corrMatrix=MLcorMat3))
(zut3 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude), fixed=list(phi=0.1,rho=1,shape=1,longdep=0.5)),
mod2=list(status2 ~ 1+ Cauchy(1|longitude+latitude))),
data=cap_mv))
(zut4 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ Cauchy(1|longitude+latitude)),
mod1=list(migStatus ~ 1+ Cauchy(1|longitude+latitude),fixed=list(phi=0.1,rho=1,shape=1,longdep=0.5))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(zut3), logLik(zut4)))<1e-8)
testthat::expect_true(diff(range( predict(zut1, newdata=zut1$data)-predict(zut1)))<1e-14)
cat(crayon::yellow("distMatrix with Cauchy; "))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Cauchy(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Cauchy(1|name))),
distMatrix=MLdistMat2, fixed=list(rho=1,shape=1,longdep=0.5),
data=cap_mv))
(zut2 <- fitmv(submodels=list(mod1=list(migStatus ~ 1+Cauchy(1|name), fixed=list(phi=0.1)),
mod2=list(status2 ~ 1+ Cauchy(1|name))),
distMatrix=2*MLdistMat2, fixed=list(rho=1/2,shape=1,longdep=0.5),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-9)
}
{cat(crayon::yellow("IMRF; "))# fit IMRF
# There is a check of simulate of a large IMRF + Matern hurdle model
{ # create IMRF model
## Creating the mesh
oldMDCopt <- options(Matrix.warnDeprecatedCoerce = 0) # INLA issue
mesh <- INLA::inla.mesh.2d(loc = blackcap[, c("longitude", "latitude")],
cutoff=30,
max.edge = c(3, 20))
mesh$n ## 40
matern <- INLA::inla.spde2.matern(mesh)
options(oldMDCopt)
}
{ # Aim was first to check mv IMRF predictions, but putting a (G)LM as first model caught many problems
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
fixed=list(phi=c(0.02,0.02)),
data=cap_mv))
p1 <- predict(zut1)
p2 <- predict(zut1, newdata=zut1$data) # note difference in frame attribute. What do we want? _FIXME_
testthat::expect_true(diff(range(p2-p1))<1e-14)
}
{ # independent-fit test
# a bit tricky bc nloptr does not find the univariate optimum without some help (... tiny blackcap data)
# and in mv case, it finds it by default or not depending on order...
(mod1 <- fitme(migStatus ~ means + IMRF(1|longitude+latitude, model=matern), init=list(lambda=0.02),
data=cap_mv, verbose=c(TRACE=FALSE)))
(mod2 <- fitme(status2 ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1,rho=0.0544659,nu=0.6285603), data=cap_mv,
verbose=c(TRACE=FALSE)))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern)),
mod2=list(status2 ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1))),
data=cap_mv, fixed=list(rho=0.0544659,nu=0.6285603), init=list(lambda=c(0.02,NA)))) # Matern in spprec...
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ Matern(1|longitude+latitude), fixed=list(phi=0.1)),
mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern))),
data=cap_mv, fixed=list(rho=0.0544659,nu=0.6285603), init=list(lambda=c(NA,0.02)))) # Matern in spprec...
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-08)
}
{ # permutation test
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern)),
mod2=list(status2 ~ means+ IMRF(1|longitude+latitude, model=matern))),
data=cap_mv))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ means+ IMRF(1|longitude+latitude, model=matern)),
mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern))),
data=cap_mv)) # sensitive to .solve_crossr22( ., use_crossr22=TRUE)
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-07)
# The following with status2 ~ 1+ ... rather than status2 ~ means+ ... makes lambda diverge (jointly with kappa: does not occur for fixed kappa).
# Then all comparisons are sensitive to floating point precision (and no stringent permutation test is passed).
# but let us keep the zut2_testr22 computation as it has been useful in the past to 'torture-test' .calc_r22() algorithms
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern)),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
init=list(lambda=c(0.02),phi=c(0.02,0.02)),
data=cap_mv))
crit <- diff(range(logLik(zut1), -6.48825776602))
try(testthat::test_that(paste0("criterion was ",signif(crit,6)," from -6.48825776602"), # affected by use_ZA_L or .calc_r22() ... and minKappa...
testthat::expect_true(crit<1e-08)))
(zut2_testr22 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern)),
mod1=list(migStatus ~ means + IMRF(1|longitude+latitude, model=matern))),
init=list(lambda=c(0.02),phi=c(0.02,0.02)),
data=cap_mv)) # has been sensitive to .solve_crossr22( ., use_crossr22=TRUE)
crit <- diff(range(-6.48825776601 , logLik(zut2_testr22)))
try(testthat::test_that(paste0("criterion was ",signif(crit,6)," from -6.48830317593"), # affected by use_ZA_L or .calc_r22() ... and minKappa...
testthat::expect_true(crit<1e-08)))
}
if (spaMM.getOption("example_maxtime")>87) { cat(crayon::yellow("multIMRF indep-fit tests; "))
(mrf1fixx <- fitme(migStatus ~ 1 + (1|pos) +
multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(phi=1,lambda=c("1"=0.5),
hyper=list("1"=list(hy_kap=0.1,hy_lam=1)))) )
(mrf2 <- fitme(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern),
data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + (1|pos) + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
fixed=list(phi=1,lambda=c("1"=0.5),
hyper=list("1"=list(hy_kap=0.1,hy_lam=1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fixx)+logLik(mrf2)))<1e-7) # precision modified (effect of changing bobyqa's rhoend?)
#
(mrf1fix <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(hyper=list("1"=list(hy_kap=0.1,hy_lam=1)))) )
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
fixed=list(hyper=list("1"=list(hy_kap=0.1,hy_lam=1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fix)+logLik(mrf2)))<1e-10)
(mrf1fix <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(hyper=list("1"=list(hy_kap=0.1)))) )
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
fixed=list(hyper=list("1"=list(hy_kap=0.1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fix)+logLik(mrf2)))<1e-10)
(mrf1fix <- fitme(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2),
data=blackcap, fixed=list(hyper=list("1"=list(hy_lam=1)))) )
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
fixed=list(hyper=list("1"=list(hy_lam=1))),
data=cap_mv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mrf1fix)+logLik(mrf2)))<1e-10)
{ # particularly slow
{
# now there is a global maximum at phi=0 which may be missed when phi outer optim is used,
# But not missed when default augZXy is used. augZXy is used by default for mrf1, but not in mv => discrepancy
spaMM.options(allow_augZXy=FALSE)
logLik(mrf1F <- fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv)) # bad
logLik(fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv, init=list(phi=(1e-4)))) # good
spaMM.options(allow_augZXy=NULL) #
# This may lose accuracy when .calc_r22() is modified! :
logLik(mrf1T <- fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv)) # good; use y-augmented matrix
logLik(mrf1I <- fitme(migStatus ~ 1 +multIMRF(1|longitude+latitude,margin=5,levels=2),data=cap_mv, init.HLfit=list(phi=1))) # good; use y-augmented matrix !
# so to compare with mv fit we must provide a low init phi to the submodel
}
(mrf2 <- fitme(status2 ~ 1, data=cap_mv))
(zut0 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1)),
data=cap_mv, init=list(phi=list("1"=1e-4))))
testthat::expect_true(diff(range(logLik(zut0), logLik(mrf1T)+logLik(mrf2)))<1e-5)
(mrf3 <- fitme(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern), data=cap_mv))
# inits to try to speed the fit () (and lambda=c("3"=.) works as it should)
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern))),
data=cap_mv, init=list(lambda=c("3"=5),phi=list("1"=1e-4,"2"=0.05)),verbose=c(TRACE=0.5)))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ IMRF(1|longitude+latitude, model=matern)),
mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(lambda=c("1"=5),phi=list("1"=0.05,"2"=1e-4))))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mrf1T)+logLik(mrf3)))<1e-5)
}
} else cat(crayon::bgGreen("\n multIMRF indep-fit tests are slow (~87s). Run them once in a while; "))
{ cat(crayon::yellow("multIMRF permutation tests; ")) # ~ 13s
# reason for init phi as above: to avoid a local maximum
(zut1 <- fitmv(submodels=list(mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod2=list(status2 ~ 1+ multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(phi=list("1"=1e-4,"2"=1e-4))))
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ multIMRF(1|longitude+latitude,margin=5,levels=2)),
mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(phi=list("1"=1e-4,"2"=1e-4))))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
if (FALSE) {
# The following code fits a single (repeated) kappa value, and only one lambda hyperparam.
# So the single IMRF from multIMRF(. ,levels=1) is possibly
# recognized as first level of multIMRF(. ,levels=2). Which may be nice, or may not be the intent. (__FIXME__).
# (What does .calc_normalized_ZAlist() do ?)
# Adding a fictitious argument bla=666 has no effect as it is not retained in the expanded formula.
(zut2 <- fitmv(submodels=list(mod2=list(status2 ~ 1+ multIMRF(1|longitude+latitude,margin=5,levels=1)),
mod1=list(migStatus ~ 1 + multIMRF(1|longitude+latitude,margin=5,levels=2))),
data=cap_mv, init=list(phi=list("1"=1e-4,"2"=1e-4))))
}
}
}
{ cat(crayon::yellow("adjacency; "))
data("scotlip")
(mod1 <- fitme(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)),
adjMatrix=Nmatrix, family=poisson(), data=scotlip) )
scotmv <- scotlip
set.seed(123)
scotmv$cases2 <- simulate(mod1)
scotmv$code2 <- scotmv$gridcode
(mod2 <- fitme(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)),
adjMatrix=Nmatrix, family=poisson(), data=scotmv) )
{ # independent-fit test
# a bit tricky bc nloptr does not find the univariate optimum without some help (... tiny blackcap data)
# and in mv case, it finds it by default or not depending on order...
(zut1 <- fitmv(submodels=list(mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson()),
mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|code2)+offset(log(expec)), family=poisson())),
data=scotmv,covStruct=list(adjMatrix=Nmatrix,adjMatrix=Nmatrix))) # Matern in spprec...
(zut2 <- fitmv(submodels=list(mod1=list(cases2 ~ I(prop.ag/10)+adjacency(1|code2)+offset(log(expec)), family=poisson()),
mod2=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson())),
data=scotmv,covStruct=list(adjMatrix=Nmatrix,adjMatrix=Nmatrix))) # Matern in spprec...
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2), logLik(mod1)+logLik(mod2)))<1e-04) # __F I X M E__ was more accurate by outer lambda
}
{ # permutation test
# here to the order affects nloptr... tiny blackcap data again. Note lambda divergence.
# we can avoid that by either providing fixed or init phi (the latter slow)
(zut1 <- fitmv(submodels=list(mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson()),
mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson())),
data=scotmv,adjMatrix=Nmatrix))
(zut2 <- fitmv(submodels=list(mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson()),
mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)), family=poisson())),
data=scotmv,adjMatrix=Nmatrix))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
#
(zut1 <- fitmv(submodels=list(mod1=list(cases ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec))),
mod2=list(cases2 ~ I(prop.ag/10)+adjacency(1|gridcode)+offset(log(expec)))),
data=scotmv,adjMatrix=Nmatrix))
confint(zut1, parm="(Intercept)_1")
}
}
{
data("sleepstudy",package = "lme4")
(mod0 <- fitme(Reaction ~ Days + AR1(1|Days), data = sleepstudy, fixed=list(ARphi=0.5,lambda=1000)))
# bc the unconstrained fit has essentially lambda=0, which cause later ambiguities
# Don't give to much importance to the result: the fixed Days capture the trend and the 'correlated' rnaefs are not really correlated
# plot(ranef(mod0)[[1]])
sleepmv <- sleepstudy
set.seed(124)
sleepmv$Reaction <- simulate(mod0)
sleepmv$reac2 <- simulate(mod0)
sleepmv$days2 <- sleepmv$Days
cat(crayon::yellow("non-standard REML (both ways); "))
(re1 <- fitme(Reaction ~ 1 + (1|Days), REMLformula=~Days, data = sleepmv, method="REML"))
(re2 <- fitme(reac2 ~ Days + (1|days2), REMLformula=~1, data = sleepmv, method="REML"))
(rezut1 <- fitmv(submodels=list(mod1=list(Reaction ~ 1 + (1|Days), REMLformula=~Days),
mod2=list(reac2 ~ Days + (1|days2), REMLformula=~1)),
data=sleepmv, method="REML"))
testthat::expect_true(diff(range(logLik(rezut1), logLik(re1)+logLik(re2)))<1e-08)
cat(crayon::yellow("AR1; "))
(mod1 <- fitme(Reaction ~ Days + AR1(1|Days), data = sleepmv))
(mod2 <- fitme(reac2 ~ Days + AR1(1|Days), data = sleepmv))
(zut1 <- fitmv(submodels=list(mod1=list(Reaction ~ Days + AR1(1|Days)),
mod2=list(reac2 ~ Days + AR1(1|days2))),
data=sleepmv))
testthat::expect_true(diff(range(logLik(zut1), logLik(mod1)+logLik(mod2)))<1e-08)
(zut1 <- fitmv(submodels=list(mod1=list(Reaction ~ Days + AR1(1|Days)),
mod2=list(reac2 ~ Days + AR1(1|Days))),
data=sleepmv))
(zut2 <- fitmv(submodels=list(mod2=list(reac2 ~ Days + AR1(1|Days)),
mod1=list(Reaction ~ Days + AR1(1|Days))),
data=sleepmv))
testthat::expect_true(diff(range(logLik(zut1), logLik(zut2)))<1e-08)
}
if(FALSE) { cat(crayon::yellow("simulation study; "))
{
set.seed(123)
replic <- function(nind=1000L, cor=-0.5, lambda=c(1,0.2), binsize=1L, method=if (binsize<4L) {"PQL/L"} else {"ML"}, return.fit=FALSE, ...) {
cat(".")
u1 <- rnorm(nind)
u2 <- rnorm(nind)
lam1 <- lambda[1]
lam2 <- lambda[2]
L <- chol(matrix(c(lam1,sqrt(lam1*lam2)*cor,sqrt(lam1*lam2)*cor,lam2),ncol=2))
v <- t(L) %*% rbind(u1,u2)
eta1 <- v[1,]
eta2 <- 1+v[2,]
surv <- rbinom(length(eta1),binsize,inv.logit(eta1)) # was it boot::inv.logit ?? Also called once more below.
feco <- rpois(length(eta2), lambda = exp(eta2))
lfh <- data.frame(id=seq_along(eta1), surv=surv, feco=feco, binsize=binsize)
(fitlfh <- fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id), family=poisson())),
data=lfh, method=method, ...))
if (return.fit) {
return(fitlfh)
} else {
corr <- VarCorr(fitlfh)[2,"Corr."]
cat(corr)
corr
}
}
replicate(3L, replic(nind=1000L))
replicate(10L, replic(nind=3000L, binsize=10)) # much better...
replicate(10L, replic(nind=300L, binsize=10)) # much better...
PQLLdist <- replicate(50L, replic(nind=300L, binsize=2)) # passable
PQLdist <- replicate(50L, replic(nind=300L, binsize=2, method="PQL")) # passable
# one can guess a *small* advantage of REML:
plot(ecdf(PQLdist))
plot(ecdf(PQLLdist), add=TRUE, col="red")
#
PQLdist <- replicate(10L, replic(nind=1000L, binsize=1, method="PQL")) # binary...
PQLdist <- replicate(10L, replic(nind=1000L, cor=0, binsize=1, method="PQL")) # binary...
PQLdist <- replicate(10L, replic(nind=300L, cor=-0.5, binsize=1, method="PQL")) # binary & small sample size -> extreme negative
PQLdist <- replicate(10L, replic(nind=10000L, cor=-0.5, binsize=1, method="PQL")) # large sample size does not improve this
(zut <- replic(nind=300L, cor=-0.5, binsize=1, method="PQL", return.fit=TRUE, verbose=c(TRACE=TRUE)))
replic(nind=1000L, binsize=1, return.fit=TRUE)
plot(eta1,eta2)
cov(t(v))
plot(eta1,surv)
plot(eta2,feco)
plot(eta1,feco)
cov(cbind(surv,feco))
}
{
set.seed(123)
replic2 <- function(nind=1000L, lambda=0.2, binsize=1L, method=if (binsize<4L) {"PQL/L"} else {"ML"}, return.fit=FALSE, ...) {
cat(".")
u1 <- sqrt(lambda)*rnorm(nind)
eta_f <- 1+u1
feco <- rpois(length(eta_f), lambda = exp(eta_f))
surv <- rbinom(length(eta_f),binsize,inv.logit(-(log(feco)-1)))
lfh <- data.frame(id=seq_along(eta_f), id2=seq_along(eta_f), surv=surv, feco=feco, binsize=binsize)
(fitlfh <- fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id), family=poisson())),
data=lfh, method=method, ...))
if (return.fit) {
return(fitlfh)
} else {
corr <- VarCorr(fitlfh)[2,"Corr."]
cat(corr)
corr
}
}
replicate(10L, replic2(nind=1000L))
replicate(10L, replic2(nind=1000L, binsize=10))
replicate(10L, replic2(nind=3000L))
replicate(10L, replic2(nind=3000L, binsize=10)) #
(zut <- replic2(nind=3000L, binsize=10, return.fit=TRUE))
with(zut$data, plot(surv,feco))
fitme(cbind(surv,binsize-surv) ~ 1+(1|id), family=binomial(),
data=zut$data)
fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(1|id), family=binomial()),
list(feco ~ 1+(1|id2), family=poisson())),
data=zut$data)
fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id), family=poisson())),
data=zut$data)
(debu <- fitmv(submodels=list(list(cbind(surv,binsize-surv) ~ 1+(mv(1,2)|id), family=binomial()),
list(feco ~ 1+(mv(1,2)|id2), family=poisson())),
data=zut$data[1:20,]))
plot(eta1,eta2)
cov(t(v))
plot(eta1,surv)
plot(eta2,feco)
plot(eta1,feco)
cov(cbind(surv,feco))
}
{
simfun <- function(nind=570L, cor=-0.5, lambda=c(0.2,0.1), return.fit=FALSE, ...) {
cat(".")
u <- rnorm(2*nind)
lam1 <- lambda[1]
lam2 <- lambda[2]
L <- t(chol(matrix(c(lam1,sqrt(lam1*lam2)*cor,sqrt(lam1*lam2)*cor,lam2),ncol=2)))
v <- L %*% matrix(u,nrow=2)
if (return.fit) plot(t(v))
lfh <- data.frame(id=seq_len(nind), id2=seq_len(nind), feco= rpois(nind, lambda = exp(1+v[1,])),
growth=rgamma(nind,shape=1/0.2, scale=0.2*exp(1+v[2,]))) # mean=exp(1+v[2,]), var= 0.2*mean^2
(fitlfh <- fitmv(submodels=list(list(feco ~ 1+(0+mv(1,2)|id), family=poisson()),
list(growth ~ 1+(0+mv(1,2)|id), family=Gamma(log))),
data=lfh, method=method, ...))
if (return.fit) {
return(fitlfh)
} else {
corr <- VarCorr(fitlfh)[2,"Corr."]
cat(corr)
corr
}
}
set.seed(123)
(zut <- simfun(cor=-0.5,return.fit=TRUE))
plot(zut$data)
replicate(3L, simfun(3000))
replicate(3L, simfun(10000))
}
}
if (FALSE) {
library(aster)
data(echinacea)
asNA.acyclic <- function(data, pred) {
order_pred <- order(pred)
predvars <- names(pred)
ord_desc_vars <- names(pred[order_pred])
for (ordered_it in order_pred) {
descvar <- ord_desc_vars[ordered_it]
predecessor <- pred[descvar]
if (predecessor>0L) {
predvar <- predvars[predecessor]
# Predecessor variable being 0 (no survival, or no flowers) means
# that descendant variable cannot be observed:
data[is.na(data[[predvar]]) | data[[predvar]]==0L, descvar] <- NA
}
}
return(data)
}
# Note the names, important here:
varpred <- c(ld02=0, ld03=1, ld04=2, fl02=1, fl03=2, fl04=3,
hdct02=4, hdct03=5, hdct04=6)
NAechin <- asNA.acyclic(echinacea, pred=varpred)
yearvars <- c("ld02", "ld03", "ld04")
byyear <- reshape(NAechin, varying = list(yearvars), direction = "long",
timevar = "varld",times = as.factor(yearvars), v.names = "respld")
yearvars <- c("fl02", "fl03", "fl04")
flbyyear <- reshape(NAechin, varying = list(yearvars), direction = "long",
timevar = "varfl",times = as.factor(yearvars), v.names = "respfl")
yearvars <- c("hdct02", "hdct03", "hdct04")
hdctbyyear <- reshape(NAechin, varying = list(yearvars), direction = "long",
timevar = "varhdct",times = as.factor(yearvars), v.names = "resphdct")
byyear$varfl <- flbyyear$varfl
byyear$varhdct <- hdctbyyear$varhdct
byyear$respfl <- flbyyear$respfl
byyear$resphdct <- hdctbyyear$resphdct
byyear <- byyear[ ,c(12,4:6,10,13,14,11,15,16)]
head(byyear)
## ----asMM, size="small"---------------------------------------------------------------------------------------------------------------------------------------
asMM <- fitmv(submodels=list(
list(respld ~ varld + nsloc + ewloc+(1|id), family=binomial()),
list(respfl ~ varfl + nsloc + ewloc+(1|id), family=binomial()),
list(resphdct ~ varhdct + nsloc + ewloc+pop, family=Tpoisson())),
data=byyear)
if (requireNamespace("multcomp", quietly = TRUE)) {
library(multcomp)
#summary(glht(asMM,mcp("varld" = "Tukey"), coef.=fixef.HLfit)) # documented limitation
}
}
spaMM.options(spaMM_tol=spaMM_tol)
summary(warnings())
Try the spaMM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.