Package 'misty'

Description

This function performs an one-way between-subject analysis of variance (ANOVA) including Tukey HSD post hoc tests for multiple comparison and provides descriptive statistics, effect size measures, and a plot showing bars representing means for each group and error bars for difference-adjusted confidence intervals.

Usage

aov.b(formula, data, hypo = FALSE, descript = FALSE, effsize = FALSE,
      weighted = TRUE, correct = FALSE, posthoc = FALSE, conf.level = 0.95,
      digits = 2, p.digits = 3, as.na = NULL, plot = FALSE, bar = TRUE,
      point = FALSE, ci = TRUE, jitter = FALSE, adjust = TRUE,
      filename = NULL, width = NA, height = NA, dpi = 600,
      write = NULL, append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

Confidence Intervals

Cumming and Finch (2005) pointed out that when 95% confidence intervals (CI) for two separately plotted means overlap, it is still possible that the CI for the difference would not include zero. Baguley (2012) proposed to adjust the width of the CIs by the factor of $\sqrt{2}$ to reflect the correct width of the CI for a mean difference:

$\hat{\mu}_{j} \pm t_{n - 1, 1 - \alpha/2} \frac{\sqrt{2}}{2} \hat{\sigma}_{{\hat{\mu}}_j}$

These difference-adjusted CIs around the individual means can be interpreted as if it were a CI for their difference. Note that the width of these intervals is sensitive to differences in the variance and sample size of each sample, i.e., unequal population variances and unequal $n$ alter the interpretation of difference-adjusted CIs.

Value

Returns an object of class misty.object, which is a list with following entries:

Author(s)

Takuya Yanagida [email protected]

References

Baguley, T. S. (2012a). Serious stats: A guide to advanced statistics for the behavioral sciences. Palgrave Macmillan.

Cumming, G., and Finch, S. (2005) Inference by eye: Confidence intervals, and how to read pictures of data. American Psychologist, 60, 170–80.

Rasch, D., Kubinger, K. D., & Yanagida, T. (2011). Statistics in psychology - Using R and SPSS. John Wiley & Sons.

See Also

aov.w, test.t, test.z, test.levene, aov.b, cohens.d, ci.mean.diff, ci.mean

Examples

#————————————————————————————————————————————————————————————————————————————
# Between-Subject Analysis of Variance

# Example 1a: Between-Subject ANOVA
aov.b(hp ~ gear, data = mtcars)

# Example 1b: Between-Subject ANOVA
# Print descriptive statistics and Tukey HSD post hoc test
aov.b(hp ~ gear, data = mtcars, descript = TRUE, posthoc = TRUE)

# Example 1c: Between-Subject ANOVA, print eta-squared and omega-squared
aov.b(hp ~ gear, data = mtcars, effsize = TRUE)

#————————————————————————————————————————————————————————————————————————————
# Plot

# Example 2a: Plot results, default setting
aov.b(hp ~ gear, data = mtcars, plot = TRUE)

# Example 2b: Plot results
# No bars, draw points representing means and jittered data points
aov.b(hp ~ gear, data = mtcars, plot = TRUE, bar = FALSE, point = TRUE, jitter = TRUE)

# Example 2c: Plot results using the plot() function, use additional arguments
# see Details in the help page of the function plot.misty.object
object <- aov.b(hp ~ gear, data = mtcars)
plot(object, jitter = TRUE, jitter.alpha = 0.4, title = "Between-Subject ANOVA")

#————————————————————————————————————————————————————————————————————————————
# Create Plot Manually

# Load ggplot2 package
library(ggplot2)

# Create misty object
object <- aov.b(hp ~ gear, data = mtcars)

# Example 3: Plot
ggplot(object$result$descript, aes(group, y)) +
  geom_bar(aes(group, m), stat = "summary", fun = "mean") +
  geom_jitter(data = object$data, aes(group, y), alpha = 0.1, width = 0.05,
             height = 0, size = 1.25) +
  geom_point(aes(group, m), stat = "identity", size = 3) +
  geom_errorbar(aes(group, m, ymin = low, ymax = upp), width = 0.1) +
  theme_bw()

#————————————————————————————————————————————————————————————————————————————
# Write Results and Save Plot

## Not run: 

# Example 4a: Write results into a text file
aov.b(hp ~ gear, data = mtcars, write = "ANOVA.txt")

# Example 4b: Write results into an Excel file
aov.b(hp ~ gear, data = mtcars, write = "ANOVA.xlsx")

# Example 4c: Save plot as PNG fine
aov.b(hp ~ gear, data = mtcars, plot = TRUE, filename = "ANOVA.png",
      width = 6, height = 5)

## End(Not run)

Description

This function performs an one-way repeated measures analysis of variance (within subject ANOVA) including paired-samples t-tests for multiple comparison and provides descriptive statistics, effect size measures, and a plot showing error bars for difference-adjusted Cousineau-Morey within-subject confidence intervals with jittered data points including subject-specific lines.

Usage

aov.w(formula, data, print = c("all", "none", "GG", "HF", "LB"),
      hypo = FALSE, epsilon = FALSE, descript = FALSE, effsize = FALSE,
      posthoc = FALSE, na.omit = TRUE, conf.level = 0.95,
      p.adj = c("none", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY", "fdr"),
      digits = 2, p.digits = 3, as.na = NULL, plot = FALSE, point = TRUE,
      line = TRUE, ci = TRUE, jitter = FALSE, adjust = TRUE, filename = NULL,
      width = NA, height = NA, dpi = 600, write = NULL, append = TRUE,
      check = TRUE, output = TRUE)

Arguments

Details

Sphericity

The F-Test of the repeated measures ANOVA is based on the assumption of sphericity, which is defined as the assumption that the variance of differences between repeated measures are equal in the population. The Mauchly's test is commonly used to test this hypothesis. However, test of assumptions addresses an irrelevant hypothesis because what matters is the degree of violation rather than its presence (Baguley, 2012a). Moreover, the test is not recommended because it lacks statistical power (Abdi, 2010). Instead, the Box index of sphericity ($\varepsilon$) should be used to assess the degree of violation of the sphericity assumption. The $\varepsilon$ parameter indicates the degree to which the population departs from sphericity with $\varepsilon = 1$ indicating that sphericity holds. As the departure becomes more extreme, $\varepsilon$ approaches its lower bound $\hat{\varepsilon}_{lb}$:

$\hat{\varepsilon}_{lb} = \frac{1}{J - 1}$

where $J$ is the number of levels of the within-subject factor. Box (1954a, 1954b) suggested a measure for sphericity, which applies to a population covariance matrix. Greenhouse and Geisser (1959) proposed an estimate for $\varepsilon$ known as $\hat{\varepsilon}_{gg}$ that can be computed from the sample covariance matrix, whereas Huynh and Feldt (1976) proposed an alternative estimate $\hat{\varepsilon}_{hf}$. These estimates can be used to correct the effect and error df of the F-test. Simulation studies showed that $\hat{\varepsilon}_{gg} \leq \hat{\varepsilon}_{hf}$ and that $\hat{\varepsilon}_{gg}$ tends to be conservative underestimating $\varepsilon$, whereas $\hat{\varepsilon}_{hf}$ tends to be liberal overestimating $\varepsilon$ and occasionally exceeding one. Baguley (2012a) recommended to compute the average of the conservative estimate $\hat{\varepsilon}_{gg}$ and the liberal estimate $\hat{\varepsilon}_{hf}$ to assess the sphericity assumption. By default, the function prints results depending on the average $\hat{\varepsilon}_{gg}$ and $\hat{\varepsilon}_{hf}$:

• If the average is less than 0.75 results of the F-Test based on Greenhouse-Geiser correction factor ($\hat{\varepsilon}_{gg}$) is printed.

• If the average is less greater or equal 0.75, but less than 0.95 results of the F-Test based on Huynh-Feldt correction factor ($\hat{\varepsilon}_{hf}$) is printed.

• If the average is greater or equal 0.95 results of the F-Test without any corrections are printed.

Missing Data

The function uses listwise deletion by default to deal with missing data. However, the function also allows to use all available observations by conducting the repeated measures ANOVA in long data format when specifying na.omit = FALSE. Note that in the presence of missing data, the F-Test without any sphericity corrections may be reliable, but it is not clear whether results based on Greenhouse-Geiser or Huynh-Feldt correction are trustworthy given that pairwise deletion is used for estimating the variance-covariance matrix when computing $\hat{\varepsilon}_{gg}$ and the total number of subjects regardless of missing values (i.e., complete and incomplete cases) are used for computing $\hat{\varepsilon}_{hf}$.

Within-Subject Confidence Intervals

The function provides a plot showing error bars for difference-adjusted Cousineau-Morey confidence intervals (Baguley, 2012b). The intervals matches that of a CI for a difference, i.e., non-overlapping CIs corresponds to an inferences of no statistically significant difference. The Cousineau-Morey confidence intervals without adjustment can be used by specifying adjust = FALSE.

Value

Returns an object of class misty.object, which is a list with following entries:

Author(s)

Takuya Yanagida [email protected]

References

Abdi, H. (2010). The Greenhouse-Geisser correction. In N. J. Salkind (Ed.) Encyclopedia of Research Design (pp. 630-634), Sage. https://dx.doi.org/10.4135/9781412961288

Baguley, T. S. (2012a). Serious stats: A guide to advanced statistics for the behavioral sciences. Palgrave Macmillan.

Baguley, T. (2012b). Calculating and graphing within-subject confidence intervals for ANOVA. Behavior Research Methods, 44, 158-175. https://doi.org/10.3758/s13428-011-0123-7

Bakerman, R. (2005). Recommended effect size statistics for repeated measures designs. Behavior Research Methods, 37, 179-384. https://doi.org/10.3758/BF03192707

Box, G. E. P. (1954a) Some Theorems on Quadratic Forms Applied in the Study of Analysis of Variance Problems, I. Effects of Inequality of Variance in the One-way Classification. Annals of Mathematical Statistics, 25, 290–302.

Box, G. E. P. (1954b) Some Theorems on Quadratic Forms Applied in the Study of Analysis of Variance Problems, II. Effects of Inequality of Variance and of Correlation between Errors in the Two-way Classification. Annals of Mathematical Statistics, 25, 484–98.

Greenhouse, S. W., and Geisser, S. (1959). On methods in the analysis of profile data.Psychometrika, 24, 95–112. https://doi.org/10.1007/BF02289823

Huynh, H., and Feldt, L. S. (1976). Estimation of the box correction for degrees of freedom from sample data in randomized block and splitplot designs. Journal of Educational Statistics, 1, 69–82. https://doi.org/10.2307/1164736

Olejnik, S., & Algina, J. (2000). Measures of effect size for comparative studies: Applications, interpretations, and limitations. Contemporary Educational Psychology, 25, 241–286. https://doi.org/10.1006/ceps.2000.1040

Rasch, D., Kubinger, K. D., & Yanagida, T. (2011). Statistics in psychology - Using R and SPSS. John Wiley & Sons.

See Also

aov.b, test.t, test.z, cohens.d, ci.mean.diff, ci.mean

Examples

#————————————————————————————————————————————————————————————————————————————
# Within-Subject Analysis of Variance

# Data frame
dat <- data.frame(time1 = c(3, 2, 1, 4, 5, 2, 3, 5, 6, 7),
                  time2 = c(4, 3, 6, 5, 8, 6, 7, 3, 4, 5),
                  time3 = c(1, 2, 2, 3, 6, 5, 1, 2, 4, 6))

# Example 1a: Within-Subject ANOVA
aov.w(cbind(time1, time2, time3) ~ 1, data = dat)

# Example 1b: Within-Subject ANOVA
# Print descriptive statistics and Paired-samples t-tests
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, descript = TRUE, posthoc = TRUE)

# Example 1c: Within-Subject ANOVA, print eta-squared and omega-squared
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, effsize = TRUE)

