Package 'misty'

Description

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

Usage

aov.b(formula, data, posthoc = FALSE, conf.level = 0.95, hypo = TRUE,
      descript = TRUE, effsize = FALSE, weighted = FALSE, correct = FALSE,
      plot = FALSE, point.size = 4, adjust = TRUE, error.width = 0.1,
      xlab = NULL, ylab = NULL, ylim = NULL, breaks = ggplot2::waiver(),
      jitter = TRUE, jitter.size = 1.25, jitter.width = 0.05,
      jitter.height = 0, jitter.alpha = 0.1, title = "",
      subtitle = "Confidence Interval", digits = 2, p.digits = 4,
      as.na = NULL, write = NULL, append = TRUE, check = TRUE,
      output = TRUE, ...)

Arguments

Details

Post Hoc Test

Tukey HSD post hoc test reports Cohen's d based on the non-weighted standard deviation (i.e., weighted = FALSE) when requesting an effect size measure (i.e., effsize = TRUE) following the recommendation by Delacre et al. (2021).

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.

Delacre, M., Lakens, D., Ley, C., Liu, L., & Leys, C. (2021). Why Hedges' g*s based on the non-pooled standard deviation should be reported with Welch's t-test. https://doi.org/10.31234/osf.io/tu6mp

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, test.welch, cohens.d, ci.mean.diff, ci.mean

Examples

dat <- data.frame(group = c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3),
                  y = c(3, 1, 4, 2, 5, 3, 2, 3, 6, 6, 3, NA))

# Example 1: Between-subject ANOVA
aov.b(y ~ group, data = dat)

# Example 2: Between-subject ANOVA
# print effect size measures
aov.b(y ~ group, data = dat, effsize = TRUE)

# Example 3: Between-subject ANOVA
# do not print hypotheses and descriptive statistics,
aov.b(y ~ group, data = dat, descript = FALSE, hypo = FALSE)

## Not run: 
# Example 4: Write results into a text file
aov.b(y ~ group, data = dat, write = "ANOVA.txt")

# Example 5: Between-subject ANOVA
# plot results
aov.b(y ~ group, data = dat, plot = TRUE)

# Load ggplot2 package
library(ggplot2)

# Example 6: Save plot, ggsave() from the ggplot2 package
ggsave("Between-Subject_ANOVA.png", dpi = 600, width = 4.5, height = 6)

# Example 7: Between-subject ANOVA
# extract plot
p <- aov.b(y ~ group, data = dat, output = FALSE)$plot
p

# Extract data
plotdat <- aov.b(y ~ group, data = dat, output = FALSE)$data

# Draw plot in line with the default setting of aov.b()
ggplot(plotdat, aes(group, y)) +
  geom_jitter(alpha = 0.1, width = 0.05, height = 0, size = 1.25) +
  geom_point(stat = "summary", fun = "mean", size = 4) +
  stat_summary(fun.data = "mean_cl_normal", geom = "errorbar", width = 0.20) +
  scale_x_discrete(name = NULL) +
  labs(subtitle = "Two-Sided 95
  theme_bw() + theme(plot.subtitle = element_text(hjust = 0.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", "LB", "GG", "HF"),
      posthoc = FALSE, conf.level = 0.95,
      p.adj = c("none", "bonferroni", "holm", "hochberg", "hommel", "BH", "BY", "fdr"),
      hypo = TRUE, descript = TRUE, epsilon = TRUE, effsize = FALSE,
      na.omit = TRUE, plot = FALSE, point.size = 4, adjust = TRUE,
      error.width = 0.1, xlab = NULL, ylab = NULL, ylim = NULL,
      breaks = ggplot2::waiver(), jitter = TRUE, line = TRUE,
      jitter.size = 1.25, jitter.width = 0.05, jitter.height = 0,
      jitter.alpha = 0.1, title = "", subtitle = "Confidence Interval",
      digits = 2, p.digits = 4, as.na = NULL, 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

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 1: Repeated measures ANOVA
aov.w(cbind(time1, time2, time3) ~ 1, data = dat)

# Example 2: Repeated measures ANOVA
# print results based on all sphericity corrections
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, print = "all")

# Example 3: Repeated measures ANOVA
# print effect size measures
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, effsize = TRUE)

# Example 4: Repeated measures ANOVA
# do not print hypotheses and descriptive statistics,
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, descript = FALSE, hypo = FALSE)

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

