3.3-11

alexanderrobitzsch · Feb 28, 2019 · 555bc2d · 555bc2d
1 parent 30f928b
commit 555bc2d
Show file tree

Hide file tree

Showing 21 changed files with 284 additions and 100 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: sirt
 Type: Package
 Title: Supplementary Item Response Theory Models
-Version: 3.3-6
-Date: 2019-02-27 22:42:09
+Version: 3.3-11
+Date: 2019-02-28 17:52:57
 Author: Alexander Robitzsch [aut,cre] (<https://orcid.org/0000-0002-8226-3132>)
 Maintainer: Alexander Robitzsch <[email protected]>
 Description: 
@@ -25,11 +25,11 @@ Description:
 Depends:
     R (>= 3.1)
 Imports:
-    CDM (>= 6.0), coda, graphics, lavaan, lavaan.survey, 
-    MASS, methods, mirt, mvtnorm, Rcpp, sfsmisc, stats, 
-    survey, TAM (>= 2.5-5), utils
+    CDM (>= 6.0), coda, graphics, lavaan, MASS, methods, mirt, 
+    mvtnorm, Rcpp, stats, TAM (>= 2.5-5), utils
 Suggests:
-    igraph, Matrix, miceadds, pbivnorm, psych, sm
+    igraph, Matrix, lavaan.survey, miceadds, pbivnorm, psych, 
+    sfsmisc, sm, survey
 LinkingTo: 
     Rcpp, RcppArmadillo
 URL: 

diff --git a/NAMESPACE b/NAMESPACE
@@ -143,7 +143,6 @@ importFrom(lavaan, parameterEstimates)
 importFrom(lavaan, parameterTable)
 importFrom(lavaan, standardizedSolution)
 importFrom(lavaan, sem)
-importFrom(lavaan.survey, lavaan.survey)
 importFrom(MASS, ginv)
 importFrom(mirt, extract.item)
 importFrom(mirt, mirt)
@@ -153,8 +152,8 @@ importFrom(mirt, probtrace)
 importFrom(mvtnorm, dmvnorm)
 importFrom(mvtnorm, pmvnorm)
 # importFrom(mvtnorm, rmvnorm) # replaced by CDM::CDM_rmvnorm
-importFrom(sfsmisc, QUnif)
-importFrom(survey, svydesign)
+
+
 
  #*** not imported anymore
  # import(qgraph)
@@ -169,14 +168,16 @@ importFrom(survey, svydesign)
 # require_namespace_msg("another_package")
 
 # importFrom(igraph, graph.edgelist)
+# importFrom(lavaan.survey, lavaan.survey)
 # importFrom(Matrix, bdiag)
 # importFrom(pbivnorm, pbivnorm)
 # importFrom(psych, cor.smooth)
 # importFrom(psych, fa)
 # importFrom(psych, omega)
+# importFrom(sfsmisc, QUnif)
 # importFrom(sm, h.select)
 # importFrom(sm, sm.binomial)
-
+# importFrom(survey, svydesign)
 
 
 ####################################################
@@ -241,6 +242,7 @@ export(lsdm)
 export(lsem.estimate)
 export(lsem.MGM.stepfunctions)
 export(lsem.permutationTest)
+export(lsem_local_weights)
 export(marginal.truescore.reliability)
 export(mcmc.2pno)
 export(mcmc.2pno.ml)

diff --git a/R/RcppExports.R b/R/RcppExports.R
@@ -1,5 +1,5 @@
 ## File Name: RcppExports.R
-## File Version: 3.003006
+## File Version: 3.003011
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 

diff --git a/R/invariance_alignment_cfa_config.R b/R/invariance_alignment_cfa_config.R
@@ -1,5 +1,5 @@
 ## File Name: invariance_alignment_cfa_config.R
-## File Version: 0.16
+## File Version: 0.17
 
 
 invariance_alignment_cfa_config <- function(dat, group, weights=NULL, ...)
@@ -25,7 +25,7 @@ invariance_alignment_cfa_config <- function(dat, group, weights=NULL, ...)
             weights_gg <- weights[ind_gg]
         }
         res <- invariance_alignment_cfa_config_estimate(dat_gg=dat_gg,
-                        weights=weights_gg, ...)
+                        weights_gg=weights_gg, ...)
         nu[gg, ind_gg] <- res$nu
         lambda[gg, ind_gg] <- res$lambda
         err_var[gg, ind_gg] <- res$err_var