#————————————————————————————————————————————————————————————————————————————
# Plot

# Example 2a: Plot results, default setting
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, plot = TRUE)

# Example 2b: Plot results
# No bars, jittered data points without subject-specific lines
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, plot = TRUE,
      line = FALSE, jitter = TRUE)

# Example 2c: Plot results using the plot() function, use additional arguments
# see Details in the help page of the function plot.misty.object
object <- aov.w(cbind(time1, time2, time3) ~ 1, data = dat)
plot(object, jitter = TRUE, jitter.alpha = 0.4, title = "Within-Subject ANOVA")

#————————————————————————————————————————————————————————————————————————————
# Create Plot Manually

# Load ggplot2 package
library(ggplot2)

# Create misty object
object <- aov.w(cbind(time1, time2, time3) ~ 1, data = dat)

# Compute Means and difference-adjusted confidence intervals
ci.table <- ci.mean.w(object$data$wide, adjust = TRUE, output = FALSE)$result

# Example 3: Plot
ggplot(object$data$long, aes(time, y, group = 1)) +
  geom_errorbar(data = ci.table, aes(variable, m, ymin = low, ymax = upp),
                width = 0.1) +
  geom_point(data = ci.table, aes(variable, m), stat = "identity", size = 3) +
  geom_line(data = ci.table, aes(variable, m), stat = "identity") +
  geom_jitter(alpha = 0.2, width = 0.05, height = 0, size = 1.25) +
  theme_bw()

#————————————————————————————————————————————————————————————————————————————
# Write Results and Save Plot

## Not run: 

# Example 4a: Write results into a text file
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, write = "RM-ANOVA.txt")

# Example 4b: Write results into an Excel file
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, write = "RM-ANOVA.xlsx")

# Example 4c: Save plot as PNG fine
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, plot = TRUE,
            filename = "RM-ANOVA.png", width = 6, height = 5)

## End(Not run)

Description

This wrapper function creates a Blimp input file, runs the input file by using the blimp.run() function, and prints the Blimp output file by using the blimp.print() function.

Usage

