frm.Rd
The factored regression model (FRM) allows the estimation of the linear regression model (with normally distributed residuals) and the generalized logistic regression model (logistic regression for dichotomous outcomes). Missing values in covariates are handled by posing a conditional univariate distribution for each covariate. The approach follows Ibrahim (1990), Ibrahim, Chen, Lipsitz and Herring (2005), Lee and Mitra (2016), and Zhang and Wang (2017) and is applied in Luedtke, Robitzsch, and West (2020a, 2020b). Latent variables and covariates with measurement error or multiple indicators can also be handled within this framework (see Examples 3, 4 and 5).
Missing values are handled by numerical integration in frm_em
(see also Allison, 2012). The user has to specify an integration grid for each
variable (defined in argument nodes
for each model).
Standard error estimates in frm_em
are obtained by a numerical differentiation of
the Fisher score function (see Jamshidian & Jennrich, 2000).
The function frm_fb
employs a fully Bayesian approach with noninformative
prior distribution. This function imputes missing values in the models from
the posterior distributions. Imputed datasets can be extracted by the
function frm2datlist
.
The current functionality only support missing values on continuous covariates (accommodating skewness and only positive values), dichotomous covariates and ordinal covariates.
Multilevel models (using model="mlreg"
) for normally distributed
(outcome="normal"
) and ordinal variables (outcome="probit"
)
as well as variables at higher levels (using argument variable_level
)
are accommodated.
The handling of nominal covariates will be included in future mdmb package versions.
# Factored regression model: Numerical integration with EM algorithm frm_em(dat, dep, ind, weights=NULL, verbose=TRUE, maxiter=500, conv_dev=1e-08, conv_parm=1e-05, nodes_control=c(11,5), h=1e-04, use_grad=2, update_model=NULL) # S3 method for frm_em coef(object, ...) # S3 method for frm_em logLik(object, ...) # S3 method for frm_em summary(object, digits=4, file=NULL, ...) # S3 method for frm_em vcov(object, ...) # Factored regression model: Fully Bayesian estimation frm_fb(dat, dep, ind, weights=NULL, verbose=TRUE, data_init=NULL, iter=500, burnin=100, Nimp=10, Nsave=3000, refresh=25, acc_bounds=c(.45,.50), print_iter=10, use_gibbs=TRUE, aggregation=TRUE) # S3 method for frm_fb coef(object, ...) # S3 method for frm_fb plot(x, idparm=NULL, ask=TRUE, ... ) # S3 method for frm_fb summary(object, digits=4, file=NULL, ...) # S3 method for frm_fb vcov(object, ...) frm2datlist(object, as_mids=FALSE) # create list of imputed datasets
dat | Data frame |
---|---|
dep | List containing model specification for dependent variable. The list has arguments (see Examples)
|
ind | List containing a list of univariate conditional models for covariates.
See |
weights | Optional vector of sampling weights |
verbose | Logical indicating whether convergence progress should be displayed. |
maxiter | Maximum number of iterations |
conv_dev | Convergence criterion for deviance |
conv_parm | Convergence criterion for regression coefficients |
nodes_control | Control arguments if nodes are not provided by the user. The first value denote the number of nodes, while the second value denotes the spread of the node distribution defined as the factor of the standard deviation of the observed data. |
h | Step width for numerical differentiation for calculating the covariance matrix |
use_grad | Computation method for gradient in |
update_model | Optional proviously fitted model can be used as an input |
data_init | Initial values for dataset |
iter | Number of iterations |
burnin | Number of burnin iterations |
Nimp | Number of imputed datasets |
Nsave | (Maximum) Number of values to be saved for MCMC chain |
refresh | Number of imputations after which proposal distribution should be updated in Metropolis-Hastings step |
acc_bounds | Bounds for acceptance rates for parameters |
print_iter | Number of imputation after which iteration progress should be displayed |
use_gibbs | Logical indicating whether Gibbs sampling instead of
Metropolis-Hastings sampling should be used. Can be only
applied for |
aggregation | Logical indicating whether complete dataset should be
used for computing the predictive distribution of missing values.
|
object | Object of corresponding class |
x | Object of corresponding class |
digits | Number of digits in |
file | File to which the |
idparm | Indices for parameters to be plotted |
ask | Logical indicating whether the user is asked before new plot |
as_mids | Logical indicating whether multiply imputed datasets
should be converted into objects of class |
... | Further arguments to be passed |
The function allows for fitting a factored regression model. Consider the case of three variables \(Y\), \(X\) and \(Z\). A factored regression model consists of a sequence of univariate conditional models \(P(Y|X,Z)\), \(P(X|Z)\) and \(P(Z)\) such that the joint distribution can be factorized as $$ P(Y,X,Z)=P( Y|X,Z) P(X|Z) P(Z) $$ Each of the three variables can contain missing values. Missing values are integrated out by posing a distributional assumption for each variable with missing values.
For frm_em
it is a list containing the following values
Estimated coefficients
Covariance matrix
Summary parameter table
List with all estimated coefficients
Log likelihood value
Individual likelihood
Data frame with included latent values for each variable with missing values
Standard errors for coefficients
Information matrix
Convergence indicator
Number of iterations
Information criteria
List with model specifications including dep
and ind
Predictor matrix
Matrix containing all variables appearing in statistical models
Descriptive statistics of variables
Results from fitted models
Summary table with informations about MCMC algorithm
Sampled parameter values saved as class mcmc
for
analysis in coda package
Object containing informations of sampled parameters
Object containing informations of imputed datasets
Allison, P. D. (2012). Handling missing data by maximum likelihood. SAS Global Forum 2012.
Bartlett, J. W., & Morris, T. P. (2015) Multiple imputation of covariates by substantive-model compatible fully conditional specification. Stata Journal, 15(2), 437-456.
Bartlett, J. W., Seaman, S. R., White, I. R., Carpenter, J. R., & Alzheimer's Disease Neuroimaging Initiative (2015). Multiple imputation of covariates by fully conditional specification: Accommodating the substantive model. Statistical Methods in Medical Research, 24(4), 462-487. doi: 10.1177/0962280214521348
Erler, N. S., Rizopoulos, D., Rosmalen, J. V., Jaddoe, V. W., Franco, O. H., & Lesaffre, E. M. (2016). Dealing with missing covariates in epidemiologic studies: A comparison between multiple imputation and a full Bayesian approach. Statistics in Medicine, 35(17), 2955-2974. doi: 10.1002/sim.6944
Ibrahim, J. G. (1990). Incomplete data in generalized linear models. Journal of the American Statistical Association, 85(411), 765-769. doi: 10.1080/01621459.1990.10474938
Ibrahim, J. G., Chen, M. H., Lipsitz, S. R., & Herring, A. H. (2005). Missing-data methods for generalized linear models: A comparative review. Journal of the American Statistical Association, 100(469), 332-346. doi: 10.1198/016214504000001844
Jamshidian, M., & Jennrich, R. I. (2000). Standard errors for EM estimation. Journal of the Royal Statistical Society (Series B), 62(2), 257-270. doi: 10.1111/1467-9868.00230
Keller, B. T., & Enders, C. K. (2018). Blimp user's manual. Los Angeles, CA.
http://www.appliedmissingdata.com/multilevel-imputation.html
Lee, M. C., & Mitra, R. (2016). Multiply imputing missing values in data sets with mixed measurement scales using a sequence of generalised linear models. Computational Statistics & Data Analysis, 95(24), 24-38. doi: 10.1016/j.csda.2015.08.004
Luedtke, O., Robitzsch, A., & West, S. (2020a). Analysis of interactions and nonlinear effects with missing data: A factored regression modeling approach using maximum likelihood estimation. Multivariate Behavioral Research, 55(3), 361-381. doi: 10.1080/00273171.2019.1640104
Luedtke, O., Robitzsch, A., & West, S. (2020b). Regression models involving nonlinear effects with missing data: A sequential modeling approach using Bayesian estimation. Psychological Methods, 25(2), 157-181. doi: 10.1037/met0000233
Zhang, Q., & Wang, L. (2017). Moderation analysis with missing data in the predictors. Psychological Methods, 22(4), 649-666. doi: 10.1037/met0000104
The coef
and vcov
methods can be used to extract coefficients and
the corresponding covariance matrix, respectively. Standard errors for a fitted
object mod
can be extracted by making use of the survey
package and the statement survey::SE(mod)
.
See also the icdGLM package for estimation of generalized linear models with incomplete discrete covariates.
The imputation of covariates in substantive models with interactions or nonlinear terms can be also conducted with the JointAI package which is a wrapper to the JAGS software (see Erler et al., 2016). This package is also based on a sequential modelling approach.
The jomo package also accommodates substantive models (jomo::jomo.lm
)
based on a joint modeling framework.
Substantive model compatible imputation based on fully conditional specification can be found in the smcfcs package (see Bartlett et al., 2015; Bartlett & Morris, 2015) or the Blimp stand-alone software (Keller & Enders, 2018).
if (FALSE) { ############################################################################# # EXAMPLE 1: Simulated example linear regression with interaction effects ############################################################################# # The interaction model stats::lm( Y ~ X + Z + X:Z) is of substantive interest. # There can be missing values in all variables. data(data.mb01) dat <- data.mb01$missing #****************************************** # Model 1: ML approach #--- specify models # define integration nodes xnodes <- seq(-4,4,len=11) # nodes for X ynodes <- seq(-10,10,len=13) # nodes for Y. These ynodes are not really necessary for this dataset because # Y has no missing values. # define model for dependent variable Y dep <- list("model"="linreg", "formula"=Y ~ X*Z, "nodes"=ynodes ) # model P(X|Z) ind_X <- list( "model"="linreg", "formula"=X ~ Z, "nodes"=xnodes ) # all models for covariates ind <- list( "X"=ind_X ) #--- estimate model with numerical integration mod1 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind ) summary(mod1) # extract some informations coef(mod1) vcov(mod1) logLik(mod1) #****************************************** # Model 2: Fully Bayesian approach / Multiple Imputation #--- define models dep <- list("model"="linreg", "formula"=Y ~ X*Z ) ind_X <- list( "model"="linreg", "formula"=X ~ Z ) ind_Z <- list( "model"="linreg", "formula"=Z ~ 1 ) ind <- list( "X"=ind_X, Z=ind_Z) #--- estimate model mod2 <- mdmb::frm_fb(dat, dep, ind, burnin=200, iter=1000) summary(mod2) #* plot parameters plot(mod2) #--- create list of multiply imputed datasets datlist <- mdmb::frm2datlist(mod2) # convert into object of class mids imp2 <- miceadds::datlist2mids(datlist) # estimate linear model on multiply imputed datasets mod2c <- with(imp2, stats::lm( Y ~ X*Z ) ) summary( mice::pool(mod2c) ) #****************************************** # Model 3: Multiple imputation in jomo package library(jomo) # impute with substantive model containing interaction effects formula <- Y ~ X*Z imp <- jomo::jomo.lm( formula=formula, data=dat, nburn=100, nbetween=100) # convert to object of class mids datlist <- miceadds::jomo2mids( imp ) # estimate linear model mod3 <- with(datlist, lm( Y ~ X*Z ) ) summary( mice::pool(mod3) ) ############################################################################# # EXAMPLE 2: Simulated example logistic regression with interaction effects ############################################################################# # Interaction model within a logistic regression Y ~ X + Z + X:Z # Y and Z are dichotomous variables. # attach data data(data.mb02) dat <- data.mb02$missing #****************************************** # Model 1: ML approach #--- specify model # define nodes xnodes <- seq(-5,5,len=15) # X - normally distributed variable ynodes <- c(0,1) # Y and Z dichotomous variable # model P(Y|X,Z) for dependent variable dep <- list("model"="logistic", "formula"=Y ~ X*Z, "nodes"=ynodes ) # model P(X|Z) ind_X <- list( "model"="linreg", "formula"=X ~ Z, "nodes"=xnodes ) # model P(Z) ind_Z <- list( "model"="logistic", "formula"=Z ~ 1, "nodes"=ynodes ) ind <- list( "Z"=ind_Z, "X"=ind_X ) #--- estimate model with numerical integration mod1 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind ) summary(mod1) #****************************************** # Model 2: Fully Bayesian approach #--- specify model dep <- list("model"="logistic", "formula"=Y ~ X*Z ) ind_X <- list( "model"="linreg", "formula"=X ~ Z ) ind_Z <- list( "model"="logistic", "formula"=Z ~ 1 ) ind <- list( "Z"=ind_Z, "X"=ind_X ) #--- Bayesian estimation mod2 <- mdmb::frm_fb(dat=dat, dep=dep, ind=ind, burnin=500, iter=1000 ) summary(mod2) ############################################################################# # EXAMPLE 3: Confirmatory factor analysis ############################################################################# # A latent variable can be considered as missing data and the 'frm_em' function # is used to estimate the latent variable model. #--- simulate data N <- 1000 set.seed(91834) # latent variable theta <- stats::rnorm(N) # simulate items y1 <- 1.5 + 1*theta + stats::rnorm(N, sd=.7 ) y2 <- 1.9 + .7*theta + stats::rnorm(N, sd=1 ) y3 <- .9 + .7*theta + stats::rnorm(N, sd=.2 ) dat <- data.frame(y1,y2,y3) dat$theta <- NA #****************************************** # Model 1: ML approach #--- define model nodes <- seq(-4,4,len=21) ind_y1 <- list("model"="linreg", "formula"=y1 ~ offset(1*theta), "nodes"=nodes ) ind_y2 <- list( "model"="linreg", "formula"=y2 ~ theta, "nodes"=nodes, "coef_inits"=c(NA,1) ) ind_y3 <- list( "model"="linreg", "formula"=y3 ~ theta, "nodes"=nodes, "coef_inits"=c(1,1) ) dep <- list( "model"="linreg", "formula"=theta ~ 0, "nodes"=nodes ) ind <- list( "y1"=ind_y1, "y2"=ind_y2, "y3"=ind_y3) #*** estimate model with mdmb::frm_em mod1 <- mdmb::frm_em(dat, dep, ind) summary(mod1) #*** estimate model in lavaan library(lavaan) lavmodel <- " theta=~ 1*y1 + y2 + y3 theta ~~ theta " mod1b <- lavaan::cfa( model=lavmodel, data=dat ) summary(mod1b) # compare likelihood logLik(mod1) logLik(mod1b) ############################################################################# # EXAMPLE 4: Rasch model ############################################################################# #--- simulate data set.seed(91834) N <- 500 # latent variable theta0 <- theta <- stats::rnorm(N) # number of items I <- 7 dat <- sirt::sim.raschtype( theta, b=seq(-1.5,1.5,len=I) ) colnames(dat) <- paste0("I",1:I) dat$theta <- NA #****************************************** # Model 1: ML approach #--- define model nodes <- seq(-4,4,len=13) dep <- list("model"="linreg", "formula"=theta ~ 0, "nodes"=nodes ) ind <- list() for (ii in 1:I){ ind_ii <- list( "model"="logistic", formula= stats::as.formula( paste0("I",ii, " ~ offset(1*theta)") ) ) ind[[ii]] <- ind_ii } names(ind) <- colnames(dat)[-(I+1)] #--- estimate Rasch model with mdmb::frm_em mod1 <- mdmb::frm_em(dat, dep, ind ) summary(mod1) #--- estmate Rasch model with sirt package library(sirt) mod2 <- sirt::rasch.mml2( dat[,-(I+1)], theta.k=nodes, use.freqpatt=FALSE) summary(mod2) #** compare estimated parameters round(cbind(coef(mod1), c( mod2$sd.trait, -mod2$item$thresh[ seq(I,1)] ) ), 3) ############################################################################# # EXAMPLE 5: Regression model with measurement error in covariates ############################################################################# #--- simulate data set.seed(768) N <- 1000 # true score theta <- stats::rnorm(N) # heterogeneous error variance var_err <- stats::runif(N, .5, 1) # simulate observed score x <- theta + stats::rnorm(N, sd=sqrt(var_err) ) # simulate outcome y <- .3 + .7 * theta + stats::rnorm( N, sd=.8 ) dat0 <- dat <- data.frame( y=y, x=x, theta=theta ) #*** estimate model with true scores (which are unobserved in real datasets) mod0 <- stats::lm( y ~ theta, data=dat0 ) summary(mod0) #****************************************** # Model 1: Model-based approach #--- specify model dat$theta <- NA nodes <- seq(-4,4,len=15) dep <- list( "model"="linreg", "formula"=y ~ theta, "nodes"=nodes, "coef_inits"=c(NA, .4 ) ) ind <- list() ind[["theta"]] <- list( "model"="linreg", "formula"=theta ~ 1, "nodes"=nodes ) ind[["x"]] <- list( "model"="linreg", "formula"=x ~ 0 + offset(theta), "nodes"=nodes ) # assumption of heterogeneous known error variance ind[["x"]]$sigma_fixed <- sqrt( var_err ) #--- estimate regression model mod1 <- mdmb::frm_em(dat, dep, ind ) summary(mod1) #****************************************** # Model 2: Fully Bayesian estimation #--- specify model dep <- list( "model"="linreg", "formula"=y ~ theta ) ind <- list() ind[["theta"]] <- list( "model"="linreg", "formula"=theta ~ 1 ) ind[["x"]] <- list( "model"="linreg", "formula"=x ~ 0 + offset(theta) ) # assumption of heterogeneous known error variance ind[["x"]]$sigma_fixed <- sqrt( var_err ) data_init <- dat data_init$theta <- dat$x # estimate model mod2 <- mdmb::frm_fb(dat, dep, ind, burnin=200, iter=1000, data_init=data_init) summary(mod2) plot(mod2) ############################################################################# # EXAMPLE 6: Non-normally distributed covariates: # Positive values with Box-Cox transformation ############################################################################# # simulate data with chi-squared distributed covariate from # regression model Y ~ X set.seed(876) n <- 1500 df <- 2 x <- stats::rchisq( n, df=df ) x <- x / sqrt( 2*df) y <- 0 + 1*x R2 <- .25 # explained variance y <- y + stats::rnorm(n, sd=sqrt( (1-R2)/R2 * 1) ) dat0 <- dat <- data.frame( y=y, x=x ) # simulate missing responses prop_miss <- .5 cor_miss <- .7 resp_tend <- cor_miss*(dat$y-mean(y) )/ stats::sd(y) + stats::rnorm(n, sd=sqrt( 1 - cor_miss^2) ) dat[ resp_tend < stats::qnorm( prop_miss ), "x" ] <- NA summary(dat) #-- complete data mod0 <- stats::lm( y ~ x, data=dat0 ) summary(mod0) #-- listwise deletion mod1 <- stats::lm( y ~ x, data=dat ) summary(mod1) # normal distribution assumption for frm # define models dep <- list("model"="linreg", "formula"=y ~ x ) # normally distributed data ind_x1 <- list( "model"="linreg", "formula"=x ~ 1 ) # Box-Cox normal distribution ind_x2 <- list( "model"="bctreg", "formula"=x ~ 1, nodes=c( seq(0.05, 3, len=31), seq( 3.5, 9, by=.5 ) ) ) ind1 <- list( "x"=ind_x1 ) ind2 <- list( "x"=ind_x2 ) #--- incorrect normal distribution assumption mod1 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind1 ) summary(mod1) #--- model chi-square distribution of predictor with Box-Cox transformation mod2 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind2 ) summary(mod2) ############################################################################# # EXAMPLE 7: Latent interaction model ############################################################################# # A latent interaction model Y ~ FX + FZ is of interest. Y is directly observed, # FX and FZ are both indirectly observed by three items #--- simulate data N <- 1000 set.seed(987) # latent variable FX <- stats::rnorm(N) FZ <- stats::rnorm(N) # simulate items x1 <- 1.5 + 1*FX + stats::rnorm(N, sd=.7 ) x2 <- 1.9 + .7*FX + stats::rnorm(N, sd=1 ) x3 <- .9 + .7*FX + stats::rnorm(N, sd=.2 ) z1 <- 1.5 + 1*FZ + stats::rnorm(N, sd=.7 ) z2 <- 1.9 + .7*FZ + stats::rnorm(N, sd=1 ) z3 <- .9 + .7*FZ + stats::rnorm(N, sd=.2 ) dat <- data.frame(x1,x2,x3,z1,z2,z3) dat$FX <- NA dat$FZ <- NA dat$y <- 2 + .5*FX + .3*FZ + .4*FX*FZ + rnorm( N, sd=1 ) # estimate interaction model with ML #--- define model nodes <- seq(-4,4,len=11) ind_x1 <- list("model"="linreg", "formula"=x1 ~ offset(1*FX), "nodes"=nodes ) ind_x2 <- list( "model"="linreg", "formula"=x2 ~ FX, "nodes"=nodes, "coef_inits"=c(NA,1) ) ind_x3 <- list( "model"="linreg", "formula"=x3 ~ FX, "nodes"=nodes, "coef_inits"=c(1,1) ) ind_FX <- list( "model"="linreg", "formula"=FX ~ 0, "nodes"=nodes ) ind_z1 <- list("model"="linreg", "formula"=z1 ~ offset(1*FZ), "nodes"=nodes ) ind_z2 <- list( "model"="linreg", "formula"=z2 ~ FZ, "nodes"=nodes, "coef_inits"=c(NA,1) ) ind_z3 <- list( "model"="linreg", "formula"=z3 ~ FZ, "nodes"=nodes, "coef_inits"=c(1,1) ) ind_FZ <- list( "model"="linreg", "formula"=FZ ~ 0 + FX, "nodes"=nodes ) ind <- list( "x1"=ind_x1, "x2"=ind_x2, "x3"=ind_x3, "FX"=ind_FX, "z1"=ind_z1, "z2"=ind_z2, "z3"=ind_z3, "FX"=ind_FZ ) dep <- list( "model"="linreg", formula=y ~ FX+FZ+FX*FZ, "coef_inits"=c(1,.2,.2,0) ) #*** estimate model with mdmb::frm_em mod1 <- mdmb::frm_em(dat, dep, ind) summary(mod1) ############################################################################# # EXAMPLE 8: Non-ignorable data in Y ############################################################################# # regression Y ~ X in which Y is missing depending Y itself library(mvtnorm) cor_XY <- .4 # correlation between X and Y prop_miss <- .5 # missing proportion cor_missY <- .7 # correlation with missing propensity N <- 3000 # sample size #----- simulate data set.seed(790) Sigma <- matrix( c(1, cor_XY, cor_XY, 1), 2, 2 ) mu <- c(0,0) dat <- mvtnorm::rmvnorm( N, mean=mu, sigma=Sigma ) colnames(dat) <- c("X","Y") dat <- as.data.frame(dat) #-- generate missing responses on Y depending on Y itself y1 <- dat$Y miss_tend <- cor_missY * y1 + rnorm(N, sd=sqrt( 1 - cor_missY^2) ) dat$Y[ miss_tend < quantile( miss_tend, prop_miss ) ] <- NA #--- ML estimation under assumption of ignorability nodes <- seq(-5,5,len=15) dep <- list("model"="linreg", "formula"=Y ~ X, "nodes"=nodes ) ind_X <- list( "model"="linreg", "formula"=X ~ 1, "nodes"=nodes ) ind <- list( "X"=ind_X ) mod1 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind) summary(mod1) #--- ML estimation under assumption with specifying a model for non-ignorability # for response indicator resp_Y dat$resp_Y <- 1* ( 1 - is.na(dat$Y) ) dep <- list("model"="linreg", "formula"=Y ~ X, "nodes"=nodes ) ind_X <- list( "model"="linreg", "formula"=X ~ 1, "nodes"=nodes ) ind_respY <- list( "model"="logistic", "formula"=resp_Y ~ Y, "nodes"=nodes ) ind <- list( "X"=ind_X, "resp_Y"=ind_respY ) mod2 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind) summary(mod2) ############################################################################# # EXAMPLE 9: Ordinal variables: Graded response model ############################################################################# #--- simulate data N <- 2000 set.seed(91834) # latent variable theta <- stats::rnorm(N) # simulate items y1 <- 1*theta + stats::rnorm(N) y2 <- .7*theta + stats::rnorm(N) y3 <- .7*theta + stats::rnorm(N) # discretize variables y1 <- as.numeric( cut( y1, breaks=c(-Inf, -.5, 0.4, 1, Inf ) ) ) - 1 y2 <- as.numeric( cut( y2, breaks=c(-Inf, 0.3, 1, Inf ) ) ) - 1 y3 <- as.numeric( cut( y3, breaks=c(-Inf, .2, Inf ) ) ) - 1 # define dataset dat <- data.frame(y1,y2,y3) dat$theta <- NA #****************************************** # Model 1: Fully Bayesian estimation #--- define model ind_y1 <- list( "model"="oprobit", "formula"=y1 ~ offset(1*theta) ) ind_y2 <- list( "model"="oprobit", "formula"=y2 ~ theta ) ind_y3 <- list( "model"="oprobit", "formula"=y3 ~ theta ) dep <- list( "model"="linreg", "formula"=theta ~ 0 ) ind <- list( "y1"=ind_y1, "y2"=ind_y2, "y3"=ind_y3) # initial data data_init <- dat data_init$theta <- as.numeric( scale(dat$y1) ) + stats::rnorm(N, sd=.4 ) #-- estimate model iter <- 3000; burnin <- 1000 mod1 <- mdmb::frm_fb(dat=dat, dep=dep, ind=ind, data_init=data_init, iter=iter, burnin=burnin) summary(mod1) plot(mod1) ############################################################################# # EXAMPLE 10: Imputation for missig predictors in models with interaction # effects in multilevel regression models ############################################################################# library(miceadds) data(data.mb04, package="mdmb") dat <- data.mb04 #*** model specification mcmc_iter <- 4 # number of MCMC iterations for model parameter sampling model_formula <- y ~ cwc(x, idcluster) + gm(x, idcluster) + w + w*cwc(x, idcluster) + w*gm(x, idcluster) + ( 1 + cwc(x, idcluster) | idcluster) dep <- list("model"="mlreg", "formula"=model_formula, R_args=list(iter=mcmc_iter, outcome="normal") ) ind_x <- list( "model"="mlreg", "formula"=x ~ w + (1|idcluster), R_args=list(iter=mcmc_iter), sampling_level="idcluster" ) # group means of x are involved in the outcome model. Therefore, Metropolis-Hastings # sampling of missing values in x should be conducted at the level of clusters, # i.e. specifying sampling_level ind <- list("x"=ind_x) # --- estimate model mod1 <- mdmb::frm_fb(dat, dep, ind, aggregation=TRUE) # argument aggregation is necessary because group means are involved in regression formulas #------------- #*** imputation of a continuous level-2 variable w # create artificially some missings on w dat[ dat$idcluster %%3==0, "w" ] <- NA # define level-2 model with argument variable_level ind_w <- list( "model"="linreg", "formula"=w ~ 1, "variable_level"="idcluster" ) ind <- list( x=ind_x, w=ind_w) #* conduct imputations mod2 <- mdmb::frm_fb(dat, dep, ind, aggregation=TRUE) summary(mod2) #--- Model 1 with user-defined prior distributions for covariance matrices model_formula <- y ~ cwc(x, idcluster) + gm(x, idcluster) + w + w*cwc(x, idcluster) + w*gm(x, idcluster) + ( 1 + cwc(x, idcluster) | idcluster) # define scale degrees of freedom (nu) and scale matrix (S) for inverse Wishart distribution psi_nu0_list <- list( -3 ) psi_S0_list <- list( diag(0,2) ) dep <- list("model"="mlreg", "formula"=model_formula, R_args=list(iter=mcmc_iter, outcome="normal", psi_nu0_list=psi_nu0_list, psi_S0_list=psi_S0_list ) ) # define nu and S parameters for covariate model psi_nu0_list <- list( .4 ) psi_S0_list <- list( matrix(.2, nrow=1, ncol=1) ) ind_x <- list( "model"="mlreg", "formula"=x ~ w + (1|idcluster), R_args=list(iter=mcmc_iter, psi_nu0_list=psi_nu0_list, psi_S0_list=psi_S0_list), sampling_level="idcluster" ) ind <- list("x"=ind_x) # --- estimate model mod3 <- mdmb::frm_fb(dat, dep, ind, aggregation=TRUE) ############################################################################# # EXAMPLE 11: Bounded variable combined with Yeo-Johnson transformation ############################################################################# #*** simulate data set.seed(876) n <- 1500 x <- mdmb::ryjt_scaled( n, location=-.2, shape=.8, lambda=.9, probit=TRUE) R2 <- .25 # explained variance y <- 1*x + stats::rnorm(n, sd=sqrt( (1-R2)/R2 * stats::var(x)) ) dat0 <- dat <- data.frame( y=y, x=x ) # simulate missing responses prop_miss <- .5 cor_miss <- .7 resp_tend <- cor_miss*(dat$y-mean(y) )/ stats::sd(y) + stats::rnorm(n, sd=sqrt( 1 - cor_miss^2) ) dat[ resp_tend < stats::qnorm(prop_miss), "x" ] <- NA summary(dat) #*** define models dep <- list("model"="linreg", "formula"=y ~ x ) # distribution according to Yeo-Johnson transformation ind_x1 <- list( "model"="yjtreg", "formula"=x ~ 1 ) # distribution according to Probit Yeo-Johnson transformation ind_x2 <- list( "model"="yjtreg", "formula"=x ~ 1, R_args=list("probit"=TRUE ) ) ind1 <- list( "x"=ind_x1 ) ind2 <- list( "x"=ind_x2 ) #--- complete data mod0 <- stats::lm( y~x, data=dat0) summary(mod0) #--- Yeo-Johnson normal distribution (for unbounded variables) mod1 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind1 ) summary(mod1) #--- Probit Yeo-Johnson normal distribution (for bounded variable on (0,1)) mod2 <- mdmb::frm_em(dat=dat, dep=dep, ind=ind2) summary(mod2) #--- same model, but MCMC estimation mod3 <- mdmb::frm_fb(dat, dep, ind=ind2, burnin=2000, iter=5000) summary(mod3) plot(mod3) ############################################################################# # EXAMPLE 12: Yeo-Johnson transformation with estimated degrees of freedom ############################################################################# #*** simulate data set.seed(876) n <- 1500 x <- mdmb::ryjt_scaled( n, location=-.2, shape=.8, lambda=.9, df=10 ) R2 <- .25 # explained variance y <- 1*x + stats::rnorm(n, sd=sqrt( (1-R2)/R2 * stats::var(x)) ) dat0 <- dat <- data.frame( y=y, x=x ) # simulate missing responses prop_miss <- .5 cor_miss <- .7 resp_tend <- cor_miss*(dat$y-mean(y) )/ stats::sd(y) + stats::rnorm(n, sd=sqrt( 1-cor_miss^2) ) dat[ resp_tend < stats::qnorm(prop_miss), "x" ] <- NA summary(dat) #*** define models dep <- list("model"="linreg", "formula"=y ~ x ) # specify distribution with estimated degrees of freedom ind_x <- list( "model"="yjtreg", "formula"=x ~ 1, R_args=list(est_df=TRUE ) ) ind <- list( "x"=ind_x ) #--- Yeo-Johnson t distribution mod1 <- mdmb::frm_fb(dat=dat, dep=dep, ind=ind, iter=3000, burnin=1000 ) summary(mod1) }