diff --git a/R/lsem.estimate.R b/R/lsem.estimate.R
@@ -1,14 +1,13 @@
 ## File Name: lsem.estimate.R
-## File Version: 0.894
+## File Version: 0.917
 
 # estimate LSEM model
 lsem.estimate <- function( data, moderator, moderator.grid,
         lavmodel, type="LSEM", h=1.1,
         residualize=TRUE, fit_measures=c("rmsea","cfi","tli","gfi","srmr"),
-        standardized=FALSE,
-        standardized_type="std.all",
+        standardized=FALSE, standardized_type="std.all",
         lavaan_fct="sem", sufficient_statistics=FALSE,
-        eps=1E-8, verbose=TRUE, ... )
+        use_lavaan_survey=FALSE, pseudo_weights=0, eps=1E-8, verbose=TRUE, ... )
 {
 
     CALL <- match.call()
@@ -20,6 +19,9 @@ lsem.estimate <- function( data, moderator, moderator.grid,
             stop("standardized=TRUE cannot be applied for type='MGM'")
         }
     }
+    if (sufficient_statistics){ use_lavaan_survey <- FALSE }
+    if (pseudo_weights>0){ use_lavaan_survey <- FALSE }
+
     # group moderator if type="MGM"
     out <- lsem_group_moderator( data=data, type=type, moderator.grid=moderator.grid,
                 moderator=moderator, residualize=residualize, h=h )
@@ -39,12 +41,11 @@ lsem.estimate <- function( data, moderator, moderator.grid,
     # unweighted fit of lavaan model
     dat <- data
     lavmodel__ <- lavmodel
-    if (lavaan_fct=="sem"){
-        lavfit <- lavaan::sem(model=lavmodel__, data=dat,  ... )
-    }
-    if (lavaan_fct=="lavaan"){
-        lavfit <- lavaan::lavaan(model=lavmodel__, data=dat,  ... )
-    }
+
+    #* fit lavaan model
+    lavaan_est_fun <- lsem_define_lavaan_est_fun(lavaan_fct=lavaan_fct)
+    lavfit <- lavaan_est_fun(model=lavmodel__, data=dat,  ... )
+
     # extract variables which are in model and data frame
     partable <- pars <- lavaan::parameterEstimates(lavfit)
     variables_model <- intersect( union( partable$lhs, partable$rhs ), colnames(dat) )
@@ -65,7 +66,9 @@ lsem.estimate <- function( data, moderator, moderator.grid,
                     moderator.grid=moderator.grid, verbose=verbose, pars=pars,
                     standardized=standardized, variables_model=variables_model,
                     sufficient_statistics=sufficient_statistics,
-                    lavaan_fct=lavaan_fct, lavmodel=lavmodel, ... )
+                    lavaan_fct=lavaan_fct, lavmodel=lavmodel,
+                    use_lavaan_survey=use_lavaan_survey, pseudo_weights=pseudo_weights,
+                    ... )
     parameters <- out2$parameters
     rownames(parameters) <- paste0( parameters$par, "__", parameters$grid_index )
 
