#############################################################################
# EXAMPLE 1: Imputed TIMSS dataset
#############################################################################
data(data.timss1)
data(data.timssrep)
# create BIFIE.dat object
bifieobj <- BIFIEsurvey::BIFIE.data( data.list=data.timss1, wgt=data.timss1[[1]]$TOTWGT,
wgtrep=data.timssrep[, -1 ] )
#****************************
#*** Model 1: Weighted means (as a toy example)
userfct <- function(X,w){
pars <- c( stats::weighted.mean( X[,1], w ),
stats::weighted.mean(X[,2], w ) )
return(pars)
}
res1 <- BIFIEsurvey::BIFIE.by( bifieobj, vars=c("ASMMAT", "migrant", "books"),
userfct=userfct, userparnames=c("MW_MAT", "MW_Migr"),
group="female" )
summary(res1)
# evaluate function in pure R implementation using the use_Rcpp argument
res1b <- BIFIEsurvey::BIFIE.by( bifieobj, vars=c("ASMMAT", "migrant", "books" ),
userfct=userfct, userparnames=c("MW_MAT", "MW_Migr"),
group="female", use_Rcpp=FALSE )
summary(res1b)
#--- statistical inference for a derived parameter (see ?BIFIE.derivedParameters)
# define gender difference for mathematics score (divided by 100)
derived.parameters <- list(
"gender_diff"=~ 0 + I( ( MW_MAT_female1 - MW_MAT_female0 ) / 100 )
)
# inference derived parameter
res1d <- BIFIEsurvey::BIFIE.derivedParameters( res1,
derived.parameters=derived.parameters )
summary(res1d)
if (FALSE) {
#****************************
#**** Model 2: Robust linear model
# (1) start from scratch to formulate the user function for X and w
dat1 <- bifieobj$dat1
vars <- c("ASMMAT", "migrant", "books" )
X <- dat1[,vars]
w <- bifieobj$wgt
library(MASS)
# ASMMAT ~ migrant + books
mod <- MASS::rlm( X[,1] ~ as.matrix( X[, -1 ] ), weights=w )
coef(mod)
# (2) define a user function "my_rlm"
my_rlm <- function(X,w){
mod <- MASS::rlm( X[,1] ~ as.matrix( X[, -1 ] ), weights=w )
return( coef(mod) )
}
# (3) estimate model
res2 <- BIFIEsurvey::BIFIE.by( bifieobj, vars, userfct=my_rlm,
group="female", group_values=0:1)
summary(res2)
# estimate model without computing standard errors
res2a <- BIFIEsurvey::BIFIE.by( bifieobj, vars, userfct=my_rlm,
group="female", se=FALSE)
summary(res2a)
# define a user function with formula language
my_rlm2 <- function(X,w){
colnames(X) <- vars
X <- as.data.frame(X)
mod <- MASS::rlm( ASMMAT ~ migrant + books, weights=w, data=X)
return( coef(mod) )
}
# estimate model
res2b <- BIFIEsurvey::BIFIE.by( bifieobj, vars, userfct=my_rlm2,
group="female", group_values=0:1)
summary(res2b)
#****************************
#**** Model 3: Number of unique values for variables in BIFIEdata
#*** define variables for which the number of unique values should be calculated
vars <- c( "female", "books","ASMMAT" )
#*** define a user function extracting these unqiue values
userfct <- function(X,w){
pars <- apply( X, 2, FUN=function(vv){
length( unique(vv)) } )
# Note that weights are (of course) ignored in this function
return(pars)
}
#*** extract number of unique values
res3 <- BIFIEsurvey::BIFIE.by( bifieobj, vars=vars, userfct=userfct,
userparnames=paste0( vars, "_Nunique"), se=FALSE )
summary(res3)
## Statistical Inference for User Definition Function
## parm Ncases Nweight est
## 1 female_Nunique 4668 78332.99 2.0
## 2 books_Nunique 4668 78332.99 5.0
## 3 ASMMAT_Nunique 4668 78332.99 4613.4
# number of unique values in each of the five imputed datasets
res3$output$parsrepM
## [,1] [,2] [,3] [,4] [,5]
## [1,] 2 2 2 2 2
## [2,] 5 5 5 5 5
## [3,] 4617 4619 4614 4609 4608
#****************************
#**** Model 4: Estimation of a lavaan model with BIFIE.by
#* estimate model in lavaan
data0 <- data.timss1[[1]]
# define lavaan model
lavmodel <- "
ASSSCI ~ likesc
ASSSCI ~~ ASSSCI
likesc ~ female
likesc ~~ likesc
female ~~ female
"
mod0 <- lavaan::lavaan(lavmodel, data=data0, sampling.weights="TOTWGT")
summary(mod0, stand=TRUE, fit.measures=TRUE)
#* construct input for BIFIE.by
vars <- c("ASSSCI","likesc","female","TOTWGT")
X <- data0[,vars]
mod0 <- lavaan::lavaan(lavmodel, data=X, sampling.weights="TOTWGT")
w <- data0$TOTWGT
#* define user function
userfct <- function(X,w){
X1 <- as.data.frame(X)
colnames(X1) <- vars
X1$studwgt <- w
mod0 <- lavaan::lavaan(lavmodel, data=X1, sampling.weights="TOTWGT")
pars <- coef(mod0)
# extract some fit statistics
pars2 <- lavaan::fitMeasures(mod0)
pars <- c(pars, pars2[c("cfi","tli")])
return(pars)
}
#* test function
res0 <- userfct(X,w)
userparnames <- names(res0)
#* estimate lavaan model with replicated sampling weights
res1 <- BIFIEsurvey::BIFIE.by( bifieobj, vars=vars, userfct=userfct,
userparnames=userparnames, use_Rcpp=FALSE )
summary(res1)
}
Run the code above in your browser using DataLab