# NOT RUN {
library(lme4)
#############################################################################
# EXAMPLE 1: Dataset data.bifie01 | TIMSS 2011
#############################################################################
data(data.bifie01)
dat <- data.bifie01
set.seed(987)
# create dataset with replicate weights and plausible values
bdat1 <- BIFIEsurvey::BIFIE.data.jack( data=dat, jktype="JK_TIMSS", jkzone="JKCZONE",
jkrep="JKCREP", wgt="TOTWGT", pv_vars=c("ASMMAT","ASSSCI") )
# create dataset without plausible values and ignoring weights
bdat2 <- BIFIEsurvey::BIFIE.data.jack( data=dat, jktype="JK_RANDOM", ngr=10 )
#=> standard errors from ML estimation
#***********************************************
# Model 1: Random intercept model
#--- Model 1a: without weights, first plausible value
mod1a <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat2, dep="ASMMAT01",
formula.fixed=~ 1, formula.random=~ 1, idcluster="idschool",
wgtlevel2="one", se=FALSE )
summary(mod1a)
#--- Model 1b: estimation in lme4
mod1b <- lme4::lmer( ASMMAT01 ~ 1 + ( 1 | idschool), data=dat, REML=FALSE)
summary(mod1b)
#--- Model 1c: Like Model 1a but for five plausible values and ML inference
mod1c <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat1, dep="ASMMAT",
formula.fixed=~ 1, formula.random=~ 1, idcluster="idschool",
wgtlevel2="one", se=FALSE )
summary(mod1c)
#--- Model 1d: weights and sampling design and all plausible values
mod1d <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat1, dep="ASMMAT",
formula.fixed=~ 1, formula.random=~ 1, idcluster="idschool",
wgtlevel2="SCHWGT" )
summary(mod1d)
#***********************************************
# Model 2: Random slope model
#--- Model 2a: without weights
mod2a <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat2, dep="ASMMAT01",
formula.fixed=~ female + ASBG06A, formula.random=~ ASBG06A,
idcluster="idschool", wgtlevel2="one", se=FALSE )
summary(mod2a)
#--- Model 2b: estimation in lme4
mod2b <- lme4::lmer( ASMMAT01 ~ female + ASBG06A + ( 1 + ASBG06A | idschool),
data=dat, REML=FALSE)
summary(mod2b)
#--- Model 2c: weights and sampling design and all plausible values
mod2c <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat1, dep="ASMMAT",
formula.fixed=~ female + ASBG06A, formula.random=~ ASBG06A,
idcluster="idschool", wgtlevel2="SCHWGT", maxiter=500, se=FALSE)
summary(mod2c)
#--- Model 2d: Uncorrelated intecepts and slopes
# constraint for zero covariance between intercept and slope
recov_constraint <- matrix( c(1,2,0), ncol=3 )
mod2d <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat2, dep="ASMMAT01",
formula.fixed=~ female + ASBG06A, formula.random=~ ASBG06A,
idcluster="idschool", wgtlevel2="one", se=FALSE,
recov_constraint=recov_constraint )
summary(mod2d)
#--- Model 2e: Fixed entries in the random effects covariance matrix
# two constraints for random effects covariance
# Cov(Int, Slo)=0 # zero slope for intercept and slope
# Var(Slo)=10 # slope variance of 10
recov_constraint <- matrix( c(1,2,0,
2,2,10), ncol=3, byrow=TRUE)
mod2e <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat2, dep="ASMMAT01",
formula.fixed=~ female + ASBG06A, formula.random=~ ASBG06A,
idcluster="idschool", wgtlevel2="one", se=FALSE,
recov_constraint=recov_constraint )
summary(mod2e)
#############################################################################
# SIMULATED EXAMPLE 2: Two-level regression with random slopes
#############################################################################
#--- (1) simulate data
set.seed(9876)
NC <- 100 # number of clusters
Nj <- 20 # number of persons per cluster
iccx <- .4 # intra-class correlation predictor
theta <- c( 0.7, .3 ) # fixed effects
Tmat <- diag( c(.3, .1 ) ) # variances of random intercept and slope
sig2 <- .60 # residual variance
N <- NC*Nj
idcluster <- rep( 1:NC, each=Nj )
dat1 <- data.frame("idcluster"=idcluster )
dat1$X <- rep( stats::rnorm( NC, sd=sqrt(iccx) ), each=Nj ) +
stats::rnorm( N, sd=sqrt( 1 - iccx) )
dat1$Y <- theta[1] + rep( stats::rnorm(NC, sd=sqrt(Tmat[1,1] ) ), each=Nj ) +
theta[2] + rep( stats::rnorm(NC, sd=sqrt(Tmat[2,2])), each=Nj )) * dat1$X +
stats::rnorm(N, sd=sqrt(sig2) )
#--- (2) create design object
bdat1 <- BIFIEsurvey::BIFIE.data.jack( data=dat1, jktype="JK_GROUP", jkzone="idcluster")
summary(bdat1)
#*** Model 1: Random slope model (ML standard errors)
#- estimation using BIFIE.twolevelreg
mod1a <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat1, dep="Y",
formula.fixed=~ 1+X, formula.random=~ 1+X, idcluster="idcluster",
wgtlevel2="one", se=FALSE )
summary(mod1a)
#- estimation in lme4
mod1b <- lme4::lmer( Y ~ X + ( 1+X | idcluster), data=dat1, REML=FALSE )
summary(mod1b)
#- using Jackknife for inference
mod1c <- BIFIEsurvey::BIFIE.twolevelreg( BIFIEobj=bdat1, dep="Y",
formula.fixed=~ 1+X, formula.random=~ 1+X, idcluster="idcluster",
wgtlevel2="one", se=TRUE )
summary(mod1c)
# extract coefficients
coef(mod1a)
coef(mod1c)
# covariance matrix
vcov(mod1a)
vcov(mod1c)
# }
Run the code above in your browser using DataLab