# Example 6: Repeated measures ANOVA
# plot results
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, plot = TRUE)

# Load ggplot2 package
library(ggplot2)

# Save plot, ggsave() from the ggplot2 package
ggsave("Repeated_measures_ANOVA.png", dpi = 600, width = 4.5, height = 4)

# Example 7: Repeated measures ANOVA
# extract plot
p <- aov.w(cbind(time1, time2, time3) ~ 1, data = dat, output = FALSE)$plot
p

# Extract data
plotdat <- aov.w(cbind(time1, time2, time3) ~ 1, data = dat, output = FALSE)$data

# Draw plot in line with the default setting of aov.w()
ggplot(plotdat$long, aes(time, y, group = 1L)) +
geom_point(aes(time, y, group = id),
          alpha = 0.1, position = position_dodge(0.05)) +
geom_line(aes(time, y, group = id),
         alpha = 0.1, position = position_dodge(0.05)) +
geom_point(data = plotdat$ci, aes(variable, m), stat = "identity", size = 4) +
stat_summary(aes(time, y), fun = mean, geom = "line") +
geom_errorbar(data = plotdat$ci, aes(variable, m, ymin = low, ymax = upp), width = 0.1) +
theme_bw() + xlab(NULL) +
labs(subtitle = "Two-Sided 95
theme(plot.subtitle = element_text(hjust = 0.5),
     plot.title = element_text(hjust = 0.5))

## End(Not run)

Description

The function as.na replaces user-specified values in the argument na in a vector, factor, matrix, array, list, or data frame with NA, while the function na.as replaces NA in a vector, factor, matrix or data frame with user-specified values in the argument na.

Usage

as.na(..., data = NULL, na, replace = TRUE, check = TRUE)

na.as(..., data = NULL, na, replace = TRUE, as.na = NULL, check = TRUE)

Arguments

Value

Returns a vector, factor, matrix, array, data frame, or list specified in the argument ... or a data frame specified in data with variables specified in ... replaced.

Author(s)

Takuya Yanagida [email protected]

References

Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth & Brooks/Cole.

See Also

na.auxiliary, na.coverage, na.descript, na.indicator, na.pattern, na.prop, na.test

Examples

#-------------------------------------------------------------------------------
# Numeric vector
num <- c(1, 3, 2, 4, 5)

# Example 1: Replace 2 with NA
as.na(num, na = 2)

# Example 2: Replace 2, 3, and 4 with NA
as.na(num, na = c(2, 3, 4))

# Example 3: Replace NA with 2
na.as(c(1, 3, NA, 4, 5), na = 2)

#-------------------------------------------------------------------------------
# Character vector
chr <- c("a", "b", "c", "d", "e")

# Example 4: Replace "b" with NA
as.na(chr, na = "b")

# Example 5: Replace "b", "c", and "d" with NA
as.na(chr, na = c("b", "c", "d"))

# Example 6: Replace NA with "b"
na.as(c("a", NA, "c", "d", "e"), na = "b")

#-------------------------------------------------------------------------------
# Factor
fac <- factor(c("a", "a", "b", "b", "c", "c"))

# Example 7: Replace "b" with NA
as.na(fac, na = "b")

# Example 8: Replace "b" and "c" with NA
as.na(fac, na = c("b", "c"))

# Example 9: Replace NA with "b"
na.as(factor(c("a", "a", NA, NA, "c", "c")), na = "b")

#-------------------------------------------------------------------------------
# Matrix
mat <- matrix(1:20, ncol = 4)

# Example 10: Replace 8 with NA
as.na(mat, na = 8)

# Example 11: Replace 8, 14, and 20 with NA
as.na(mat, na = c(8, 14, 20))

# Example 12: Replace NA with 2
na.as(matrix(c(1, NA, 3, 4, 5, 6), ncol = 2), na = 2)

#-------------------------------------------------------------------------------
# Array

# Example 13: Replace 1 and 10 with NA
as.na(array(1:20, dim = c(2, 3, 2)), na = c(1, 10))

#-------------------------------------------------------------------------------
# List

# Example 14:  Replace 1 with NA
as.na(list(x1 = c(1, 2, 3, 1, 2, 3),
           x2 = c(2, 1, 3, 2, 1),
           x3 = c(3, 1, 2, 3)), na = 1)

#-------------------------------------------------------------------------------
# Data frame
df <- data.frame(x1 = c(1, 2, 3),
                 x2 = c(2, 1, 3),
                 x3 = c(3, 1, 2))

# Example 15a: Replace 1 with NA
as.na(df, na = 1)

# Example 15b: Alternative specification using the 'data' argument
as.na(., data = df, na = 1)

# Example 16: Replace 1 and 3 with NA
as.na(df, na = c(1, 3))

# Example 17a: Replace 1 with NA in 'x2'
as.na(df$x2, na = 1)

# Example 17b: Alternative specification using the 'data' argument
as.na(x2, data = df, na = 1)

# Example 18: Replace 1 with NA in 'x2' and 'x3'
as.na(x2, x3, data = df, na = 1)

# Example 19: Replace 1 with NA in 'x1', 'x2', and 'x3'
as.na(x1:x3, data = df, na = 1)

# Example 20: Replace NA with -99
na.as(data.frame(x1 = c(NA, 2, 3),
                 x2 = c(2, NA, 3),
                 x3 = c(3, NA, 2)), na = -99)

# Example 2: Recode by replacing 30 with NA and then replacing NA with 3
na.as(data.frame(x1 = c(1, 2, 30),
                 x2 = c(2, 1, 30),
                 x3 = c(30, 1, 2)), na = 3, as.na = 30)

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", "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

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 a 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",
           nrow = NULL, ncol = NULL, 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"),
           file = "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", file = "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(c(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(c(-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", "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

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 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 cluster (CWC, i.e., group mean centering).

Usage

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

Arguments

Details

Single-Level Data

Predictor variables in single-level data can only be centered at the grand mean (CGM) by specifying type = "CGM":

$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

Level-1 (L1) predictor variables in two-level data can be centered at the grand mean (CGM) 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.

L1 predictor variables are centered at the group mean (CWC) by specifying type = "CWC" (Default):

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

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

Level-2 (L1) predictor variables in two-level data can only be centered at the grand mean:

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

where $x_{.j}$ is the predictor value of Level 2 cluster $j$ and $\bar{x}_{..}$ is the average Level-2 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

Level-1 (L1) predictor variables in three-level data can be 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 Level-2 cluster $j$ within Level-3 cluster $k$ and $\bar{x}_{...}$ is the average $x$ score.

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

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

or within Level-3 cluster (cwc.mean = "L3", see Enders, 2013):

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

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

Level-2 (L2) predictor variables in three-level data can be centered at the grand mean (CGM) by specifying type = "CGM":

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

where $x_{.jk}$ is the predictor value of Level-2 cluster $j$ within Level-3 cluster $k$ and $\bar{x}_{...}$ is the average Level-2 cluster score.

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

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

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

Level-3 (L3) predictor variables in three-level data can only be centered at the grand mean:

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

where $x_{..k}$ is the predictor value of Level-3 cluster $k$ and $\bar{x}_{...}$ is the average Level-3 cluster score. Note that the cluster membership variable needs to be specified when centering a L3 predictor variable in three-level data.

Value

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

Author(s)

Takuya Yanagida [email protected]

References

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

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

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, rwg.lindell, item.scores.

Examples

#----------------------------------------------------------------------------
# Predictor Variables in Single-Level Data

# Example 1a: Center predictor 'disp' at the grand mean
center(mtcars$disp)

# Example 1b: Alternative specification using the 'data' argument
center(disp, data = mtcars)

# Example 2a: Center predictors 'disp' and 'hp' at the grand mean and append to 'mtcars'
cbind(mtcars, center(mtcars[, c("disp", "hp")]))

# Example 2b: Alternative specification using the 'data' argument
center(disp, hp, data = mtcars)

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

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

#----------------------------------------------------------------------------
# Predictor Variables in Two-Level Data

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

# Example 5a: Center L1 predictor 'y1' within cluster
center(Demo.twolevel$y1, cluster = Demo.twolevel$cluster)

# Example 5b: Alternative specification using the 'data' argument
center(y1, data = Demo.twolevel, cluster = "cluster")

# Example 6: Center L2 predictor 'w2' at the grand mean
center(w1, data = Demo.twolevel, cluster = "cluster")

# Example 6: Center L1 predictor 'y1' within cluster and L2 predictor 'w1' at the grand mean
center(y1, w1, data = Demo.twolevel, cluster = "cluster")

#----------------------------------------------------------------------------
# Predictor Variables in Three-Level Data

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

# Example 7a: Center L1 predictor 'y1' within L2 cluster
center(y1, data = Demo.threelevel, cluster = c("cluster3", "cluster2"))

# Example 7b: Center L1 predictor 'y1' within L3 cluster
center(y1, data = Demo.threelevel, cluster = c("cluster3", "cluster2"), cwc.mean = "L3")

# Example 7b: Center L1 predictor 'y1' within L2 cluster and L2 predictor 'w1' within L3 cluster
center(y1, w1, data = Demo.threelevel, cluster = c("cluster3", "cluster2"))

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

# Estimate linear mixed-effects model with continuous predictors using lme4 package
mod.lmer <- lme4::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 with continuous predictors using nlme package
mod.lme <- nlme::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 with continuous predictors using glmmTMB package
mod.glmmTMB1 <- glmmTMB::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 with continuous predictors using lme4 package
mod.glmer <- lme4::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 with continuous predictors using glmmTMB package
mod.glmmTMB2 <- glmmTMB::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)