@@ -87,8 +90,7 @@ lsem.estimate <- function( data, moderator, moderator.grid,
     dfr0 <- data.frame("M"=mean( data[,moderator], na.rm=TRUE ),
                 "SD"=out$sd.moderator,
                 "min"=min( data[, moderator ], na.rm=TRUE ),
-                "max"=max( data[, moderator ], na.rm=TRUE )
-                            )
+                "max"=max( data[, moderator ], na.rm=TRUE ) )
     obji <- rbind( dfr0, dfr )
     rownames(obji) <- NULL
     moderator.stat <- data.frame("variable"=c("moderator",
@@ -111,8 +113,8 @@ lsem.estimate <- function( data, moderator, moderator.grid,
                     fit_measures=fit_measures, s1=s1, s2=s2,
                     standardized=standardized,
                     standardized_type=standardized_type,
-                    lavaan_fct=lavaan_fct,
-                    type=type, CALL=CALL )
+                    lavaan_fct=lavaan_fct, use_lavaan_survey=use_lavaan_survey,
+                    pseudo_weights=pseudo_weights, type=type, CALL=CALL )
     class(res) <- "lsem"
     return(res)
 }
diff --git a/R/lsem.permutationTest.R b/R/lsem.permutationTest.R
@@ -1,5 +1,5 @@
 ## File Name: lsem.permutationTest.R
-## File Version: 0.38
+## File Version: 0.44
 
 ############################################
 # permutation test for LSEM model