blimp(x, file = "Blimp_Input.imp", data = NULL, comment = FALSE, replace.inp = TRUE,
      blimp.run = TRUE, posterior = FALSE, folder = "Posterior_",
      format = c("csv", "csv2", "excel", "rds", "workspace"), clear = TRUE,
      replace.out = c("always", "never", "modified"), Blimp = .detect.blimp(),
      result = c("all", "default", "algo.options", "data.info", "model.info",
                 "warn.mess", "error.mess", "out.model", "gen.param"),
      exclude = NULL, color = c("none", "blue", "green"),
      style = c("bold", "regular"), not.result = TRUE, write = NULL,
      append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

VARIABLES Section

The VARIABLES section used to assign names to the variables in the data set can be specified by using the data argument:

• Write Blimp Data File: In the first step, the Blimp data file is written by using the write.mplus() function, e.g. write.mplus(data1, file = "data1.dat").

• Specify Blimp Input: In the second step, the Blimp input is specified as a character string. The VARIABLES option is left out from the Blimp input text, e.g., input <- 'DATA: data1.dat;\nMODEL: y ~ x1@b1 x2@b2 d2;'.

• Run Blimp Input: In the third step, the Blimp input is run by using the blimp() function. The argument data needs to be specified given that the VARIABLES section was left out from the Blimp input text in the previous step, e.g., blimp(input, file = "Ex4.3.imp", data = data1).

Note that unlike Mplus, Blimp allows to specify a CSV data file with variable names in the first row. Hence, it is recommended to export the data from R using the write.csv() function to specify the data file in the DATA section of the Blimp input file without specifying the VARIABLES section.

Value

Returns an object of class misty.object, which is a list with following entries:

Author(s)

Takuya Yanagida

References

Keller, B. T., & Enders, C. K. (2023). Blimp user’s guide (Version 3). Retrieved from www.appliedmissingdata.com/blimp

See Also

blimp.update, blimp.run, blimp.print, blimp.plot, blimp.bayes

Examples

## Not run: 

#————————————————————————————————————————————————————————————————————————————
# Example 1: Write data, specify input without VARIABLES section, and run input

# Write Data File
# Note that row.names = FALSE needs to be specified
write.csv(data1, file = "data1.csv", row.names = FALSE)

# Specify Blimp input
input1 <- '
DATA: data1.csv;
ORDINAL: d;
MISSING: 999;
FIXED: d;
CENTER: x1 x2;
MODEL: y ~ x1 x2 d;
SEED: 90291;
BURN: 1000;
ITERATIONS: 10000;
'

# Run Blimp input
blimp(input1, file = "Ex4.3.imp")

#————————————————————————————————————————————————————————————————————————————
# Example 2: Write data, specify input with VARIABLES section, and run input

# Write Data File
write.mplus(data1, file = "data1.dat", input = FALSE)

# Specify Blimp input
input2 <- '
DATA: data1.dat;
VARIABLES: id v1 v2 v3 y x1 d x2 v4;
ORDINAL: d;
MISSING: 999;
FIXED: d;
CENTER: x1 x2;
MODEL: y ~ x1 x2 d;
SEED: 90291;
BURN: 1000;
ITERATIONS: 10000;
'

# Run Blimp input
blimp(input2, file = "Ex4.3.imp")

#————————————————————————————————————————————————————————————————————————————
# Example 3: Alternative specification using the data argument

# Write Data File
write.mplus(data1, file = "data1.dat", input = FALSE)

# Specify Blimp input
input3 <- '
DATA: data1.dat;
ORDINAL: d;
MISSING: 999;
FIXED: d;
CENTER: x1 x2;
MODEL: y ~ x1 x2 d;
SEED: 90291;
BURN: 1000;
ITERATIONS: 10000;
'

# Run Blimp input
blimp(input3, file = "Ex4.3.imp", data = data1)

## End(Not run)

Description

This function reads the posterior distribution for all parameters saved in long format in a file called posterior.* by the function blimp.run or blimp when specifying posterior = TRUE to compute point estimates (i.e., mean, median, and MAP), measures of dispersion (i.e., standard deviation and mean absolute deviation), measures of shape (i.e., skewness and kurtosis), credible intervals (i.e., equal-tailed intervals and highest density interval), convergence and efficiency diagnostics (i.e., potential scale reduction factor R-hat, effective sample size, and Monte Carlo standard error), probability of direction, and probability of being in the region of practical equivalence for the posterior distribution for each parameter. By default, the function computes the maximum of rank-normalized split-R-hat and rank normalized folded-split-R-hat, Bulk effective sample size (Bulk-ESS) for rank-normalized values using split chains, tail effective sample size (Tail-ESS) defined as the minimum of the effective sample size for 0.025 and 0.975 quantiles, the Bulk Monte Carlo standard error (Bulk-MCSE) for the median and Tail Monte Carlo standard error (Tail-MCSE) defined as the maximum of the MCSE for 0.025 and 0.975 quantiles.

Usage

blimp.bayes(x, param = NULL,
            print = c("all", "default", "m", "med", "map", "sd", "mad",
                      "skew", "kurt", "eti", "hdi",
                      "rhat", "b.ess", "t.ess", "b.mcse", "t.mcse"),
            m.bulk = FALSE, split = TRUE, rank = TRUE, fold = TRUE,
            pd = FALSE, null = 0, rope = NULL,
            ess.tail = c(0.025, 0.975), mcse.tail = c(0.025, 0.975),
            alternative = c("two.sided", "less", "greater"), conf.level = 0.95,
            digits = 2, r.digits = 3, ess.digits = 0, mcse.digits = 3,
            p.digits = 3, write = NULL, append = TRUE, check = TRUE,
            output = TRUE)

Arguments

Details

Convergence and Efficiency Diagnostics for Markov Chains

Convergence and efficiency diagnostics for Markov chains is based on following numeric measures:

• Potential Scale Reduction (PSR) factor R-hat: The PSR factor R-hat compares the between- and within-chain variance for a model parameter, i.e., R-hat larger than 1 indicates that the between-chain variance is greater than the within-chain variance and chains have not mixed well. According to the default setting, the function computes the improved R-hat as recommended by Vehtari et al. (2020) based on rank-normalizing (i.e., rank = TRUE) and folding (i.e., fold = TRUE) the posterior draws after splitting each MCMC chain in half (i.e., split = TRUE). The traditional R-hat used in Blimp can be requested by specifying split = TRUE, rank = FALSE, and fold = FALSE. Note that the traditional R-hat can catch many problems of poor convergence, but fails if the chains have different variances with the same mean parameter or if the chains have infinite variance with one of the chains having a different location parameter to the others (Vehtari et al., 2020). According to Gelman et al. (2014) a R-hat value of 1.1 or smaller for all parameters can be considered evidence for convergence. The Stan Development Team (2024) recommends running at least four chains and a convergence criterion of less than 1.05 for the maximum of rank normalized split-R-hat and rank normalized folded-split-R-hat. Vehtari et al. (2020), however, recommended to only use the posterior samples if R-hat is less than 1.01 because the R-hat can fall below 1.1 well before convergence in some scenarios (Brooks & Gelman, 1998; Vats & Knudon, 2018).

• Effective Sample Size (ESS): The ESS is the estimated number of independent samples from the posterior distribution that would lead to the same precision as the autocorrelated samples at hand. According to the default setting, the function computes the ESS based on rank-normalized split-R-hat and within-chain autocorrelation. The function provides the estimated Bulk-ESS (B.ESS) and the Tail-ESS (T.ESS). The Bulk-ESS is a useful measure for sampling efficiency in the bulk of the distribution (i.e, efficiency of the posterior mean), and the Tail-ESS is useful measure for sampling efficiency in the tails of the distribution (e.g., efficiency of tail quantile estimates). Note that by default, the Tail-ESS is the minimum of the effective sample sizes for 2.5% and 97.5% quantiles (tail = c(0.025, 0.975)). According to Kruschke (2015), a rank-normalized ESS greater than 400 is usually sufficient to get a stable estimate of the Monte Carlo standard error. However, a ESS of at least 1000 is considered optimal (Zitzmann & Hecht, 2019).

• Monte Carlo Standard Error (MCSE): The MCSE is defined as the standard deviation of the chains divided by their effective sample size and reflects uncertainty due to the stochastic algorithm of the Markov Chain Monte Carlo method. The function provides the estimated Bulk-MCSE (B.MCSE) for the margin of error when using the MCMC samples to estimate the posterior mean and the Tail-ESS (T.MCSE) for the margin of error when using the MCMC samples for interval estimation.

Value

Returns an object of class misty.object, which is a list with following entries:

Note

This function is a modified copy of functions provided in the rstan package by Stan Development Team (2024) and bayestestR package by Makowski et al. (2019).

Author(s)

Takuya Yanagida

References

Brooks, S. P. and Gelman, A. (1998). General Methods for Monitoring Convergence of Iterative Simulations. Journal of Computational and Graphical Statistics, 7(4): 434–455. MR1665662.

Gelman, A., & Rubin, D.B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457-472. https://doi.org/10.1214/ss/1177011136

Keller, B. T., & Enders, C. K. (2023). Blimp user’s guide (Version 3). Retrieved from www.appliedmissingdata.com/blimp

Kruschke, J. (2015). Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. Academic Press.

Makowski, D., Ben-Shachar, M., & Lüdecke, D. (2019). bayestestR: Describing effects and their uncertainty, existence and significance within the Bayesian framework. Journal of Open Source Software, 4(40), 1541. https://doi.org/10.21105/joss.01541

Stan Development Team (2024). RStan: the R interface to Stan. R package version 2.32.6. https://mc-stan.org/.

Vats, D. and Knudson, C. (2018). Revisiting the Gelman-Rubin Diagnostic. arXiv:1812.09384.

Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., & Bürkner, P.-C. (2020). Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC. Bayesian analysis, 16(2), 667-718. https://doi.org/110.1214/20-BA1221

Zitzmann, S., & Hecht, M. (2019). Going beyond convergence in Bayesian estimation: Why precision matters too and how to assess it. Structural Equation Modeling, 26(4), 646–661. https://doi.org/10.1080/10705511.2018.1545232

See Also

blimp, blimp.update, blimp.run, blimp.plot,blimp.print, blimp.plot,

Examples

## Not run: 

#————————————————————————————————————————————————————————————————————————————
# Blimp Example 4.3: Linear Regression

# Example 1a: Default setting, specifying name of the folder
blimp.bayes("Posterior_Ex4.3")

# Example 1b: Default setting, specifying the posterior file
blimp.bayes("Posterior_Ex4.3/posterior.csv")

# Example 2a: Print all summary measures, convergence, and efficiency diagnostics
blimp.bayes("Posterior_Ex4.3", print = "all")

# Example 3a: Print default measures plus MAP
blimp.bayes("Posterior_Ex4.3", print = c("default", "map"))

# Example 4: Print traditional R-hat in line with Blimp
blimp.bayes("Posterior_Ex4.3", split = TRUE, rank = FALSE, fold = FALSE)

# Example 5: Print probability of direction and the probability of
# being ROPE [-0.1, 0.1]
blimp.bayes("Posterior_Ex4.3", pd = TRUE, rope = c(-0.1, 0.1))

# Example 6: Write Results into a text file
blimp.bayes("Posterior_Ex4.3", write = "Bayes_Summary.txt")

# Example 7b: Write Results into an Excel file
blimp.bayes("Posterior_Ex4.3", write = "Bayes_Summary.xlsx")

## End(Not run)

Description

This function reads the posterior distribution including burn-in and post-burn-in phase for all parameters saved in long format in a file called posterior.* by the function blimp.run or blimp when specifying posterior = TRUE to display trace plots and posterior distribution plots.

Usage

blimp.plot(x, plot = c("none", "trace", "post"), param = NULL, labels = TRUE,
           burnin = TRUE, point = c("all", "none", "m", "med", "map"),
           ci = c("none", "eti", "hdi"), conf.level = 0.95, hist = TRUE,
           density = TRUE, area = TRUE, alpha = 0.4, fill = "gray85",
           facet.nrow = NULL, facet.ncol = NULL,
           facet.scales = c("fixed", "free", "free_x", "free_y"),
           xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL,
           xbreaks = ggplot2::waiver(), ybreaks = ggplot2::waiver(),
           xexpand = ggplot2::waiver(), yexpand = ggplot2::waiver(),
           palette = "Set 2", binwidth = NULL, bins = NULL,
           density.col = "#0072B2", shape = 21,
           point.col = c("#CC79A7", "#D55E00", "#009E73"),
           linewidth = 0.6, linetype = "dashed", line.col = "black",
           plot.margin = NULL, legend.title.size = 10, legend.text.size = 10,
           legend.box.margin = NULL, saveplot = c("all", "none", "trace", "post"),
           filename = "Blimp_Plot.pdf", file.plot = c("_TRACE", "_POST"),
           width = NA, height = NA, units = c("in", "cm", "mm", "px"), dpi = 600,
           check = TRUE)

Arguments

Value

Returns an object of class misty.object, which is a list with following entries:

Author(s)

Takuya Yanagida

References

Keller, B. T., & Enders, C. K. (2023). Blimp user’s guide (Version 3). Retrieved from www.appliedmissingdata.com/blimp

See Also

blimp, blimp.update, blimp.run, blimp.print, blimp.plot, blimp.bayes

Examples

## Not run: 

#————————————————————————————————————————————————————————————————————————————
# Blimp Example 4.3: Linear Regression

#···················
# Trace Plots

# Example 1a: Default setting, specifying name of the folder
blimp.plot("Posterior_Ex4.3")

# Example 1b: Default setting, specifying the posterior file
blimp.plot("Posterior_Ex4.3/posterior.csv")

# Example 1c: Print parameters 2, 3, 4, and 5
blimp.plot("Posterior_Ex4.3", param = 2:5)

# Example 1e: Arrange panels in three columns
blimp.plot("Posterior_Ex4.3", ncol = 3)

# Example 1f: Specify "Pastel 1" palette for the hcl.colors function
blimp.plot("Posterior_Ex4.3", palette = "Pastel 1")

#···················
# Posterior Distribution Plots

# Example 2a: Default setting, i.e., posterior median and equal-tailed interval
blimp.plot("Posterior_Ex4.3", plot = "post")

# Example 2b: Display posterior mean and maximum a posteriori
blimp.plot("Posterior_Ex4.3", plot = "post", point = c("m", "map"))

# Example 2c: Display maximum a posteriori and highest density interval
blimp.plot("Posterior_Ex4.3", plot = "post", point = "map", ci = "hdi")

# Example 2d: Do not display any point estimates and credible interval
blimp.plot("Posterior_Ex4.3", plot = "post", point = "none", ci = "none")

# Example 2d: Do not display histograms
blimp.plot("Posterior_Ex4.3", plot = "post", hist = FALSE)

#···················
# Save Plots

# Example 3a: Save all plots in pdf format
blimp.plot("Posterior_Ex4.3", saveplot = "all")

# Example 3b: Save all plots in png format with 300 dpi
blimp.plot("Posterior_Ex4.3", saveplot = "all", filename = "Blimp_Plot.png", dpi = 300)

# Example 3a: Save posterior distribution plot, specify width and height of the plot
blimp.plot("Posterior_Ex4.3", plot = "none", saveplot = "post",
           width = 7.5, height = 7)

#————————————————————————————————————————————————————————————————————————————
# Plot from misty.object

# Create misty.object
object <- blimp.plot("Posterior_Ex4.3", plot = "none")

# Trace plot
blimp.plot(object, plot = "trace")

# Posterior distribution plot
blimp.plot(object, plot = "post")

#————————————————————————————————————————————————————————————————————————————
# Create Plots Manually

# Load ggplot2 package
library(ggplot2)

# Create misty object
object <- blimp.plot("Posterior_Ex4.3", plot = "none")

#···················
# Example 4: Trace Plots

# Extract data
data.trace <- object$data$trace

# Plot
ggplot(data.trace, aes(x = iter, y = value, color = chain)) +
  annotate("rect", xmin = 0, xmax = 1000, ymin = -Inf, ymax = Inf,
           alpha = 0.4, fill = "gray85") +
  geom_line() +
  facet_wrap(~ param, ncol = 2, scales = "free") +
  scale_x_continuous(name = "", expand = c(0.02, 0)) +
  scale_y_continuous(name = "", expand = c(0.02, 0)) +
  scale_colour_manual(name = "Chain",
                      values = hcl.colors(n = 2, palette = "Set 2")) +
  theme_bw() +
  guides(color = guide_legend(nrow = 1, byrow = TRUE)) +
  theme(plot.margin = margin(c(4, 15, -10, 0)),
        legend.position = "bottom",
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.box.margin = margin(c(-16, 6, 6, 6)),
        legend.background = element_rect(fill = "transparent"))

#···················
# Example 5: Posterior Distribution Plots

# Extract data
data.post <- object$data$post

# Plot
ggplot(data.post, aes(x = value)) +
  geom_histogram(aes(y = after_stat(density)), color = "black", alpha = 0.4,
                 fill = "gray85") +
  geom_density(color = "#0072B2") +
  geom_vline(data = data.frame(param = levels(data.post$param),
                               stat = tapply(data.post$value, data.post$param, median)),
             aes(xintercept = stat, color = "Median"), linewidth = 0.6) +
  geom_vline(data = data.frame(param = levels(data.post$param),
                               low = tapply(data.post$value, data.post$param,
                                            function(y) quantile(y, probs = 0.025))),
             aes(xintercept = low), linetype = "dashed", linewidth = 0.6) +
  geom_vline(data = data.frame(param = levels(data.post$param),
                               upp = tapply(data.post$value, data.post$param,
                                            function(y) quantile(y, probs = 0.975))),
             aes(xintercept = upp), linetype = "dashed", linewidth = 0.6) +
  facet_wrap(~ param, ncol = 2, scales = "free") +
  scale_x_continuous(name = "", expand = c(0.02, 0)) +
  scale_y_continuous(name = "Probability Density, f(x)",
                     expand = expansion(mult = c(0L, 0.05))) +
  scale_color_manual(name = "Point Estimate", values = c(Median = "#D55E00")) +
  labs(caption = "95% Equal-Tailed Interval") +
  theme_bw() +
  theme(plot.margin = margin(4, 15, -8, 4),
        plot.caption = element_text(hjust = 0.5, vjust = 7),
        legend.position = "bottom",
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10),
        legend.box.margin = margin(-30, 6, 6, 6),
        legend.background = element_rect(fill = "transparent"))