Description

This function computes statistical measures for leverage, distance, and influence for linear models estimated by using the lm() function. Mahalanobis distance and hat values are computed for quantifying leverage, standardized leverage-corrected residuals and studentized leverage-corrected residuals are computed for quantifying distance, and Cook's distance and DfBetas are computed for quantifying influence.

Usage

check.outlier(model, check = TRUE, ...)

Arguments

Details

In regression analysis, an observation can be extreme in three major ways (see Darlington & Hayes, p. 484): (1) An observation has high leverage if it has a atypical pattern of values on the predictors, (2) an observation has high distance if its observed outcome value $Y_i$ has a large deviation from the predicted value $\hat{Y}_i$, and (3) an observation has high influence if its inclusion substantially changes the estimates for the intercept and/or slopes.

Value

Returns a data frame with following entries:

Author(s)

Takuya Yanagida [email protected]

References

Darlington, R. B., &, Hayes, A. F. (2017). Regression analysis and linear models: Concepts, applications, and implementation. The Guilford Press.

See Also

check.collin, lm

Examples

# Example 1: Regression model and measures for leverage, distance, and influence
mod.lm <- lm(mpg ~ cyl + disp + hp, data = mtcars)
check.outlier(mod.lm)

# Merge result table with the data
dat1 <- cbind(mtcars, check.outlier(mod.lm))

Description

This function performs residual diagnostics for linear models estimated by using the lm() function for detecting nonlinearity (partial residual or component-plus-residual plots), nonconstant error variance (predicted values vs. residuals plot), and non-normality of residuals (Q-Q plot and histogram with density plot).

Usage

check.resid(model, type = c("linear", "homo", "normal"),
            resid = c("unstand", "stand", "student"),
            point.shape = 21, point.fill = "gray80", point.size = 1,
            line1 = TRUE, line2 = TRUE,
            line.type1 = "solid", line.type2 = "dashed",
            line.width1 = 1, line.width2 = 1,
            line.color1 = "#0072B2", line.color2 = "#D55E00",
            bar.width = NULL, bar.n = 30, bar.color = "black",
            bar.fill = "gray95", strip.size = 11,
            label.size = 10, axis.size = 10,
            xlimits = NULL, ylimits = NULL,
            xbreaks = ggplot2::waiver(), ybreaks = ggplot2::waiver(),
            check = TRUE, plot = TRUE)

Arguments