@@ -31,6 +31,8 @@ lsem.permutationTest <- function( lsem.object, B=1000, residualize=TRUE,
     arglist$standardized <- lsem.object$standardized
     arglist$standardized_type <- lsem.object$standardized_type
     arglist$lavaan_fct <- lsem.object$lavaan_fct
+    use_lavaan_survey <- arglist$use_lavaan_survey <- lsem.object$use_lavaan_survey
+    arglist$pseudo_weights <- lsem.object$pseudo_weights
 
     data1 <- data0 <- object$data
     parameters <- object$parameters
@@ -128,6 +130,7 @@ lsem.permutationTest <- function( lsem.object, B=1000, residualize=TRUE,
                      moderator.grid=object$moderator.grid,
                      h=object$h, bw=object$bw, N=object$N,
                      nonconverged_rate=nonconverged_rate,
+                     use_lavaan_survey=use_lavaan_survey,
                      B=B, s1=s1, s2=s2, lavmodel=object$lavmodel, CALL=CALL
                             )
     class(res) <- "lsem.permutationTest"

diff --git a/R/lsem_define_lavaan_est_fun.R b/R/lsem_define_lavaan_est_fun.R
@@ -0,0 +1,16 @@
+## File Name: lsem_define_lavaan_est_fun.R
+## File Version: 0.03
+
+lsem_define_lavaan_est_fun <- function(lavaan_fct)
+{
+    if (lavaan_fct=="sem"){
+        lavaan_est_fun <- lavaan::sem
+    }
+    if (lavaan_fct=="lavaan"){
+        lavaan_est_fun <- lavaan::lavaan
+    }
+    if (lavaan_fct=="cfa"){
+        lavaan_est_fun <- lavaan::cfa
+    }
+    return(lavaan_est_fun)
+}
diff --git a/R/lsem_fitsem.R b/R/lsem_fitsem.R
@@ -1,18 +1,18 @@
 ## File Name: lsem_fitsem.R
-## File Version: 0.44
+## File Version: 0.472
 
-##############################################################
 lsem_fitsem <- function( dat, weights, lavfit,
             fit_measures, NF, G, moderator.grid, verbose,
             pars, standardized, variables_model,
             sufficient_statistics, lavaan_fct, lavmodel,
-            ... )
+            use_lavaan_survey=TRUE, pseudo_weights=0, ... )
 {
 
     parameters <- NULL
     fits <- NULL
+    survey.fit <- lavfit
     pars0 <- pars
-    env_temp <- environment()
+    lavaan_est_fun <- lsem_define_lavaan_est_fun(lavaan_fct=lavaan_fct)
 
     if (verbose){
         cat( "** Fit lavaan model\n")
@@ -25,41 +25,42 @@ lsem_fitsem <- function( dat, weights, lavfit,
     }
 
     for (gg in 1:G){
-        # gg <- 1
         dat$weight <- weights[,gg]
         #***** fit the model using weighted data
         if (! sufficient_statistics){
-            datsvy <- survey::svydesign(id=~index, weights=~weight, data=dat)
-            # assign(x="lavmodel__", value=lavmodel, pos=1)
-            assign_args <- list( x="lavmodel__", value=lavmodel, pos=1)
-            res0 <- do.call( what="assign", args=assign_args)
-            survey.fit <- lavaan.survey::lavaan.survey(lavaan.fit=lavfit,
-                                survey.design=datsvy )
+            if (use_lavaan_survey){
+                survey.fit <- lsem_fitsem_raw_data_lavaan_survey(dat=dat,
+                                    lavmodel=lavmodel, lavfit=lavfit)
+            }
+            if (! use_lavaan_survey){
+                res <- lsem_fitsem_raw_data_define_pseudo_weights(dat=dat,
+                            pseudo_weights=pseudo_weights)
+                dat1 <- res$dat
+                sampling_weights <- res$sampling_weights
+                # use starting values
+                partable <- lavaan::parameterTable(object=survey.fit)
+                partable$start <- partable$est
+                survey.fit <- lavaan_est_fun(model=partable, data=dat1,
+                                sampling.weights=sampling_weights, ... )
+            }
         }
         #***** fit the model using sufficient statistics
         if (sufficient_statistics){
             res <- lsem_weighted_mean( x=dat[, variables_model], weights=dat$weight )
             wmean <- res$mean
             res <- lsem_weighted_cov( x=dat[, variables_model], weights=dat$weight )
             wcov <- res$cov
-            Nobs <- round( res$Nobs )
-            if (lavaan_fct=="sem"){
-                survey.fit <- lavaan::sem(model=lavmodel, sample.cov=wcov,
+            Nobs <- round(res$Nobs)
+            survey.fit <- lavaan_est_fun(model=lavmodel, sample.cov=wcov,
                                 sample.mean=wmean, sample.nobs=Nobs, ... )
-            }
-            if (lavaan_fct=="lavaan"){
-                survey.fit <- lavaan::lavaan(model=lavmodel, sample.cov=wcov,
-                                sample.mean=wmean, sample.nobs=Nobs, ... )
-            }
         }
 
-        dfr.gg <- pars <- lavaan::parameterEstimates(survey.fit)
+        dfr.gg <- pars <- lavaan::parameterEstimates(object=survey.fit)
         if (standardized){
-            sol <- lavaan::standardizedSolution( survey.fit )
+            sol <- lavaan::standardizedSolution(object=survey.fit)
             colnames(sol)[ which( colnames(sol)=="est.std" ) ] <- "est"
             sol$lhs <- paste0( "std__", sol$lhs)
             pars <- sirt_rbind_fill( x=pars, y=sol )
-            # pars <- plyr::rbind.fill( pars, sol )
             dfr.gg <- pars
         }
         pars <- paste0( pars$lhs, pars$op, pars$rhs )
@@ -68,10 +69,10 @@ lsem_fitsem <- function( dat, weights, lavfit,
         dfr.gg <- dfr.gg[ ind, ]
         dfr.gg <- data.frame("grid_index"=gg, "moderator"=moderator.grid[gg],
                           "par"=pars0, "parindex"=1:NP, dfr.gg    )
+        est_fit <- lavaan::fitMeasures(object=survey.fit, fit.measures=fit_measures )
         dfr.gg0 <- data.frame("grid_index"=gg, "moderator"=moderator.grid[gg],
                           "par"=fit_measures, "parindex"=NP + 1:NF,
-                          "est"=lavaan::fitMeasures(survey.fit, fit.measures=fit_measures ),
-                          "op"="fit" )
+                          "est"=est_fit, "op"="fit" )
         vars <- setdiff( colnames(dfr.gg), colnames(dfr.gg0) )
         for (vv in vars){ dfr.gg0[,vv] <- NA }
         dfr.gg <- rbind( dfr.gg, dfr.gg0[, colnames(dfr.gg) ] )
@@ -93,6 +94,5 @@ lsem_fitsem <- function( dat, weights, lavfit,
     res <- list( parameters=parameters )
     return(res)
 }
-#######################################################################
 
 lsem.fitsem <- lsem_fitsem
diff --git a/R/lsem_fitsem_raw_data_define_pseudo_weights.R b/R/lsem_fitsem_raw_data_define_pseudo_weights.R
@@ -0,0 +1,18 @@
+## File Name: lsem_fitsem_raw_data_define_pseudo_weights.R
+## File Version: 0.01
+
+lsem_fitsem_raw_data_define_pseudo_weights <- function(dat, pseudo_weights)
+{
+    sampling_weights <- "weight"
+    if (pseudo_weights>0){
+        weights <- dat$weight
+        weights <- round( weights * pseudo_weights )
+        N <- nrow(dat)
+        ind <- rep(1:N, weights)
+        dat <- dat[ind,]
+        sampling_weights <- NULL
+    }
+    #-- output
+    res <- list(dat=dat, sampling_weights=sampling_weights)
+    return(res)
+}
diff --git a/R/lsem_fitsem_raw_data_lavaan_survey.R b/R/lsem_fitsem_raw_data_lavaan_survey.R
@@ -0,0 +1,15 @@
+## File Name: lsem_fitsem_raw_data_lavaan_survey.R
+## File Version: 0.03
+
+lsem_fitsem_raw_data_lavaan_survey <- function(dat, lavmodel, lavfit)
+{
+    TAM::require_namespace_msg("lavaan.survey")
+    TAM::require_namespace_msg("survey")
+    datsvy <- survey::svydesign(id=~index, weights=~weight, data=dat)
+    # assign(x="lavmodel__", value=lavmodel, pos=1)
+    assign_args <- list( x="lavmodel__", value=lavmodel, pos=1)
+    res0 <- do.call( what="assign", args=assign_args)
+    survey.fit <- lavaan.survey::lavaan.survey(lavaan.fit=lavfit,
+                                survey.design=datsvy )
+    return(survey.fit)
+}
diff --git a/R/lsem_local_weights.R b/R/lsem_local_weights.R
@@ -0,0 +1,36 @@
+## File Name: lsem_local_weights.R
+## File Version: 0.05
+
+lsem_local_weights <- function(data.mod, moderator.grid, h)
+{
+    eps <- 1E-8
+    N <- length(data.mod)
+    # select nearest neighbor in moderator group for calculating residuals
+    G <- length(moderator.grid)
+    modgrid_index <- rep(1,N)
+    for (gg in 2:G){
+        modgrid_index <- ifelse( abs( data.mod - moderator.grid[ modgrid_index ] ) <
+                    abs( data.mod - moderator.grid[ gg ] ),
+                    modgrid_index, gg )
+    }
+    # compute weights for every grid point gg
+    weights <- matrix( NA, nrow=N, ncol=G )
+    sd.moderator <- stats::sd( data.mod, na.rm=TRUE)
+    bw <- h * sd.moderator * N^(-1/5)
+    moderator.density <- stats::density( data.mod, from=min(moderator.grid),
+                            to=max(moderator.grid ), n=G )$y
+    moderator.density <- data.frame( moderator=moderator.grid,
+                            wgt=moderator.density / sum(moderator.density) )
+
+    for (gg in 1:G){
+        xgg <- moderator.grid[gg]
+        wgt <- stats::dnorm( data.mod, mean=xgg, sd=bw ) /
+                    stats::dnorm( xgg, mean=xgg, sd=bw )
+        weights[,gg] <- ifelse( wgt < eps, eps, wgt )
+    }
+    #--- output
+    res <- list(weights=weights, N=N, G=G, modgrid_index=modgrid_index,
+                sd.moderator=sd.moderator, bw=bw, moderator.density=moderator.density)
+    return(res)
+}
+