BIFIEsurvey (version 1.5-0)

BIFIE.twolevelreg: Two Level Regression

Description

This function computes the hierarchical two level model with random intercepts and random slopes. The full maximum likelihood estimation is conducted by means of an EM algorithm (Raudenbush & Bryk, 2002).

Usage

BIFIE.twolevelreg( BIFIEobj , dep , formula.fixed , formula.random , idcluster , 
   wgtlevel2 = NULL , wgtlevel1 = NULL , group=NULL , group_values=NULL , se=TRUE ,
   globconv = 1E-6 , maxiter = 1000 )    

## S3 method for class 'BIFIE.twolevelreg':
summary(object,digits=4,...)

## S3 method for class 'BIFIE.twolevelreg':
coef(object,...)

## S3 method for class 'BIFIE.twolevelreg':
vcov(object,...)

Arguments

BIFIEobj

Object of class BIFIEdata

dep

String for the dependent variable in the regression model

formula.fixed

An Rformula for fixed effects

formula.random

An Rformula for random effects

idcluster

Cluster identifier

wgtlevel2

Name of Level 2 weight variable

wgtlevel1

Name of Level 1 weight variable. This is optional. If it is not provided, wgtlevel is calculated from the total weight and wgtlevel2.

group

Optional grouping variable

group_values

Optional vector of grouping values. This can be omitted and grouping values will be determined automatically.

se

Optional logical indicating whether statistical inference based on replication should be employed.

globconv

Convergence criterion for maximum parameter change

maxiter

Maximum number of iterations

object

Object of class BIFIE.twolevelreg

digits

Number of digits for rounding output

...

Further arguments to be passed

Value

• A list with following entries
• stat
  Data frame with coefficients and different sources of variance.
• output
  Extensive output with all replicated statistics
• ...
  More values

Details

The implemented random slope model can be written as $$y_{ij} = \bold{X}_{ij} \bold{\gamma} + \bold{Z}_{ij} \bold{u}_j + \varepsilon_{ij}$$ where $y_{ij}$ is the dependent variable, $\bold{X}_{ij}$ includes the fixed effects predictors (specified by formula.fixed) and $\bold{Z}_{ij}$ includes the random effects predictors (specified by formula.random). The random effects $\bold{u}_j$ follow a multivariate normal distribution. The function also computes a variance decomposition of explained variance due to fixed and random effects for the within and the between level. This variance decomposition is conducted for the predictor matrices $\bold{X}$ and $\bold{Z}$. It is assumed that $\bold{X}_{ij} = \bold{X}_j^B + \bold{X}_{ij}^W$. The different sources of variance are computed by formulas as proposed in Snijders and Bosker (2012, Ch. 7).

References

Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical linear models: Applications and data analysis methods. Thousand Oaks: Sage. Snijders, T. A. B., & Bosker, R. J. (2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling. Thousand Oaks: Sage.

See Also

lme4::lmer

Examples

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 <- 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 <- 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 <- 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: weights and sampling design and all plausible values
mod1c <- BIFIE.twolevelreg( BIFIEobj=bdat1 , dep="ASMMAT" , 
                formula.fixed = ~ 1 , formula.random = ~ 1 , idcluster="idschool" , 
                wgtlevel2 = "SCHWGT" )
summary(mod1c)

#***********************************************
# Model 2: Random slope model

#--- Model 1a: without weights
mod2a <- 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 1c: weights and sampling design and all plausible values
mod2c <- BIFIE.twolevelreg( BIFIEobj=bdat1 , dep="ASMMAT" , 
                formula.fixed = ~  female +  ASBG06A , formula.random = ~ ASBG06A , 
                idcluster="idschool" , wgtlevel2 = "SCHWGT" , maxiter=500 , se=FALSE)
summary(mod2c)

#############################################################################
# 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( rnorm( NC , sd = sqrt(iccx) ) , each=Nj ) + rnorm( N , sd = sqrt( 1 - iccx) )
dat1$Y <- theta[1] + rep( rnorm(NC , sd = sqrt(Tmat[1,1] ) ) , each=Nj ) + 
        ( theta[2] + rep( rnorm(NC , sd = sqrt(Tmat[2,2] ) ) , each=Nj ) ) * dat1$X +
                rnorm(N , sd= sqrt(sig2) )

#--- (2) create design object 
bdat1 <- 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 <- 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 <- 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)
# variance matrix
vcov(mod1a)
vcov(mod1c)