## End(Not run)

Description

This function prints the result sections of a Blimp output file (.blimp-out) on the R console. By default, the function prints selected result sections, i.e., Algorithmic Options Specified, Data Information, Model Information, Warning Messages, Outcome Model Estimates, and Generated Parameters.

Usage

blimp.print(x,
            result = c("all", "default", "algo.options", "data.info",
            "model.info", "warn.mess", "error.mess", "out.model", "gen.param"),
            exclude = NULL, color = c("none", "blue", "green"),
            style = c("bold", "regular"), not.result = TRUE,
            write = NULL, append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

Result Sections

Following result sections can be selected by using the result argument or excluded by using the exclude argument:

• "algo.options" for the ALGORITHMIC OPTIONS SPECIFIED section

• "simdat.summary" for the SIMULATED DATA SUMMARIES section

• "order.simdat" for the VARIABLE ORDER IN SIMULATED DATA section

• "burnin.psr" for the BURN-IN POTENTIAL SCALE REDUCTION (PSR) OUTPUT section

• "mh.accept" for the METROPOLIS-HASTINGS ACCEPTANCE RATES section

• "data.info" for the DATA INFORMATION section

• "var.imp" for the VARIABLES IN IMPUTATION MODEL section

• "model.info" for the MODEL INFORMATION section

• "param.label" for the PARAMETER LABELS section

• "warn.mess" for the WARNING MESSAGES section

• "fit" for the MODEL FIT section

• "cor.resid" for the CORRELATIONS AMONG RESIDUALS section

• "out.model" for the OUTCOME MODEL ESTIMATES section

• "pred.model" for the PREDICTOR MODEL ESTIMATES section

• "gen.param" for the GENERATED PARAMETERS section

• "order.impdat" for the VARIABLE ORDER IN IMPUTED DATA section

Note that all result sections are requested by specifying result = "all". The result argument is also used to select one (e.g., result = "algo.options") or more than one result sections (e.g., result = c("algo.options", "fit")), or to request result sections in addition to the default setting (e.g., result = c("default", "fit")). The exclude argument is used to exclude result sections from the output (e.g., exclude = "algo.options").

Value

Returns an object of class misty.object, which is a list with following entries:

Author(s)

Takuya Yanagida

References

Keller, B. T., & Enders, C. K. (2023). Blimp user’s guide (Version 3). Retrieved from www.appliedmissingdata.com/blimp

See Also

blimp, blimp.update, blimp.run, blimp.plot, blimp.bayes

Examples

## Not run: 

#————————————————————————————————————————————————————————————————————————————
# Blimp Example 4.3: Linear Regression

# Example 1a: Default setting
blimp.print("Ex4.3.blimp-out")

# Example 1c: Print OUTCOME MODEL ESTIMATES only
blimp.print("Ex4.3.blimp-out", result = "out.model")

# Example 1d: Print MODEL FIT in addition to the default setting
blimp.print("Ex4.3.blimp-out", result = c("default", "fit"))

# Example 1e: Exclude DATA INFORMATION section
blimp.print("Ex4.3.blimp-out", exclude = "data.info")

# Example 1f: Print all result sections, but exclude MODEL FIT section
blimp.print("Ex4.3.blimp-out", result = "all", exclude = "fit")

# Example 1g: Print result section in a different order
blimp.print("Ex4.3.blimp-out", result = c("model.info", "fit", "algo.options"))

#————————————————————————————————————————————————————————————————————————————
# misty.object of type 'blimp.print'

# Example 2
# Create misty.object
object <- blimp.print("Ex4.3.blimp-out", output = FALSE)

# Print misty.object
blimp.print(object)

#————————————————————————————————————————————————————————————————————————————
# Write Results

# Example 3: Write Results into a text file
blimp.print("Ex4.3.blimp-out", write = "Output_4-3.txt")

## End(Not run)

Description

This function runs a group of Blimp models (.imp files) located within a single directory or nested within subdirectories.

Usage

blimp.run(target = getwd(), recursive = FALSE,
          replace.out = c("always", "never", "modified"), posterior = FALSE,
          folder = "Posterior_", format = c("csv", "csv2", "xlsx", "rds", "RData"),
          clear = FALSE, Blimp = .detect.blimp(), check = TRUE)

Arguments

Value

None.

Note

This function is based on the detect_blimp() and rblimp() function in the rblimp package by Brian T.Keller (2024).

Author(s)

Takuya Yanagida

References

Keller, B. T., & Enders, C. K. (2023). Blimp user’s guide (Version 3). Retrieved from www.appliedmissingdata.com/blimp

Keller B (2024). rblimp: Integration of Blimp Software into R. R package version 0.1.31. https://github.com/blimp-stats/rblimp

See Also

blimp, blimp.update, blimp.print, blimp.plot, blimp.bayes

Examples

## Not run: 

# Example 1: Run Blimp models located within the current working directory
blimp.run()

# Example 2: Run Blimp models located nested within subdirectories
blimp.run(recursive = TRUE)

# Example 3: Run Blimp input file
blimp.run("Ex4.1a.imp")

# Example 4: Run Blimp input files
blimp.run(c("Ex4.1a.imp", "Ex4.1b.imp"))

# Example 5: Run Blimp models, save posterior distribution in a R workspace
blimp.run(posterior = TRUE, format = "workspace")

## End(Not run)

Description

This function updates specific input command sections of a misty.object of type blimp to create an updated Blimp input file, run the updated input file by using the blimp.run() function, and print the updated Blimp output file by using the blimp.print() function.

Usage

blimp.update(x, update, file = "Blimp_Input_Update.imp", comment = FALSE,
             replace.inp = TRUE, blimp.run = TRUE, posterior = FALSE,
             folder = "Posterior_",
             format = c("csv", "csv2", "xlsx", "rds", "RData"),
             clear = TRUE, replace.out = c("always", "never", "modified"),
             Blimp = .detect.blimp(),
             result = c("all", "default", "algo.options",  "data.info",
                        "model.info", "warn.mess", "out.model", "gen.param"),
             exclude = NULL, color = c("none", "blue", "violet"),
             style = c("bold", "regular"), not.result = TRUE,
             write = NULL, append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

Bimp Input Sections

The function is used to update following Blimp input sections:

• DATA

• VARIABLES

• CLUSTERID

• ORDINAL

• NOMINAL

• COUNT

• WEIGHT

• MISSING

• LATENT

• RANDOMEFFECT

• TRANSFORM

• BYGROUP

• FIXED

• CENTER

• MODEL

• SIMPLE

• PARAMETERS

• TEST

• FCS

• SIMUALTE

• SEED

• BURN

• ITERATIONS

• CHAINS

• NIMPS

• THIN

• OPTIONS

• OUTPUT

• SAVE

The ---; Specification

The ---; specification is used to remove entire sections (e.g., CENTER: ---;) from the Blimp input. Note that ---; including the semicolon ; needs to be specified, i.e., --- without the semicolon ; will result in an error message.

Value

Returns an object of class misty.object, which is a list with following entries:

Author(s)

Takuya Yanagida

References

Keller, B. T., & Enders, C. K. (2023). Blimp user’s guide (Version 3). Retrieved from www.appliedmissingdata.com/blimp

See Also

blimp.run, blimp.print, blimp.plot, blimp.bayes

Examples

## Not run: 

#————————————————————————————————————————————————————————————————————————————
# Example 1a: Update BURN and ITERATIONS section

# Specify Blimp input
input <- '
DATA: data1.csv;
ORDINAL: d;
MISSING: 999;
FIXED: d;
CENTER: x1 x2;
MODEL: y ~ x1 x2 d;
SEED: 90291;
BURN: 1000;
ITERATIONS: 10000;
'

# Run Blimp input
mod0 <- blimp(input, file = "Ex4.3.imp", clear = FALSE)

# Update sections
update1 <- '
BURN: 5000;
ITERATIONS: 20000;
'

# Run updated Blimp input
mod1 <- blimp.update(mod0, update1, file = "Ex4.3_update1.imp")

#————————————————————————————————————————————————————————————————————————————
# Example 1b: Remove CENTER section

# Remove section
update2 <- '
CENTER: ---;
'

# Run updated Blimp input
mod2 <- blimp.update(mod1, update2, file = "Ex4.3_update2.imp")

## End(Not run)

Description

This function performs the model-based Bollen-Stine Bootstrapping with incomplete data of the chi-square statistic. By default, the function performs model-based bootstrapping based on transformation method 2 in Savalei and Yuan (2009).

Usage

boot.bs(object = NULL, data = NULL, model = NULL, sigma = NULL, mu = NULL,
        group = NULL, chisq = NULL, em.cov = NULL, trans = c(1, 2), R = 500,
        return = c("transdat", "bootsamp", "output"), seed = NULL,
        progress = TRUE, digits = 2, p.digits = 3, plot = FALSE, filename = NULL,
        width = NA, height = NA, dpi = 600, write = NULL, append = TRUE,
        check = TRUE, output = TRUE, ...)

Arguments

Value

Returns an object of class misty.object when specifying return = "output":

When specifying return = "transdat", the transformed data and when specifying return = "bootsamp", the bootstrap samples are returned.

Note

This function is based on modified copies of the functions bsBootMiss from the semTools package by Terrence D. Jorgensen et al. (2026).

Author(s)

Takuya Yanagida

References

Bollen, K. A., & Stine, R. A. (1992). Bootstrapping goodness-of-fit measures in structural equation models. Sociological Methods & Research, 21(2), 205-229. https://doi.org/10.1177/0049124192021002004

Jorgensen, T. D., Pornprasertmanit, S., Schoemann, A. M., & Rosseel, Y. (2026). semTools: Useful tools for structural equation modeling. R package version 0.5-8. Retrieved from https://CRAN.R-project.org/package=semTools

Savalei, V., & Yuan, K.-H. (2009). On the model-based bootstrap with missing data: Obtaining a p-value for a test of exact fit. Multivariate Behavioral Research, 44(6), 741-763. https://doi.org/10.1080/00273170903333590

Examples

## Not run: 
# Load lavaan package
library(lavaan)

# Holzinger and Swineford data set
dat <- HolzingerSwineford1939

# Introduce missing data
dat$x5 <- ifelse(dat$x1 <= quantile(dat$x1, 0.3), NA, dat$x5)
dat$x9 <- ifelse(is.na(dat$x5), NA, dat$x9)

# Model specification
model <- 'visual  =~ x1 + x2 + x3
          textual =~ x4 + x5 + x6
          speed   =~ x7 + x8 + x9'

# Model estimation
fit <- sem(model, data = dat, meanstructure = TRUE, std.lv = TRUE,
           missing = "fiml", group = "school")

#————————————————————————————————————————————————————————————————————————————
# Bollen-Stine Bootstrapping with Incomplete Data

# Example 1: Default setting, transformation method 2, R = 500 replicates
# Plot bootstrap sampling distribution of the test statistic
boot.bs(fit, seed = 42, plot = TRUE)

#————————————————————————————————————————————————————————————————————————————
# Transformed Data and Bootstrap Samples

# Example 2: Return transformed data only
transdat <- boot.bs(fit, return = "transdat")

# Example 3: Return bootstrap samples only
bootsamp <- boot.bs(fit, return = "bootsamp")

#————————————————————————————————————————————————————————————————————————————
# Plot Bootstrap Sampling Distribution of Chi-Square Test Statistic

# Bollen-Stine Bootstrapping
object <- boot.bs(fit, seed = 42)

# Load ggplot2 package
library(ggplot2)

# Plot data
plotdat <- data.frame(chisq = object$boot.chisq)

# Example 3: Plot bootstrap sampling distribution, create plot manually
ggplot(plotdat, aes(chisq)) +
  geom_histogram(aes(y = after_stat(density)), color = "black", alpha = 0.4, fill = "gray85") +
  geom_density(color = "#0072B2") +
  geom_vline(aes(xintercept = object$result$chisq, color = "Observed Test Statistic")) +
  scale_x_continuous(name = expression(paste(chi^2, " Test Statistic")),
                     limits = c(0, max(c(plotdat$chisq, object$result$chisq), na.rm = TRUE))) +
  scale_y_continuous(name = "Probability Density, f(x)", expand = expansion(mult = c(0, 0.05))) +
  scale_color_manual(values = c("Observed Test Statistic" = "#CC79A7")) +
  theme_bw() +
  theme(legend.position = "bottom", legend.box.margin = margin(-12, 0, 0, 0),
        legend.title = element_blank())

## End(Not run)

Description

This function centers predictor variables in single-level data, two-level data, and three-level data at the grand mean (CGM, i.e., grand mean centering) or within clusters (CWC, i.e., group mean centering).

Usage

center(data, ..., cluster = NULL, type = c("CGM", "CWC", "latent"),
       cwc.mean = c("L2", "L3"), value = NULL, append = TRUE, name = ".c",
       as.na = NULL, check = TRUE)

Arguments

Details

Single-Level Data

Predictor variables are centered at the grand mean (CGM) by default:

$x_{i} - \bar{x}_{.}$

where $x_{i}$ is the predictor value of observation $i$ and $\bar{x}_{.}$ is the average $x$ score. Note that predictor variables can be centered on any meaningful value specifying the argument value, e.g., a predictor variable centered at 5 by applying following formula:

$x_{i} - \bar{x}_{.} + 5$

resulting in a mean of the centered predictor variable of 5.

Two-Level Data

In two-level data, there are predictor variables at Level-1 (L1) and Level-2 (L2) with L1 predictor variables centered within L2 clusters (CWC) and L2 predictors centered at the average L2 cluster scores (CGM) by default:

• Level-1 (L1) Predictor Variables:

  L1 predictor variable can be centered within L2 clusters (CWC) or at the grand-mean (CGM):

  • L1 predictor variables are centered within L2 clusters by specifying type = "CWC" (Default):

    $x_{ij} - \bar{x}_{.j}$

    where $\bar{x_{.j}}$ is the average $x$ score in cluster $j$.

  • L1 predictor variables are centered at the grand-mean by specifying type = "CGM":

    $x_{ij} - \bar{x}_{..}$

    where $x_{ij}$ is the predictor value of observation $i$ in L2 cluster $j$ and $\bar{x}_{..}$ is the average $x$ score.

• Level-2 (L2) Predictor Variables:

  L2 predictor variables are centered at the average L2 cluster score:

  $x_{.j} - \bar{x}_{..}$

  where $x_{.j}$ is the predictor value of L2 cluster $j$ and $\bar{x}_{..}$ is the average L2 cluster score. Note that the cluster membership variable needs to be specified when centering a L2 predictor variable in two-level data. Otherwise the average $x_{ij}$ individual score instead of the average $x_{.j}$ cluster score is used to center the predictor variable.

Three-Level Data

In three-level data, there are predictor variables at Level-1 (L1), Level-2 (L2), and Level-3 (L3) with L1 predictor variables centered within L2 clusters (CWC L2), L2 predictors centered within L3 clusters (CWC L3), and L3 predictors centered at the average L3 cluster scores (CGM) by default:

• Level-1 (L1) Predictor Variables:

  L1 predictor variables can be centered within L2 clusters (CWC L2), within L3 clusters (CWC L3) or at the grand-mean (CGM):

  • L1-predictor variables are centered within cluster (CWC) by specifying type = "CWC" (Default). Note that L1 predictor variables can be either centered within L2 clusters (cwc.mean = "L2", Default, see Brincks et al., 2017):

    $x_{ijk} - \bar{x}_{.jk}$

    or within L3 clusters (cwc.mean = "L3", see Enders, 2013):

    $x_{ijk} - \bar{x}_{..k}$

    where $\bar{x}_{.jk}$ is the average $x$ score in L2 cluster $j$ within Level-3 cluster $k$ and $\bar{x}_{..k}$ is the average $x$ score in L3 cluster $k$.

  • L1 predictor variables are centered at the grand mean (CGM) by specifying type = "CGM":

    $x_{ijk} - \bar{x}_{...}$

    where $x_{ijk}$ is the predictor value of observation $i$ in L2 cluster $j$ within L3 cluster $k$ and $\bar{x}_{...}$ is the average $x$ score.

• Level-2 (L2) Predictor Variables:

  L2 predictor variables can be centered within L3 clusters (CWC) or at the L2 grand-mean (CGM):

  • L2 predictor variables are centered within cluster by specifying type = "CWC" (Default):

    $x_{.jk} - \bar{x}_{..k}$

    where $\bar{x}_{..k}$ is the average $x$ score in L3 cluster $k$.

  • L2 predictor variables are centered at the grand mean by specifying type = "CGM":

    $x_{.jk} - \bar{x}_{...}$

    where $x_{.jk}$ is the predictor value of L2 cluster $j$ within L3 cluster $k$ and $\bar{x}_{...}$ is the average L2 cluster score.

• Level-3 (L3) Predictor Variables:

  L3-predictor variables are centered at the L3 grand mean:

  $x_{..k} - \bar{x}_{...}$

  where $x_{..k}$ is the predictor value of L3 cluster $k$ and $\bar{x}_{...}$ is the average L3 cluster score.

Two-Step Latent Mean Centering

The latent mean centering approach (Asparouhov & Muthén, 2019) in a two-level model decomposes the Level-1 predictor variable $x_{ij}$ as within and between compoments as follows:

$x_{ij} = x_{w,ij} + x_{b,.j}$

where $x_{w,ij}$ is the individual specific contribution and $x_{b,.j}$ is the cluster specific contribution to the predictor variable $x_{ij}$. Here, $x_{b,.j}$ can be interpreted as the intercepts and $x_{w,ij}$ can be interpreted as the residuals in the random intercept model. Note that $x_{w,ij}$ is equivalent to a L1 predictor centered within L2 clusters (CWC), while $x_{b,.j}$ is equivalent to a L2 predictor centered at the average L2 cluster scores (CGM). Latent mean centering treats $x_{b,.j}$ as unknown quantity that is estimated while taking into the sampling error in the mean estimate under the assumption of large cluster sizes in the population and less than 5% of the cluster population sampled. As a result, this approach resolves problems that occur with the traditional observed centering methods, e.g., Lüdtke's bias (Lüdtke et al., 2008) in the estimation of contextual effects or Nickell's bias (Asparouhov et al., 2018) in the estimation of the autocorrelations in time-series models.

The latent mean centering approach requires a latent variable modeling program, e.g., commercial software Mplus (Muthen & Muthen, 1998-2017) or the R package lavaan (Rosseel, 2012) and cannot be used in mixed-effects modeling programs like lme4 (Bates et al., 2015) or nlme (Pinheiro & Bates, 2000). In order to mimic the latent mean centering approach, a two-step approach is proposed in the center() function, where a random intercept model is fit to the L1 predictor variable to extract the intercepts representing $x_{b,.j}$ and residuals representing $x_{w,ij}$. These two components can be used as L1 predictor centered within clusters and L2 predictor centered at the grand mean. Note that compared to the latent mean centering approach, this two-step approach will result in bias because $x_{w,ij}$ and $x_{b,.j}$ are treated as observed instead of latent variables. However, the magnitude of the bias is unclear without conducting a simulation study. Hence, the latent mean centering using a latent variable modeling program is recommended whenever possible, while the two-step latent mean centering approach implemented in the center() function is just an 'experimental' approach that cannot be recommend at this time.

Value

Returns a numeric vector or data frame with the same length or same number of rows as data containing the centered variable(s).

Author(s)

Takuya Yanagida [email protected]

References

Asparouhov, T., Hamaker, E. L., & Muthén, B. (2017). Dynamic Structural Equation Models. Structural Equation Modeling: A Multidisciplinary Journal, 25(3), 359-388. https://doi.org/10.1080/10705511.2017.1406803

Asparouhov, T., & Muthén, B. (2019). Latent variable centering of predictors and mediators in multilevel and time-series models. Structural Equation Modeling, 26(1), 119-142. https://doi.org/10.1080/10705511.2018.1511375

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01

Brincks, A. M., Enders, C. K., Llabre, M. M., Bulotsky-Shearer, R. J., Prado, G., & Feaster, D. J. (2017). Centering predictor variables in three-level contextual models. Multivariate Behavioral Research, 52(2), 149–163. https://doi.org/10.1080/00273171.2016.1256753

Chang, C.-N., & Kwok, O.-M. (2022) Partitioning Variance for a Within-Level Predictor in Multilevel Models. Structural Equation Modeling: A Multidisciplinary Journal. Advance online publication. https://doi.org/10.1080/10705511.2022.2051175

Enders, C. K. (2013). Centering predictors and contextual effects. In M. A. Scott, J. S. Simonoff, & B. D. Marx (Eds.), The Sage handbook of multilevel modeling (pp. 89-109). Sage. https://dx.doi.org/10.4135/9781446247600

Enders, C. K., & Tofighi, D. (2007). Centering predictor variables in cross-sectional multilevel models: A new look at an old issue. Psychological Methods, 12, 121-138. https://doi.org/10.1037/1082-989X.12.2.121

Lüdtke, O., Marsh, H. W., Robitzsch, A., Trautwein, U., Asparouhov, T., & Muthén, B. (2008). The multilevel latent covariate model: A new, more reliable approach to group-level effects in contextual studies. Psychological Methods, 13(3), 203-229. https://doi.org/10.1037/a0012869

Muthén, L. K., & Muthén, B. O. (1998-2017). Mplus User’s Guide (8th ed). Muthén & Muthén.

Pinheiro, J. C., & Bates, D. M. (2000). Mixed-Effects Models in S and S-PLUS. Springer. https://doi.org/10.1007/b98882

Rights, J. D., Preacher, K. J., & Cole, D. A. (2020). The danger of conflating level-specific effects of control variables when primary interest lies in level-2 effects. British Journal of Mathematical & Statistical Psychology, 73, 194-211. https://doi.org/10.1111/bmsp.12194

Rosseel, Y. (2012). lavaan: An R Package for Structural Equation Modeling. Journal of Statistical Software, 48(2), 1-36. https://doi.org/10.18637/jss.v048.i02

Yaremych, H. E., Preacher, K. J., & Hedeker, D. (2021). Centering categorical predictors in multilevel models: Best practices and interpretation. Psychological Methods. Advance online publication. https://doi.org/10.1037/met0000434

See Also

coding, cluster.scores, rec, item.reverse, cluster.rwg, item.scores.

Examples

#————————————————————————————————————————————————————————————————————————————
# Single-Level Data

# Example 1a: Center predictor 'disp' at the grand mean
center(mtcars, disp, append = FALSE)

# Alternative specification without using the '...' argument
center(mtcars$disp)

# Example 1b: Center predictors 'disp' and 'hp' at the grand mean and append to 'mtcars'
center(mtcars, disp, hp)

# Alternative specification without using the '...' argument
cbind(mtcars, center(mtcars[, c("disp", "hp")]))

# Example 1c: Center predictor 'disp' at the value 3
center(mtcars, disp, value = 3)

# Example 1d: Center predictors 'disp' and 'hp' and label with the suffix ".v"
center(mtcars, disp, hp, name = ".v")

#————————————————————————————————————————————————————————————————————————————
# Two-Level Data

# Load data set "Demo.twolevel" in the lavaan package
data("Demo.twolevel", package = "lavaan")

#·········································
# Level-1 (L1) Predictor

# Example 2a: Center L1 predictor 'y1' within L2 clusters
center(Demo.twolevel, y1, cluster = "cluster", append = FALSE)

# Alternative specification without using the '...' argument
center(Demo.twolevel$y1, cluster = Demo.twolevel$cluster)

# Example 2b: Center L1 predictor 'y1' at the grand-mean
#             Note that cluster ID is ignored when type = "CGM"
center(Demo.twolevel, y1, cluster = "cluster", type = "CGM")

# Alternative specification
center(Demo.twolevel, y1)

#·········································
# Level-2 (L2) Predictor

# Example 2c: Center L2 predictor 'w2' at the average L2 cluster scores
#             Note that cluster ID is needed
center(Demo.twolevel, w1, cluster = "cluster")

#·········································
# L1 and L2 Predictors

# Example 2d: Center L1 predictor 'y1' within L2 clusters
#             and L2 predictor 'w1' at the average L2 cluster scores
center(Demo.twolevel, y1, w1, cluster = "cluster")

#·········································
# Two-Step Latent Mean Centering

# Example 2e: Decompose L1 predictor 'y1' as within-between components
center(Demo.twolevel, y1, cluster = "cluster", type = "latent")

# Example 2d: Decompose L1 predictor 'y1' as within-between components
#             label variables as 'l1.y1' and 'l2.y1'
center(Demo.twolevel, y1, cluster = "cluster", type = "latent", name = c("l1.y1", "l2.y1"))

## Not run: 
#————————————————————————————————————————————————————————————————————————————
# Three-Level Data

# Load data set "Demo.twolevel" in the lavaan package
data("Demo.twolevel", package = "lavaan")

# Create arbitrary three-level data
Demo.threelevel <- data.frame(Demo.twolevel, cluster2 = Demo.twolevel$cluster,
                                             cluster3 = rep(1:10, each = 250))

# Compute L3 cluster scores for the L2 predictor 'w1'
Demo.threelevel <- cluster.scores(Demo.threelevel, w1, cluster = "cluster3", name = "w1.l3")

#·········································
# Level-1 (L1) Predictor

# Example 3a: Center L1 predictor 'y1' within L2 clusters (CWC L2)
#             Note that L3 cluster IDs are ignored when type = "CWC"
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"))

# Alternative specification when L2 cluster IDs are unique across L3 clusters
center(Demo.threelevel, y1, cluster = "cluster2")

# Example 3b: Center L1 predictor 'y1' within L3 clusters (CWC L3)
#             Note that both L3 and L2 cluster IDs are needed
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"), cwc.mean = "L3")

# Example 3c: Center L1 predictor 'y1' at the grand-mean (CGM)
#             Note that the cluster argument is ignored when type = "CGM",
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"), type = "CGM")

# Alternative specification
center(Demo.threelevel, y1)

#·········································
# Level-2 (L2) Predictor

# Example 3d: Center L2 predictor 'w1' within L3 cluster
#             Note that both L3 and L2 cluster IDs are needed
center(Demo.threelevel, w1, cluster = c("cluster3", "cluster2"))

# Example 3e: Center L2 predictor 'w1' at the grand-mean (CGM)
#             Note that both L3 and L2 cluster IDs are needed
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"), type = "CGM")

#·········································
# Level-3 (L3) Predictor

# Example 3f: Center L3 predictor 'w1.l3' at the average L3 cluster scores
#             Note that L2 cluster ID is ignored
center(Demo.threelevel, w1.l3, cluster = c("cluster3", "cluster2"))

# Alternative specification
center(Demo.threelevel, w1.l3, cluster = "cluster3")

#·········································
# L1, L2, and L3 Predictors

# Example 3g: Center L1 predictor 'y1' within L2 cluster, L2 predictor 'w1' within
#            L3 clusters, and L3 predictor 'w1.l3' at the average L3 cluster scores
center(Demo.threelevel, y1, w1, w1.l3, cluster = c("cluster3", "cluster2"))

#·········································
# Two-Step Latent Mean Centering

# Load data set "Demo.twolevel" in the lavaan package
data("Demo.twolevel", package = "lavaan")

# Create arbitrary three-level data
Demo.threelevel <- data.frame(Demo.twolevel, cluster2 = Demo.twolevel$cluster,
                                             cluster3 = rep(1:10, each = 250))

# Example 3h: Decompose L1 predictor 'y1' as within-between components
center(Demo.threelevel, y1, cluster = "cluster2", type = "latent")

# Example 3i: Decompose L1 predictor 'y1' as within-between components
#             label variables as 'l1.y1' and 'l2.y1'
center(Demo.threelevel, y1, cluster = "cluster2", type = "latent",
       name = c("l1.y1", "l2.y2"))

# Example 3j: Decompose L1 predictor 'y1' as within-between components
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"), type = "latent")

# Example 3k: Decompose L1 predictor 'y1' as within-between components
#             label variables as 'l1.y1', 'l2.y1', and 'l3.y1'
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"), type = "latent",
       name = c("l1.y1", "l2.y1", "l3.y1"))

## End(Not run)

Description

This function computes tolerance, standard error inflation factor, variance inflation factor, eigenvalues, condition index, and variance proportions for linear, generalized linear, and mixed-effects models.

Usage

check.collin(model, print = c("all", "vif", "eigen"), digits = 3, p.digits = 3,
             write = NULL, append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

Collinearity diagnostics can be conducted for objects returned from the lm() and glm() function, but also from objects returned from the lmer() and glmer() function from the lme4 package, lme() function from the nlme package, and the glmmTMB() function from the glmmTMB package.

The generalized variance inflation factor (Fox & Monette, 1992) is computed for terms with more than 1 df resulting from factors with more than two levels. The generalized VIF (GVIF) is interpretable as the inflation in size of the confidence ellipse or ellipsoid for the coefficients of the term in comparison with what would be obtained for orthogonal data. GVIF is invariant to the coding of the terms in the model. In order to adjust for the dimension of the confidence ellipsoid, GVIF$^\frac{1}{2df}$ is computed. Note that the adjusted GVIF (aGVIF) is actually a generalized standard error inflation factor (GSIF). Thus, the aGIF needs to be squared before applying a common cutoff threshold for the VIF (e.g., VIF > 10). Note that the output of check.collin() function reports either the variance inflation factor or the squared generalized variance inflation factor in the column VIF, while the standard error inflation factor or the adjusted generalized variance inflation factor is reported in the column SIF.

Value

Returns an object of class misty.object, which is a list with following entries:

Note

The computation of the VIF and the GVIF is based on the vif() function in the car package by John Fox, Sanford Weisberg and Brad Price (2020), and the computation of eigenvalues, condition index, and variance proportions is based on the ols_eigen_cindex() function in the olsrr package by Aravind Hebbali (2020).

Author(s)

Takuya Yanagida [email protected]

References

Fox, J., & Monette, G. (1992). Generalized collinearity diagnostics. Journal of the American Statistical Association, 87, 178-183.

Fox, J., Weisberg, S., & Price, B. (2020). car: Companion to Applied Regression. R package version 3.0-8. https://cran.r-project.org/web/packages/car/

Hebbali, A. (2020). olsrr: Tools for building OLS regression models. R package version 0.5.3. https://cran.r-project.org/web/packages/olsrr/

See Also

check.outlier, lm

Examples

dat <- data.frame(group = c(1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4),
                  x1 = c(3, 2, 4, 9, 5, 3, 6, 4, 5, 6, 3, 5),
                  x2 = c(1, 4, 3, 1, 2, 4, 3, 5, 1, 7, 8, 7),
                  x3 = c(7, 3, 4, 2, 5, 6, 4, 2, 3, 5, 2, 8),
                  x4 = c("a", "b", "a", "c", "c", "c", "a", "b", "b", "c", "a", "c"),
                  y1 = c(2, 7, 4, 4, 7, 8, 4, 2, 5, 1, 3, 8),
                  y2 = c(0, 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1), stringsAsFactors = TRUE)

#————————————————————————————————————————————————————————————————————————————
# Linear model

# Estimate linear model with continuous predictors
mod.lm1 <- lm(y1 ~ x1 + x2 + x3, data = dat)

# Example 1: Tolerance, std. error, and variance inflation factor
check.collin(mod.lm1)

# Example 2: Tolerance, std. error, and variance inflation factor
# Eigenvalue, Condition index, and variance proportions
check.collin(mod.lm1, print = "all")

# Estimate model with continuous and categorical predictors
mod.lm2 <- lm(y1 ~ x1 + x2 + x3 + x4, data = dat)

# Example 3: Tolerance, generalized std. error, and variance inflation factor
check.collin(mod.lm2)

#————————————————————————————————————————————————————————————————————————————
# Generalized linear model

# Estimate logistic regression model with continuous predictors
mod.glm <- glm(y2 ~ x1 + x2 + x3, data = dat, family = "binomial")

# Example 4: Tolerance, std. error, and variance inflation factor
check.collin(mod.glm)

## Not run: 
#————————————————————————————————————————————————————————————————————————————
# Linear mixed-effects model

# Load lme4, nlme, and glmmTMB package
libraries(lme4, nlme, glmmTMB)

# Estimate linear mixed-effects model using lme4 package
mod.lmer <- lmer(y1 ~ x1 + x2 + x3 + (1|group), data = dat)

# Example 5: Tolerance, std. error, and variance inflation factor
check.collin(mod.lmer)

# Estimate linear mixed-effects model using nlme package
mod.lme <- lme(y1 ~ x1 + x2 + x3, random = ~ 1 | group, data = dat)

# Example 6: Tolerance, std. error, and variance inflation factor
check.collin(mod.lme)

# Estimate linear mixed-effects model using glmmTMB package
mod.glmmTMB1 <- glmmTMB(y1 ~ x1 + x2 + x3 + (1|group), data = dat)

# Example 7: Tolerance, std. error, and variance inflation factor
check.collin(mod.glmmTMB1)

#————————————————————————————————————————————————————————————————————————————
# Generalized linear mixed-effects model

# Estimate mixed-effects logistic regression model using lme4 package
mod.glmer <- glmer(y2 ~ x1 + x2 + x3 + (1|group), data = dat, family = "binomial")

# Example 8: Tolerance, std. error, and variance inflation factor
check.collin(mod.glmer)

# Estimate mixed-effects logistic regression model using glmmTMB package
mod.glmmTMB2 <- glmmTMB(y2 ~ x1 + x2 + x3 + (1|group), data = dat, family = "binomial")

# Example 9: Tolerance, std. error, and variance inflation factor
check.collin(mod.glmmTMB2)

#————————————————————————————————————————————————————————————————————————————
# Write Results

# Example 10: Write Results into a text file
check.collin(mod.lm1, write = "Diagnostics.txt")

## End(Not run)