Between-Subject Analysis of Variance

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, posthoc = FALSE, conf.level = 0.95, hypo = TRUE,
      descript = TRUE, effsize = FALSE, weighted = FALSE, correct = FALSE,
      digits = 2, p.digits = 3, as.na = NULL, plot = FALSE, bar = TRUE,
      point = FALSE, ci = TRUE, jitter = FALSE, adjust = TRUE,
      point.size = 3, errorbar.width = 0.1, jitter.size = 1.25,
      jitter.width = 0.05, jitter.height = 0, jitter.alpha = 0.1,
      xlab = NULL, ylab = "y", ylim = NULL, ybreaks = ggplot2::waiver(),
      title = NULL, subtitle = "Confidence Interval", filename = NULL,
      width = NA, height = NA, units = c("in", "cm", "mm", "px"), dpi = 600,
      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 takuya.yanagida@univie.ac.at

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

# Example 1: Between-subject ANOVA
aov.b(mpg ~ gear, data = mtcars)

# Example 2: Between-subject ANOVA
# print effect size measures
aov.b(mpg ~ gear, data = mtcars, effsize = TRUE)

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

# Example 4: Between-subject ANOVA
# plot results
aov.b(mpg ~ gear, data = mtcars, plot = TRUE)

## Not run: 
# Example 5: Write Results into a text file
aov.b(mpg ~ gear, data = mtcars, write = "ANOVA.txt")

# Example 6: Save plot
aov.b(mpg ~ gear, data = mtcars, plot = TRUE, filename = "Between-Subject_ANOVA.png",
      width = 7, height = 6)

# Example 7: Draw plot in line with the default setting of aov.b()
library(ggplot2)

object <- aov.b(mpg ~ gear, data = mtcars, output = FALSE)

ggplot(object$data, aes(group, y)) +
  geom_bar(stat = "summary", fun = "mean") +
  geom_errorbar(data = ci.mean(object$data, y, group = "group", adjust = TRUE,
                output = FALSE)$result,
                aes(group, m, ymin = low, ymax = upp), width = 0.1) +
  scale_x_discrete(name = NULL) +
  labs(subtitle = "Two-Sided Difference-Adjusted Confidence Interval") +
  theme_bw() + theme(plot.subtitle = element_text(hjust = 0.5))
## End(Not run)

Repeated Measures Analysis of Variance (Within-Subject ANOVA)

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,
      digits = 2, p.digits = 3, as.na = NULL, plot = FALSE, point = TRUE, line = TRUE,
      ci = TRUE, jitter = FALSE, adjust = TRUE, point.size = 3, line.width = 0.5,
      errorbar.width = 0.1, jitter.size = 1.25,jitter.width = 0.05, jitter.alpha = 0.1,
      xlab = NULL, ylab = "y", ylim = NULL, ybreaks = ggplot2::waiver(), title = NULL,
      subtitle = "Confidence Interval", filename = NULL, width = NA, height = NA,
      units = c("in", "cm", "mm", "px"), 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 takuya.yanagida@univie.ac.at

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 of 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)

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

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

# Example 7: Save plot
aov.w(cbind(time1, time2, time3) ~ 1, data = dat, plot = TRUE,
      filename = "Repeated_measures_ANOVA.png", width = 7, height = 6)
## End(Not run)

Create, Run, and Print Blimp Models

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)

Blimp Summary Measures, Convergence and Efficiency Diagnostics

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)

Blimp Trace Plots and Posterior Distribution Plots

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(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)

Print Blimp Output

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)

Run Blimp Models

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)

Blimp Input Updating

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)

Centering Predictor Variables in Single-Level and Multilevel Data

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, ..., 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 clusters (cwc.mean = "L2", Default, see Brincks et al., 2017):

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

or within Level-3 clusters (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 variables at Level 2 and Level 3 need 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 data containing the centered variable(s).

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

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

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

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

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

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

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

# Example 2c: Center L1 predictor 'y1' within cluster and L2 predictor 'w1' at the grand mean
center(Demo.twolevel, y1, w1, 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 3a: Center L1 predictor 'y1' within L2 cluster
center(Demo.threelevel, y1, cluster = c("cluster3", "cluster2"))

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

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

Collinearity Diagnostics

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 takuya.yanagida@univie.ac.at

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: 
# Load lme4, nlme, and glmmTMB package
libraries(lme4, nlme, glmmTMB)

#----------------------------------------------------------------------------
# Linear mixed-effects model

# 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)

Statistical Measures for Leverage, Distance, and Influence

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, append = TRUE, 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 takuya.yanagida@univie.ac.at

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: Statistical measures for leverage, distance, and influence
check.outlier(lm(mpg ~ cyl + disp + hp, data = mtcars))

# Example 2: Append statistical measures to the mtcars data frame
cbind(mtcars,
      check.outlier(lm(mpg ~ cyl + disp + hp, data = mtcars), append = FALSE))

Residual Diagnostics

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"), plot = TRUE,
            point.shape = 21, point.fill = "gray80", point.size = 1,
            line1 = TRUE, line2 = TRUE, linetype1 = "solid",
            linetype2 = "dashed", linewidth1 = 1, linewidth2 = 1,
            line.col1 = "#0072B2", line.col2 = "#D55E00", bar.width = NULL,
            bar.n = 30, bar.col = "black", bar.fill = "gray95",
            strip.text.size = 11, label.size = 10, axis.text.size = 10,
            xlim = NULL, ylim = NULL, xbreaks = ggplot2::waiver(),
            ybreaks = ggplot2::waiver(), check = TRUE)

Arguments

Details

Nonlinearity

The violation of the assumption of linearity implies that the model cannot accurately capture the systematic pattern of the relationship between the outcome and predictor variables. In other words, the specified regression surface does not accurately represent the relationship between the conditional mean values of Y and the Xs. That means the average error E(\varepsilon) is not 0 at every point on the regression surface (Fox, 2015).

In multiple regression, plotting the outcome variable Y against each predictor variable X can be misleading because it does not reflect the partial relationship between Y and X (i.e., statistically controlling for the other Xs), but rather the marginal relationship between Y and X (i.e., ignoring the other Xs). Partial residual plots or component-plus-residual plots should be used to detect nonlinearity in multiple regression. The partial residual for the jth predictor variable is defined as

e_i^{(j)} = b_jX_{ij} + e_i

The linear component of the partial relationship between Y and X_j is added back to the least-squares residuals, which may include an unmodeled nonlinear component. Then, the partial residual e_i^{(j)} is plotted against the predictor variable X_j. Nonlinearity may become apparent when a non-parametric regression smoother is applied.

By default, the function plots each predictor against the partial residuals, and draws the linear regression and the loess smooth line to the partial residual plots.

Nonconstant Error Variance

The violation of the assumption of constant error variance, often referred to as heteroscedasticity, implies that the variance of the outcome variable around the regression surface is not the same at every point on the regression surface (Fox, 2015).

Plotting residuals against the outcome variable Y instead of the predicted values \hat{Y} is not recommended because Y = \hat{Y} + e. Consequently, the linear correlation between the outcome variable Y and the residuals e is \sqrt{1 - R^2} where R is the multiple correlation coefficient. In contrast, plotting residuals against the predicted values \hat{Y} is much easier to examine for evidence of nonconstant error variance as the correlation between \hat{Y} and e is 0. Note that the least-squares residuals generally have unequal variance Var(e_i) = \sigma^2 / (1 - h_i) where h is the leverage of observation i, even if errors have constant variance \sigma^2. The studentized residuals e^*_i, however, have a constant variance under the assumption of the regression model. Residuals are studentized by dividing them by \sigma^2_i(\sqrt{(1 - h_i)} where \sigma^2_i is the estimate of \sigma^2 obtained after deleting the ith observation, and h_i is the leverage of observation i (Meuleman et al, 2015).

By default, the function plots the predicted values against the studentized residuals. It also draws a horizontal line at 0, a loess smooth lines for all residuals as well as separate loess smooth lines for positive and negative residuals.

Non-normality of Residuals

Statistical inference under the violation of the assumption of normally distributed errors is approximately valid in all but small samples. However, the efficiency of least squares is not robust because the least-squares estimator is the most efficient and unbiased estimator only when the errors are normally distributed. For instance, when error distributions have heavy tails, the least-squares estimator becomes much less efficient compared to robust estimators. In addition, error distributions with heavy-tails result in outliers and compromise the interpretation of conditional means because the mean is not an accurate measure of central tendency in a highly skewed distribution. Moreover, a multimodal error distribution suggests the omission of one or more discrete explanatory variables that naturally divide the data into groups (Fox, 2016).

By default, the function plots a Q-Q plot of the unstandardized residuals, and a histogram of the unstandardized residuals and a density plot. Note that studentized residuals follow a t-distribution with n - k - 2 degrees of freedom where n is the sample size and k is the number of predictors. However, the normal and t-distribution are nearly identical unless the sample size is small. Moreover, even if the model is correct, the studentized residuals are not an independent random sample from t_{n - k - 2}. Residuals are correlated with each other depending on the configuration of the predictor values. The correlation is generally negligible unless the sample size is small.

Value

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

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

Fox, J. (2016). Applied regression analysis and generalized linear models (3rd ed.). Sage Publications, Inc.

Meuleman, B., Loosveldt, G., & Emonds, V. (2015). Regression analysis: Assumptions and diagnostics. In H. Best & C. Wolf (Eds.), The SAGE handbook of regression analysis and causal inference (pp. 83-110). Sage.

See Also

check.collin, check.outlier

Examples

## Not run: 

#-------------------------------------------------------------------------------
# Residual diagnostics for a linear model

mod <- lm(Ozone ~ Solar.R + Wind + Temp, data = airquality)

# Example 1: Partial (component-plus-residual) plots
check.resid(mod, type = "linear")

# Example 2: Predicted values vs. residuals plot
check.resid(mod, type = "homo")

# Example 3: Q-Q plot and histogram with density plot
check.resid(mod, type = "normal")

#-------------------------------------------------------------------------------
# Extract data and ggplot2 object

object <- check.resid(mod, type = "linear", plot = FALSE)

# Data frame
object$plotdat

# ggplot object
object$plot

## End(Not run)

Colored and Styled Terminal Output Text

Description

This function adds color and style to output texts on terminals that support 'ANSI' color and highlight codes that can be printed by using the cat function.

Usage

chr.color(x, color = c("black", "red", "green", "yellow", "blue", "violet",
                       "cyan", "white", "gray1", "gray2", "gray3",
                       "b.red", "b.green", "b.yellow", "b.blue", "b.violet",
                       "b.cyan", "b.white"),
          bg = c("none", "black", "red", "green", "yellow", "blue", "violet",
                 "cyan", "white"),
          style = c("regular", "bold", "italic", "underline"), check = TRUE)

Arguments

Value

Returns a character vector.

Note

This function is based on functions provided in the crayon package by Gábor Csárdi.

Author(s)

Takuya Yanagida

References

Csárdi G (2022). crayon: Colored Terminal Output. R package version 1.5.2, https://CRAN.R-project.org/package=crayon

See Also

chr.grep, chr.grepl, chr.gsub, chr.omit, chr.trim, chr.trunc

Examples

## Not run: 

# Example 1:
cat(chr.color("Text in red.", color = "red"))

# Example 2:
cat(chr.color("Text in blue with green background.",
              color = "blue", bg = "yellow"))

# Example 3a:
cat(chr.color("Text in boldface.", style = "bold"))

# Example 3b:
cat(chr.color("Text in boldface and italic.", style = c("bold", "italic")))


## End(Not run)

Multiple Pattern Matching

Description

This function searches for matches to the character vector specified in pattern within each element of the character vector x.

Usage

chr.grep(pattern, x, ignore.case = FALSE, perl = FALSE, value = FALSE,
         fixed = FALSE, useBytes = FALSE, invert = FALSE, check = TRUE)

chr.grepl(pattern, x, ignore.case = FALSE, perl = FALSE, fixed = FALSE,
          useBytes = FALSE, check = TRUE)

Arguments

Value

Returns a integer vector with the indices of the mathces when value = FALSE, character vector containing the matching elements when value = TRUE, or a logical vector when using the chr.grepl function.

Author(s)

Takuya Yanagida

References

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

See Also

chr.color, chr.grepl, chr.gsub, chr.omit, chr.trim, chr.trunc

Examples

chr.vector <- c("James", "Mary", "Michael", "Patricia", "Robert", "Jennifer")

# Example 1: Indices of matching elements
chr.grep(c("am", "er"), chr.vector)

# Example 2: Values of matching elements
chr.grep(c("am", "er"), chr.vector, value = TRUE)

# Example 3: Matching element?
chr.grepl(c("am", "er"), chr.vector)

Multiple Pattern Matching And Replacements

Description

This function is a multiple global string replacement wrapper that allows access to multiple methods of specifying matches and replacements.

Usage

chr.gsub(pattern, replacement, x, recycle = FALSE, check = TRUE, ...)

Arguments

Value

Return a character vector of the same length and with the same attributes as x (after possible coercion to character).

Note

This function was adapted from the mgsub() function in the mgsub package by Mark Ewing (2019).

Author(s)

Mark Ewing

References

Mark Ewing (2019). mgsub: Safe, Multiple, Simultaneous String Substitution. R package version 1.7.1. https://CRAN.R-project.org/package=mgsub

See Also

chr.color, chr.grep, chr.grepl, chr.omit, chr.trim, chr.trunc

Examples

# Example 1: Replace 'the' and 'they' with 'a' and 'we'
chr.vector <- "they don't understand the value of what they seek."
chr.gsub(c("the", "they"), c("a", "we"), chr.vector)

# Example 2: Replace 'heyy' and 'ho' with 'yo'
chr.vector <- c("hey ho, let's go!")
chr.gsub(c("hey", "ho"), "yo", chr.vector, recycle = TRUE)

# Example 3: Replace with regular expressions
chr.vector <- "Dopazamine is not the same as dopachloride or dopastriamine, yet is still fake."
chr.gsub(c("[Dd]opa([^ ]*?mine)","fake"), c("Meta\1","real"), chr.vector)

Omit Strings

Description

This function omits user-specified values or strings from a numeric vector, character vector or factor.

Usage

chr.omit(x, omit = "", na.omit = FALSE, check = TRUE)

Arguments

Value

Returns a numeric vector, character vector or factor with values or strings specified in omit omitted from the vector specified in x.

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

See Also

chr.color, chr.grep, chr.grepl, chr.gsub, chr.trim, chr.trunc

Examples

#-------------------------------------------------------------------------------
# Charater vector
x.chr <- c("a", "", "c", NA, "", "d", "e", NA)

# Example 1: Omit character string ""
chr.omit(x.chr)

# Example 2: Omit character string "" and missing values (NA)
chr.omit(x.chr, na.omit = TRUE)

# Example 3: Omit character string "c" and "e"
chr.omit(x.chr, omit = c("c", "e"))

# Example 4: Omit character string "c", "e", and missing values (NA)
chr.omit(x.chr, omit = c("c", "e"), na.omit = TRUE)

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

# Example 5: Omit values 2 and 4
chr.omit(x.num, omit = c(2, 4))

# Example 6: Omit values 2, 4, and missing values (NA)
chr.omit(x.num, omit = c(2, 4), na.omit = TRUE)

#-------------------------------------------------------------------------------
# Factor
x.factor <- factor(letters[1:10])

# Example 7: Omit factor levels "a", "c", "e", and "g"
chr.omit(x.factor, omit = c("a", "c", "e", "g"))

Trim Whitespace from String

Description

This function removes whitespace from start and/or end of a string

Usage

chr.trim(x, side = c("both", "left", "right"), check = TRUE)

Arguments

Value

Returns a character vector with whitespaces removed from the vector specified in x.

Note

This function is based on the str_trim() function from the stringr package by Hadley Wickham.

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

Wickham, H. (2019). stringr: Simple, consistent wrappers for common string operations. R package version 1.4.0.

See Also

chr.color, chr.grep, chr.grepl, chr.gsub, chr.omit, chr.trunc

Examples

x <- "  string  "

# Example 1: Remove whitespace at both sides
chr.trim(x)

# Example 2: Remove whitespace at the left side
chr.trim(x, side = "left")

# Example 3: Remove whitespace at the right side
chr.trim(x, side = "right")

Truncate a Character Vector to a Maximum Width

Description

This function truncates a character vector, so that the number of characters of each element of the character vector is always less than or equal to the width specified in the argument width.

Usage

chr.trunc(x, width, side = c("right", "left", "center"), ellipsis = "...",
          check = TRUE)

Arguments

Value

Returns a truncated character vector.

Note

This function was adapted from the str_trunc() function in the stringr package by Hadley Wickham (2023).

Author(s)

Takuya Yanagida

References

Wickham H (2023). stringr: Simple, Consistent Wrappers for Common String Operations. R package version 1.5.1, https://CRAN.R-project.org/package=stringr

See Also

chr.color, chr.grep, chr.grepl, chr.gsub, chr.omit, chr.trim

Examples

# Example 1: Truncate at the right side with a max. of 10 characters
chr.trunc(row.names(mtcars), width = 10)

# Example 2: Truncate at the left side with a max. of 10 characters
chr.trunc(row.names(mtcars), width = 10, side = "left")

# Example 3: Truncate without ellipses
chr.trunc(row.names(mtcars), width = 10, ellipsis = "")

(Bootstrap) Confidence Intervals for Correlation Coefficients

Description

This function computes and plots (1) Fisher z' confidence intervals for Pearson product-moment correlation coefficients (a) without non-normality adjustment, (1b) adjusted via sample joint moments method or (1c) adjusted via approximate distribution method (Bishara et al., 2018), (2) Spearman's rank-order correlation coefficients with (2a) Fieller et al. (1957) standard error, (2b) Bonett and Wright (2000) standard error, or (2c) rank-based inverse normal transformation, (3) Kendall's Tau-b, and (4) Kendall-Stuart's Tau-c correlation coefficients with Fieller et al. (1957) standard error, optionally by a grouping and/or split variable. The function also supports five types of bootstrap confidence intervals (e.g., bias-corrected (BC) percentile bootstrap or bias-corrected and accelerated (BCa) bootstrap confidence intervals) and plots the bootstrap samples with histograms and density curves. By default, the function computes Pearson product-moment correlation coefficients adjusted via approximate distribution method.

Usage

ci.cor(data, ...,
       method = c("pearson", "spearman", "kendall-b", "kendall-c"),
       adjust = c("none", "joint", "approx"),
       se = c("fisher", "fieller", "bonett", "rin"),
       sample = TRUE, seed = NULL, maxtol = 1e-05, nudge = 0.001,
       boot = c("none", "norm", "basic", "perc", "bc", "bca"), R = 1000,
       fisher = TRUE, alternative = c("two.sided", "less", "greater"),
       conf.level = 0.95, group = NULL, split = NULL, na.omit = FALSE, digits = 2,
       as.na = NULL, plot = c("none", "ci", "boot"), point.size = 2.5,
       point.shape = 19, errorbar.width = 0.3, dodge.width = 0.5, hist = TRUE,
       binwidth = NULL, bins = NULL, hist.alpha = 0.4, fill = "gray85", density = TRUE,
       density.col = "#0072B2", density.linewidth = 0.5, density.linetype = "solid",
       point = TRUE, point.col = "#CC79A7", point.linewidth = 0.6,
       point.linetype = "solid", ci = TRUE, ci.col = "black",
       ci.linewidth = 0.6, ci.linetype = "dashed", line = TRUE, intercept = 0,
       linetype = "solid", line.col = "gray65", xlab = NULL, ylab = NULL,
       xlim = NULL, ylim = NULL, xbreaks = ggplot2::waiver(), ybreaks = ggplot2::waiver(),
       axis.title.size = 11, axis.text.size = 10, strip.text.size = 11, title = NULL,
       subtitle = NULL, group.col = NULL, plot.margin = NA,  legend.title = "",
       legend.position = c("right", "top", "left", "bottom", "none"),
       legend.box.margin = c(-10, 0, 0, 0), facet.ncol = NULL, facet.nrow = NULL,
       facet.scales = "free_y", filename = NULL, width = NA, height = NA,
       units = c("in", "cm", "mm", "px"), dpi = 600, write = NULL,
       append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

Pearson Product-Moment Correlation Coefficient

The Fisher z' confidence interval method for the Pearson product-moment correlation coefficient is based on the assumption that X and Y have a bivariate normal distribution in the population. Non-normality resulting from either high kurtosis or high absolute skewness can distort the Fisher z' confidence interval that produces a coverage rate that does not equal the one intended. The distortion is largest when population correlation is large and both variables X and Y were non-normal (Bishara et al., 2017). Note that increasing sample size improves coverage only when the population correlation is zero, while increasing the sample size worsens coverage with a non-zero population correlation (Bishara & Hittner, 2017). The ci.cor function computes the Fisher z' confidence interval without non-normality adjustment (adjust = "none"), with non-normality adjustment via sample joint moments (adjust = "joint"), or with non-normality adjustment via approximate distribution (adjust = "approx"):

• Fisher z' confidence interval method uses the r-toz' transformation for the correlation coefficient r:

  z' = 0.5\cdot \ln\left(\frac{1 + r}{1 - r}\right)

  The sampling distribution of z is approximately normal with a standard error of approximately

  \sigma_z' = \sqrt{\frac{1}{(n - 3)}}

  The two-sided 95% confidence interval is defined as

  z' \pm 1.96\cdot\sigma_{z'}

  These confidence interval bounds are transformed back to the scale of r:

  r = \frac{exp(2z') - 1}{exp(2z') + 1}

  The resulting confidence interval of the correlation coefficient is an approximation and is only accurate when X and Y have a bivariate normal distribution in the population or when the population correlation is zero.

• The Joint Moments Method multiplies the asymptotic variance of z' by \tau_f^2 (Hawkins, 1989):

  \tau_f^2 = \frac{(\mu_{40} + 2\mu_{22} + \mu_{04})\rho^2 - 4(\mu_{31} + \mu_{13})\rho + 4\mu_{22})}{4(1 - \rho^2)^2}

  where \mu_{jk} represents a population joint moment defined as

  \mu_{jk} = E[X^jY^k]

  where X and Y are assumed to be standardized (\mu_{10} = \mu_{01} = 0, \mu_{20} = \mu_{02} = 1). The standard error of z' can then be approximated as \tilde{\sigma}_{z'}:

  \tilde{\sigma}_{z'} = \tau_f\sqrt{\frac{1}{n - 3}}

  The corresponding sample moments, m_{jk} can be used to estimate \tau_f^2:

  \hat{\mu}_{jk} = m_{jk} = \frac{1}{n}\sum_{i=1}^{n}(x_i^jy_i^k)

  However, the higher-order sample joint moments may be unstable estimators of their population counterparts unless the sample size is extremely large. Thus, this estimate of \tau_f^2 may be inaccurate, leading to inaccurate confidence intervals.

• The Approximate Distribution Method estimates an approximate distribution that the sample appears to be drawn from to analytically solve for \tau_f^2 based on that distribution's parameters. The ci.cor function uses a third-order polynomial family allowing estimation of distribution parameters using marginal skewness and kurtosis that are estimated using the marginal sample skewness and kurtosis statistics (Bishara et al., 2018).

Bishara et al. (2018) conducted two Monte Carlo simulations that showed that the approximate distribution method was effective in dealing with violations of the bivariate normality assumption for a wide range of sample sizes, while the joint moments method was effective mainly when the sample size was extremely large, in the thousands. However, the third-order polynomial family used for the approximate distribution method cannot deal with absolute skewness above 4.4 or kurtosis above 43.4. Note that the approximate distribution method is accurate even when the bivariate normality assumption is satisfied, while the sample joint moments method sometimes fails to achieve the intended coverage even when the bivariate normality was satisfied.

Spearman's Rank-Order Correlation Coefficient

The confidence interval for Spearman's rank-order correlation coefficient is based on the Fisher's z method (se = "fisher"), Fieller et al. (1957) approximate standard error (se = "fieller", default), Bonett and Wright (2000) approximate standard error (se = "bonett") or rank-based inverse normal (RIN) transformation (se = "rin") :

• Fisher's z Standard Error

  \sqrt{\frac{1}{(n - 3)}}

• Fieller et al. (1957) Approximate Standard Error

  \sqrt{\frac{1.06}{(n - 3)}}

  Note that this approximation for the standard error is recommended for n > 10 and |rs| < 0.8.

• Bonett and Wright (2000) Approximate Standard Error

  \sqrt{\frac{1 + \frac{\hat{\theta}^2}{2}}{(n - 3)}}

  where \hat{\theta} is the point estimate of the Spearman's rank-order correlation coefficient. Note that this approximation for the standard error is recommended for |\tau| \le 0.9.

• Rin Transformation involves three steps. First, the variable is converted to ranks. Second, the ranks are converted to a 0-to-1 scale using a linear function. Third, this distribution is transformed via the inverse of the normal cumulative distribution function (i.e., via probit transformation). The result is an approximately normal distribution regardless of the original shape of the data, so long as ties are infrequent and n is not too small.

Kendall's Tau-b and Tau-c Correlation Coefficient

The confidence interval for Kendall's Tau-b and Tau-c correlation coefficient is based on the approximate standard error by Fieller et al. (1957):

\sigma_z' = \sqrt{\frac{0.437}{(n - 4)}}

Note that this approximation for the standard error is recommended for n > 10 and |\tau| < 0.8.

Bootstrap Confidence Intervals

The ci.cor function supports bootstrap confidence intervals (CI) for the correlation coefficient by changing the default setting boot = "none" to request one of five different types of bootstrap CI (see Efron & Tibshirani, 1993; Davidson & Hinkley, 1997):

• "norm": The bias-corrected normal approximation bootstrap CI relies on the normal distribution based on the standard deviation of the bootstrap samples \hat{\mathit{SE}}^*. The function corrects for the bootstrap bias, i.e., difference between the bootstrap estimate \hat{\theta}^* and the sample statistic \hat{\theta} centering the interval at 2\hat{\theta} - \hat{\theta}^*. The BC normal CI of intended coverage of 1 - 2(\alpha/2) is given by

  Normal: (\hat{\theta}_{low}, \hat{\theta}_{upp} = \hat{\theta} - (\hat{\theta}^* - \hat{\theta}) + z^{\alpha/2} \cdot \hat{\mathit{SE}}^*, \hat{\theta} - (\hat{\theta}^* - \hat{\theta}) + z^{1 - \alpha/2} \cdot \hat{\mathit{SE}}^*

  where z^{\alpha/2} and z^{1 - \alpha/2} denotes the \alpha and the 1 - \alpha quantile from the standard normal distribution.

• "basic": The basic bootstrap (aka reverse bootstrap percentile) CI is based on the distribution of \hat{\delta} = \hat{\theta} - \theta which is approximated with the bootstrap distribution of \hat{\delta}^* = \hat{\theta}^* - \hat{\theta}.

  Basic: (\hat{\theta}_{low}, \hat{\theta}_{upp} = \hat{\theta} - \hat{\delta}^{*1 - (\alpha/2)}, \hat{\theta} - \hat{\delta}^{*\alpha/2} = 2\hat{\theta} - \hat{\theta}^{*(1 - \alpha/2)} , 2\hat{\theta} - \hat{\theta}^{*(\alpha/2)})

• "perc": The percentile bootstrap CI is computed by ordering the bootstrap estimates \hat{\theta}^*_1, \ldots, \hat{\theta}^*_B to determine the (100(\alpha)/2)th and (100(1 - \alpha)/2)th empirical percentile with intended coverage of 1 - 2(\alpha/2):

  Percentile: (\hat{\theta}_{low}, \hat{\theta}_{upp} = \hat{\theta}^{*(1 - \alpha/2)}, \hat{\theta}^{*(\alpha/2)})

• "bc" (default): The bias-corrected (BC) percentile bootstrap CI corrects the percentile bootstrap CI for median bias of \hat{\theta^*}, i.e., the discrepancy between the median of \hat{\theta}^* and \hat{\theta} in normal units. The bias correction \hat{z}_0 is obtained from the proportion of bootstrap replications less than the sample estimate \hat{\theta}:

  \hat{z}_0 = \Phi^{-1}\left(\frac{\#{\hat{\theta}^*_b < \hat{\theta}}}{B}\right)

  where \Phi^{-1}(.) represents the inverse function of the standard normal cumulative distribution function and B is the number of bootstrap replications. The BC percentile CI of intended coverage of 1 - 2(\alpha/2) is given by

  BC: (\hat{\theta}_{low}, \hat{\theta}_{upp} = \hat{\theta}^{*(\alpha_1)}, \hat{\theta}^{*(\alpha_2)})

  where

  \alpha_1 = \Phi(2\hat{z}_0 + z^{\alpha/2})

  \alpha_2 = \Phi(2\hat{z}_0 + z^{1 - \alpha/2})

  where \Phi(.) represents the standard normal cumulative distribution function and z^{\alpha/2} is the 100(\alpha/2) percentile of a standard normal distribution.

• "bca": The bias-corrected and accelerated (BCa) bootstrap CI corrects the percentile bootstrap CI for median bias \hat{z}_0 and for acceleration or skewness \hat{a}, i.e., the rate of change of the standard error of \hat{\theta} with respect to the true parameter value \theta on a normalized scale. The standard normal approximation \hat{\theta} \sim N(\theta, \mathit{SE}^2) assumes that the standard error of \hat{\theta} is the same for all \theta. The acceleration constant \hat{a} corrects for this unrealistic assumption and can be computed by using jackknife resampling:

  \hat{a} = \frac{\sum_{i=1}^{n}(\hat{\theta}_{(.)} - \hat{\theta}_{(i)})^3}{6\{\sum_{i=1}^{n}(\theta_{(.)} - \hat{\theta}_{(i)})^2\}^{3/2}}

  where \hat{\theta}_{(i)} is the sample estimate with the ith observation deleted and \hat{\theta}_{(.)} = \sum_{i=1}^{n}\frac{\hat{\theta}_{(i)}}{n}. Note that the function uses infinitesimal jackknife instead of regular leave-one-out jackknife that down-weights each observation by an infinitesimal amount of \frac{0.001}{n} instead of removing observations. The BCa percentile CI of intended coverage of 1 - 2(\alpha/2) is given by

  BCa: (\hat{\theta}_{low}, \hat{\theta}_{upp} = \hat{\theta}^{*(\alpha_1)}, \hat{\theta}^{*(\alpha_2)})

  where

  \alpha_1 = \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + z^{\alpha/2}}{1 - \hat{a}(\hat{z}_0 + z^{\alpha/2})}\right)

  \alpha_2 = \Phi\left(\hat{z}_0 + \frac{\hat{z}_0 + z^{1 - \alpha/2}}{1 - \hat{a}(\hat{z}_0 + z^{1 - \alpha/2})}\right)

Note that Fisher transformation is applied before computing the confidence intervals to reverse-transform the limits of the interval using the inverse of the Fisher z transformation (fisher = TRUE) when specifying "norm" or "basic" for the argument boot. In addition, interpolation on the normal quantile scale is applied for "basic", "perc", "bc", and "bca" when a non-integer order statistic is required (see equation 5.8 in Davison & Hinkley, 1997).

Value

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

Note

This function is based on a modified copy of the functions provided in the supporting information in Bishara et al. (2018) for the sample joint moments method and approximate distribution method, functions provided in the supplementary materials in Bishara and Hittner (2017) for Fieller et al. (1957) and Bonett and Wright (2000) correction, and a function provided by Thom Baguley (2024) for the rank-based inverse normal (RIN) transformation. Bootstrap confidence intervals are computed using the R package boot by Angelo Canty and Brain Ripley (2024).

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

Baguley, T. (2024). CI for Spearman's rho. seriousstats. https://rpubs.com/seriousstats/616206

Bonett, D. G., & Wright, T. A. (2000). Sample size requirements for estimating Pearson, Kendall and Spearman correlations. Psychometrika, 65, 23-28. https://doi.org/10.1007/BF02294183

Bishara, A. J., & Hittner, J. B. (2012). Testing the significance of a correlation with nonnormal data: Comparison of Pearson, Spearman, transformation, and resampling approaches. Psychological Methods, 17(3), 399–417. https://doi.org/10.1037/a0028087

Bishara, A. J., & Hittner, J.B. (2017). Confidence intervals for correlations when data are not normal. Behavior Research Methods, 49, 294-309. https://doi.org/10.3758/s13428-016-0702-8

Bishara, A. J., Li, J., & Nash, T. (2018). Asymptotic confidence intervals for the Pearson correlation via skewness and kurtosis. British Journal of Mathematical and Statistical Psychology, 71(1), 167–185. https://doi.org/10.1111/bmsp.12113

Brown, M. B., & Benedetti, J. K. (1977). Sampling behavior of tests for correlation in two-way contingency tables. Journal of the American Statistical Association, 72 (358), 309-315. https://doi.org/10.1080/01621459.1977.10480995

Canty, A., & Ripley, B. (2024). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-31.

Davison, A. C., & Hinkley, D. V. (1997). Bootstrap methods and their application. Cambridge University Press.

Efron, B., & Tibshirani, R. (1993). An introduction to the bootstrap. Chapman & Hall.

Fieller, E. C., Hartley, H. O., & Pearson, E. S. (1957). Tests for rank correlation coefficients: I. Biometrika, 44, 470-481. https://doi.org/10.2307/2332878

Fisher, R. A. (1921). On the “Probable Error” of a Coefficient of Correlation Deduced from a Small Sample. Metron, 1, 3-32.

Hawkins, D. L. (1989). Using U statistics to derive the asymptotic distribution of Fisher’s Z statistic. American Statistician, 43, 235–237. https://doi.org/10.2307/2685369

Hollander, M., Wolfe, D. A., & Chicken, E. (2015). Nonparametric statistical methods. Wiley.

See Also

cor.matrix, ci.mean, ci.mean.diff, ci.prop, ci.var, ci.sd

Examples

#----------------------------------------------------------------------------
# Pearson product-moment correlation coefficient

# Example 1a: Approximate distribution method
ci.cor(mtcars, mpg, drat, qsec)

# Alternative specification without using the '...' argument
ci.cor(mtcars[, c("mpg", "drat", "qsec")])

# Example 1b: Joint moments method
ci.cor(mtcars, mpg, drat, qsec, adjust = "joint")

#----------------------------------------------------------------------------
# Spearman's rank-order correlation coefficient

# Example 2a: Fieller et al. (1957) approximate standard error
ci.cor(mtcars, mpg, drat, qsec, method = "spearman")

# Example 2b: Bonett and Wright (2000) approximate standard error
ci.cor(mtcars, mpg, drat, qsec, method = "spearman", se = "bonett")

# Example 2c: Rank-based inverse normal (RIN) transformation
ci.cor(mtcars, mpg, drat, qsec, method = "spearman", se = "rin")

#----------------------------------------------------------------------------
# Kendall's Tau

# Example 3a:  Kendall's Tau-b
ci.cor(mtcars, mpg, drat, qsec, method = "kendall-b")

# Example 3b:  Kendall's Tau-c
ci.cor(mtcars, mpg, drat, qsec, method = "kendall-c")

## Not run: 
#----------------------------------------------------------------------------
# Bootstrap Confidence Interval (CI)

# Example 4a: Bias-corrected (BC) percentile bootstrap CI
ci.cor(mtcars, mpg, drat, qsec, boot = "bc")

# Example 4b: Bias-corrected and accelerated (BCa) bootstrap CI,
# 5000 bootstrap replications, set seed of the pseudo-random number generator
ci.cor(mtcars, mpg, drat, qsec, boot = "bca", R = 5000, seed = 123)

#----------------------------------------------------------------------------
# Grouping and Split Variable

# Example 5a: Grouping variable
ci.cor(mtcars, mpg, drat, qsec, group = "vs")

# Alternative specification without using the argument '...'
ci.cor(mtcars[, c("mpg", "drat", "qsec")], group = mtcars$vs)

# Example 5b: Split variable
ci.cor(mtcars, mpg, drat, qsec, split = "am")

# Alternative specification without using the argument '...'
ci.cor(mtcars[, c("mpg", "drat", "qsec")], split = mtcars$am)

# Example 5c: Grouping and split variable
ci.cor(mtcars, mpg, drat, qsec, group = "vs", split = "am")

# Alternative specification without using the argument '...'
ci.cor(mtcars[, c("mpg", "drat", "qsec")], group = mtcars$vs, split = mtcars$am)

#----------------------------------------------------------------------------
# Write Output

# Example 6a: Text file
ci.cor(mtcars, mpg, drat, qsec, write = "CI_Cor_Text.txt")

# Example 6b: Excel file
ci.cor(mtcars, mpg, drat, qsec, write = "CI_Cor_Excel.xlsx")

#----------------------------------------------------------------------------
# Plot Confidence Intervals

# Example 7a: Pearson product-moment correlation coefficient
ci.cor(mtcars, mpg, drat, qsec, plot = "ci")

# Example 7b: Grouping variable
ci.cor(mtcars, mpg, drat, qsec, group = "vs", plot = "ci")

# Example 7c: Split variable
ci.cor(mtcars, mpg, drat, qsec, split = "am", plot = "ci")

# Example 7d: Save plot as PDF file
ci.cor(mtcars, mpg, drat, qsec, plot = "ci", saveplot = "CI_Cor.pdf",
       width = 8, height = 6)

# Example 7e: Save plot as PNG file
ci.cor(mtcars, mpg, drat, qsec, plot = "ci", saveplot = "CI_Cor.png",
       width = 8, height = 6)

#----------------------------------------------------------------------------
# Plot Bootstrap Samples

# Example 8a: Pearson product-moment correlation coefficient
ci.cor(mtcars, mpg, drat, qsec, boot = "bc", plot = "boot")

# Example 8b: Grouping variable
ci.cor(mtcars, mpg, drat, qsec, group = "vs", boot = "bc", plot = "boot")

# Example 8c: Split variable
ci.cor(mtcars, mpg, drat, qsec, split = "am", boot = "bc", plot = "boot")

# Example 8d: Save plot as PDF file
ci.cor(mpg, drat, qsec, data = mtcars, boot = "bc", plot = "boot",
       saveplot = "CI_Cor_Boot.pdf", width = 14, height = 9)

# Example 8e: Save plot as PNG file
ci.cor(mtcars, mpg, drat, qsec, boot = "bc", plot = "boot",
       saveplot = "CI_Cor_Boot.png", width = 14, height = 9)
## End(Not run)

(Bootstrap) Confidence Intervals for Arithmetic Means and Medians

Description

The function ci.mean computes and plots confidence intervals for arithmetic means with known or unknown population standard deviation or population variance and the function ci.median computes confidence intervals for medians, optionally by a grouping and/or split variable. These functions also supports six types of bootstrap confidence intervals (e.g., bias-corrected (BC) percentile bootstrap or bias-corrected and accelerated (BCa) bootstrap confidence intervals) and plots the bootstrap samples with histograms and density curves.

Usage

ci.mean(data, ..., sigma = NULL, sigma2 = NULL, adjust = FALSE,
        boot = c("none", "norm", "basic", "stud", "perc", "bc", "bca"),
        R = 1000, seed = NULL, sample = TRUE,
        alternative = c("two.sided", "less", "greater"),
        conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
        na.omit = FALSE, digits = 2, as.na = NULL,
        plot = c("none", "ci", "boot"), point.size = 2.5, point.shape = 19,
        errorbar.width = 0.3, dodge.width = 0.5, hist = TRUE,
        binwidth = NULL, bins = NULL, hist.alpha = 0.4, fill = "gray85", density = TRUE,
        density.col = "#0072B2", density.linewidth = 0.5, density.linetype = "solid",
        point = TRUE, point.col = "#CC79A7", point.linewidth = 0.6,
        point.linetype = "solid", ci = TRUE, ci.col = "black",
        ci.linewidth = 0.6, ci.linetype = "dashed", line = FALSE, intercept = 0,
        linetype = "solid", line.col = "gray65", xlab = NULL, ylab = NULL,
        xlim = NULL, ylim = NULL, xbreaks = ggplot2::waiver(),
        ybreaks = ggplot2::waiver(), axis.title.size = 11, axis.text.size = 10,
        strip.text.size = 11, title = NULL, subtitle = NULL, group.col = NULL,
        plot.margin = NA, legend.title = "",
        legend.position = c("right", "top", "left", "bottom", "none"),
        legend.box.margin = c(-10, 0, 0, 0), facet.ncol = NULL, facet.nrow = NULL,
        facet.scales = "free", filename = NULL, width = NA, height = NA,
        units = c("in", "cm", "mm", "px"), dpi = 600, write = NULL,
        append = TRUE, check = TRUE, output = TRUE)

ci.median(data, ..., boot = c("none", "norm", "basic", "stud", "perc", "bc", "bca"),
          R = 1000, seed = NULL, sample = TRUE,
          alternative = c("two.sided", "less", "greater"),
          conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
          na.omit = FALSE, digits = 2, as.na = NULL, plot = c("none", "ci", "boot"),
          point.size = 2.5, point.shape = 19, errorbar.width = 0.3, dodge.width = 0.5,
          hist = TRUE, binwidth = NULL, bins = NULL, hist.alpha = 0.4, fill = "gray85",
          density = TRUE, density.col = "#0072B2", density.linewidth = 0.5,
          density.linetype = "solid", point = TRUE, point.col = "#CC79A7",
          point.linewidth = 0.6, point.linetype = "solid", ci = TRUE, ci.col = "black",
          ci.linewidth = 0.6, ci.linetype = "dashed", line = FALSE, intercept = 0,
          linetype = "solid", line.col = "gray65", xlab = NULL, ylab = NULL,
          xlim = NULL, ylim = NULL, xbreaks = ggplot2::waiver(),
          ybreaks = ggplot2::waiver(), axis.title.size = 11, axis.text.size = 10,
          strip.text.size = 11, title = NULL, subtitle = NULL, group.col = NULL,
          plot.margin = NA,  legend.title = "",
          legend.position = c("right", "top", "left", "bottom", "none"),
          legend.box.margin = c(-10, 0, 0, 0), facet.ncol = NULL, facet.nrow = NULL,
          facet.scales = "free", filename = NULL, width = NA, height = NA,
          units = c("in", "cm", "mm", "px"), dpi = 600, write = NULL, append = TRUE,
          check = TRUE, output = TRUE)

Arguments

Value

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

Note

Bootstrap confidence intervals are computed using the R package boot by Angelo Canty and Brain Ripley (2024).

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

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

Canty, A., & Ripley, B. (2024). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-31.

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

See Also

test.z, test.t, ci.mean.diff, ci.cor, ci.prop, ci.var, ci.sd, descript

Examples

#----------------------------------------------------------------------------
# Confidence Interval (CI) for the Arithmetic Mean

# Example 1a: Two-Sided 95% CI
ci.mean(mtcars)

# Example 1b: Two-Sided 95% Difference-Adjusted CI
ci.mean(mtcars, adjust = TRUE)

# Example 1c: Two-Sided 95% CI with known population standard deviation
ci.mean(mtcars, mpg, sigma = 6)

# Alternative specification without using the '...' argument
ci.mean(mtcars$mpg, sigma = 6)

#----------------------------------------------------------------------------
# Confidence Interval (CI) for the Median

# Example 2a: Two-Sided 95% CI
ci.median(mtcars)

# Example 2b: One-Sided 99% CI
ci.median(mtcars, alternative = "less", conf.level = 0.99)

## Not run: 
#----------------------------------------------------------------------------
# Bootstrap Confidence Interval (CI)

# Example 3a: Bias-corrected (BC) percentile bootstrap CI
ci.mean(mtcars, boot = "bc")

# Example 3b: Bias-corrected and accelerated (BCa) bootstrap CI,
# 5000 bootstrap replications, set seed of the pseudo-random number generator
ci.mean(mtcars, boot = "bca", R = 5000, seed = 123)

#----------------------------------------------------------------------------
# Grouping and Split Variable

# Example 4a: Grouping variable
ci.mean(mtcars, mpg, cyl, disp, group = "vs")

# Alternative specification without using the '...' argument
ci.mean(mtcars[, c("mpg", "cyl", "disp")], group = mtcars$vs)

# Example 4b: Split variable
ci.mean(mtcars, mpg, cyl, disp, split = "am")

# Alternative specification
ci.mean(mtcars[, c("mpg", "cyl", "disp")], split = mtcars$am)

# Example 4c: Grouping and split variable
ci.mean(mtcars, mpg, cyl, disp, group = "vs", split = "am")

# Alternative specification
ci.mean(mtcars[, c("mpg", "cyl", "disp")], group = mtcars$vs, split = mtcars$am)

#----------------------------------------------------------------------------
# Write Output

# Example 5a: Text file
ci.mean(mtcars, write = "CI_Mean_Text.txt")

# Example 5b: Excel file
ci.mean(mtcars, write = "CI_Mean_Excel.xlsx")

#----------------------------------------------------------------------------
# Plot Confidence Intervals

# Example 6a: Two-Sided 95
ci.mean(mtcars, disp, hp, plot = "ci")

# Example 6b: Grouping variable
ci.mean(mtcars, disp, hp, group = "vs", plot = "ci")

# Example 6c: Split variable
ci.mean(mtcars, disp, hp, split = "am", plot = "ci")

# Example 6d: Save plot as PDF file
ci.mean(mtcars, disp, hp, plot = "ci", saveplot = "CI_Mean.pdf",
        width = 9, height = 6)

# Example 6e: Save plot as PNG file
ci.mean(mtcars, disp, hp, plot = "ci", saveplot = "CI_Mean.png",
        width = 9, height = 6)

#----------------------------------------------------------------------------
# Example 7: Plot Bootstrap Samples

# Example 7a: Two-Sided 95
ci.mean(mtcars, disp, hp, boot = "bc", plot = "boot")

# Example 7b: Grouping variable
ci.mean(mtcars, disp, hp, group = "vs", boot = "bc", plot = "boot")

# Example 7c: Split variable
ci.mean(mtcars, disp, hp, split = "am", boot = "bc", plot = "boot")

# Example 7d: Save plot as PDF file
ci.mean(mtcars, disp, hp, boot = "bc", plot = "boot", saveplot = "CI_Mean_Boot.pdf",
        width = 12, height = 7)

# Example 7e: Save plot as PNG file
ci.mean(mtcars, disp, hp, boot = "bc", plot = "boot", saveplot = "CI_Mean_Boot.png",
        width = 12, height = 7)

## End(Not run)

Confidence Interval for the Difference in Arithmetic Means

Description

This function computes a confidence interval for the difference in arithmetic means in a one-sample, two-sample and paired-sample design with known or unknown population standard deviation or population variance for one or more variables, optionally by a grouping and/or split variable.

Usage

ci.mean.diff(x, ...)

## Default S3 method:
ci.mean.diff(x, y, mu = 0, sigma = NULL, sigma2 = NULL,
             var.equal = FALSE, paired = FALSE,
             alternative = c("two.sided", "less", "greater"),
             conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
             digits = 2, as.na = NULL, write = NULL, append = TRUE,
             check = TRUE, output = TRUE, ...)

## S3 method for class 'formula'
ci.mean.diff(formula, data, sigma = NULL, sigma2 = NULL,
             var.equal = FALSE, alternative = c("two.sided", "less", "greater"),
             conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
             na.omit = FALSE, digits = 2, as.na = NULL, write = NULL,
             append = TRUE, check = TRUE, output = TRUE, ...)

Arguments

Value

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

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

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

See Also

test.z, test.t, ci.mean, ci.median, ci.prop, ci.var, ci.sd, descript

Examples

#----------------------------------------------------------------------------
# One-sample design

# Example 1a: Two-Sided 95% CI for 'mpg'
# population mean = 20
ci.mean.diff(mtcars$mpg, mu = 20)

# Example 1a: One-Sided 95% CI for 'mpg'
# population mean = 20
ci.mean.diff(mtcars$mpg, mu = 20, alternative = "greater")

#----------------------------------------------------------------------------
# Two-sample design

# Example 2a: Two-Sided 95% CI for 'mpg' by 'vs'
# unknown population variances, unequal variance assumption
ci.mean.diff(mpg ~ vs, data = mtcars)

# Example 2b: Two-Sided 95% CI for 'mpg' by 'vs'
# unknown population variances, equal variance assumption
ci.mean.diff(mpg ~ vs, data = mtcars, var.equal = TRUE)

# Example 2c: Two-Sided 95% CI for 'mpg' by 'vs'
# known population standard deviations, equal standard deviation assumption
ci.mean.diff(mpg ~ vs, data = mtcars, sigma = 4)

# Example 2d: Two-Sided 95% CI for 'mpg' by 'vs'
# known population standard deviations, unequal standard deviation assumption
ci.mean.diff(mpg ~ vs, data = mtcars, sigma = c(4, 5))

# Example 2e: Two-Sided 95% CI for 'mpg', 'cyl', and 'disp' by 'vs'
# unknown population variances, unequal variance assumption
ci.mean.diff(cbind(mpg, cyl, disp) ~ vs, data = mtcars)

# Example 2f: Two-Sided 95% CI for 'mpg', 'cyl', and 'disp' by 'vs'
# unknown population variances, unequal variance assumption,
# analysis by am separately
ci.mean.diff(cbind(mpg, cyl, disp) ~ vs, data = mtcars, group = mtcars$am)

# Example 2g: Two-Sided 95% CI for 'mpg', 'cyl', and 'disp' by 'vs'
# unknown population variances, unequal variance assumption,
# split analysis by am
ci.mean.diff(cbind(mpg, cyl, disp) ~ vs, data = mtcars, split = mtcars$am)

# Example 2h: Two-Sided 95% CI for the mean difference between 'group1' and 'group2'
# unknown population variances, unequal variance assumption
group1 <- c(3, 1, 4, 2, 5, 3, 6, 7)
group2 <- c(5, 2, 4, 3, 1)

ci.mean.diff(group1, group2)

#----------------------------------------------------------------------------
# Paired-sample design

dat.p <- data.frame(pre = c(1, 3, 2, 5, 7, 6), post = c(2, 2, 1, 6, 8, 9),
                    group = c(1, 1, 1, 2, 2, 2))

# Example 3a: Two-Sided 95% CI for the mean difference in 'pre' and 'post'
# unknown poulation variance of difference scores
ci.mean.diff(dat.p$pre, dat.p$post, paired = TRUE)

# Example 21: Two-Sided 95% CI for the mean difference in 'pre' and 'post'
# unknown poulation variance of difference scores
# analysis by group separately
ci.mean.diff(dat.p$pre, dat.p$post, paired = TRUE, group = dat.p$group)

# Example 22: Two-Sided 95% CI for the mean difference in 'pre' and 'post'
# unknown poulation variance of difference scores
# analysis by group separately
ci.mean.diff(dat.p$pre, dat.p$post, paired = TRUE, split = dat.p$group)

# Example 23: Two-Sided 95% CI for the mean difference in 'pre' and 'post'
# known population standard deviation of difference scores
ci.mean.diff(dat.p$pre, dat.p$post, sigma = 2, paired = TRUE)

Within-Subject Confidence Interval for the Arithmetic Mean

Description

This function computes difference-adjusted Cousineau-Morey within-subject confidence interval for the arithmetic mean.

Usage

ci.mean.w(data, ..., adjust = TRUE,
          alternative = c("two.sided", "less", "greater"), conf.level = 0.95,
          na.omit = TRUE, digits = 2, as.na = NULL, write = NULL, append = TRUE,
          check = TRUE, output = TRUE)

Arguments

Details

The Cousineau within-subject confidence interval (CI, Cousineau, 2005) is an alternative to the Loftus-Masson within-subject CI (Loftus & Masson, 1994) that does not assume sphericity or homogeneity of covariances. This approach removes individual differences by normalizing the raw scores using participant-mean centering and adding the grand mean back to every score:

Y_{ij}^{'} = Y_{ij} - \hat{\mu}_{i} + \hat{\mu}_{grand}

where Y_{ij}^{'} is the score of the ith participant in condition j (for i = 1 to n), \hat{\mu}_{i} is the mean of participant i across all J levels (for j = 1 to J), and \hat{\mu}_{grand} is the grand mean.

Morey (2008) pointed out that Cousineau's (2005) approach produces intervals that are consistently too narrow due to inducing a positive covariance between normalized scores within a condition introducing bias into the estimate of the sample variances. The degree of bias is proportional to the number of means and can be removed by rescaling the confidence interval by a factor of \sqrt{J - 1}/J:

\hat{\mu}_j \pm t_{n - 1, 1 - \alpha/2} \sqrt{\frac{J}{J - 1}} \hat{\sigma}^{'}_{{\hat{\mu}}_j}

where \hat{\sigma}^{'}_{{\mu}_j} is the standard error of the mean computed from the normalized scores of he jth factor level.

Baguley (2012) pointed out that the Cousineau-Morey interval is larger than that for a difference in means by a factor of \sqrt{2} leading to a misinterpretation of these intervals that overlap of 95% confidence intervals around individual means is indicates that a 95% confidence interval for the difference in means would include zero. Hence, following adjustment to the Cousineau-Morey interval was proposed:

\hat{\mu}_j \pm \frac{\sqrt{2}}{2} (t_{n - 1, 1 - \alpha/2} \sqrt{\frac{J}{J - 1}} \hat{\sigma}^{'}_{{\hat{\mu}}_j})

The adjusted Cousineau-Morey interval is informative about the pattern of differences between means and is computed by default (i.e., adjust = TRUE).

Value

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

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

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

Cousineau, D. (2005) Confidence intervals in within-subject designs: A simpler solution to Loftus and Masson’s Method. Tutorials in Quantitative Methods for Psychology, 1, 42–45. https://doi.org/10.20982/tqmp.01.1.p042

Loftus, G. R., and Masson, M. E. J. (1994). Using confidence intervals in within-subject designs. Psychonomic Bulletin and Review, 1, 476–90. https://doi.org/10.3758/BF03210951

Morey, R. D. (2008). Confidence intervals from normalized data: A correction to Cousineau. Tutorials in Quantitative Methods for Psychology, 4, 61–4. https://doi.org/10.20982/tqmp.01.1.p042

See Also

aov.w, test.z, test.t, ci.mean.diff,' ci.median, ci.prop, ci.var, ci.sd, descript

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: Difference-adjusted Cousineau-Morey confidence intervals
ci.mean.w(dat)

# Example 2: Cousineau-Morey confidence intervals
ci.mean.w(dat, adjust = FALSE)

## Not run: 
# Example 3: Write results into a text file
ci.mean.w(dat, write = "WS_Confidence_Interval.txt")
## End(Not run)

(Bootstrap) Confidence Intervals for Proportions

Description

This function computes and plots confidence intervals for proportions, optionally by a grouping and/or split variable. The function also supports three types of bootstrap confidence intervals (e.g., bias-corrected (BC) percentile bootstrap or bias-corrected and accelerated (BCa) bootstrap confidence intervals) and plots the bootstrap samples with histograms and density curves.

Usage

ci.prop(data, ..., method = c("wald", "wilson"),
        boot = c("none", "perc", "bc", "bca"), R = 1000, seed = NULL,
        alternative = c("two.sided", "less", "greater"), conf.level = 0.95,
        group = NULL, split = NULL, sort.var = FALSE, na.omit = FALSE, digits = 3,
        as.na = NULL, plot = c("none", "ci", "boot"), point.size = 2.5,
        point.shape = 19, errorbar.width = 0.3, dodge.width = 0.5, hist = TRUE,
        binwidth = NULL, bins = NULL, hist.alpha = 0.4, fill = "gray85",
        density = TRUE, density.col = "#0072B2", density.linewidth = 0.5,
        density.linetype = "solid", point = TRUE, point.col = "#CC79A7",
        point.linewidth = 0.6, point.linetype = "solid", ci = TRUE, ci.col = "black",
        ci.linewidth = 0.6, ci.linetype = "dashed",  line = FALSE, intercept = 0.5,
        linetype = "solid", line.col = "gray65",  xlab = NULL, ylab = NULL,
        xlim = NULL, ylim = NULL,  xbreaks = ggplot2::waiver(),
        ybreaks = ggplot2::waiver(),  axis.title.size = 11, axis.text.size = 10,
        strip.text.size = 11, title = NULL, subtitle = NULL, group.col = NULL,
        plot.margin = NA, legend.title = "",
        legend.position = c("right", "top", "left", "bottom", "none"),
        legend.box.margin = c(-10, 0, 0, 0), facet.ncol = NULL, facet.nrow = NULL,
        facet.scales = "free_y", filename = NULL, width = NA, height = NA,
        units = c("in", "cm", "mm", "px"), dpi = 600, write = NULL, append = TRUE,
        check = TRUE, output = TRUE)

Arguments

Details

The Wald confidence interval which is based on the normal approximation to the binomial distribution are computed by specifying method = "wald", while the Wilson (1927) confidence interval (aka Wilson score interval) is requested by specifying method = "wilson". By default, Wilson confidence interval is computed which have been shown to be reliable in small samples of n = 40 or less, and larger samples of n > 40 (Brown, Cai & DasGupta, 2001), while the Wald confidence intervals is inadequate in small samples and when p is near 0 or 1 (Agresti & Coull, 1998).

Value

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

Note

Bootstrap confidence intervals are computed using the R package boot by Angelo Canty and Brain Ripley (2024).

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

Agresti, A. & Coull, B.A. (1998). Approximate is better than "exact" for interval estimation of binomial proportions. American Statistician, 52, 119-126.

Brown, L. D., Cai, T. T., & DasGupta, A., (2001). Interval estimation for a binomial proportion. Statistical Science, 16, 101-133.

Canty, A., & Ripley, B. (2024). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-31.

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

Wilson, E. B. (1927). Probable inference, the law of succession, and statistical inference. Journal of the American Statistical Association, 22, 209-212.

See Also

ci.prop, ci.prop.diff, ci.median, ci.prop.diff, ci.cor, ci.var, ci.sd, descript

Examples

#----------------------------------------------------------------------------
# Confidence Interval (CI) for proportions

# Example 1a: Two-Sided 95% CI
ci.prop(mtcars, vs, am)

# Alternative specification without using the '...' argument
ci.prop(mtcars[, c("vs", "am")])

# Example 1b: One-Sided 95% CI using Wald method
ci.prop(mtcars, vs, am, method = "wald", alternative = "less")

## Not run: 
#----------------------------------------------------------------------------
# Bootstrap Confidence Interval (CI)

# Example 2a: Bias-corrected (BC) percentile bootstrap CI
ci.prop(mtcars, vs, am, boot = "bc")

# Example 2b: Bias-corrected and accelerated (BCa) bootstrap CI,
# 5000 bootstrap replications, set seed of the pseudo-random number generator
ci.prop(mtcars, vs, am, boot = "bca", R = 5000, seed = 123)

#----------------------------------------------------------------------------
# Grouping and Split Variable

# Example 3a: Grouping variable
ci.prop(mtcars, vs, group = "am")

# Alternative specification without using the '...' argument
ci.prop(mtcars$vs, group = mtcars$am)

# Example 3b: Split variable
ci.prop(mtcars, vs, split = "am")

# Alternative specification without using the '...' argument
ci.prop(mtcars$vs, split = mtcars$am)

# Example 3c: Grouping and split variable
ci.prop(mtcars, vs, group = "am", split = "cyl")

# Alternative specification without using the '...' argument
ci.prop(mtcars$vs, group = mtcars$am, split = mtcars$cyl)

#----------------------------------------------------------------------------
# Write Output

# Example 4a: Text file
ci.prop(mtcars, vs, am, write = "CI_Prop_Text.txt")

# Example 4b: Excel file
ci.prop(mtcars, vs, am, write = "CI_Prop_Excel.xlsx")

#----------------------------------------------------------------------------
# Plot Confidence Intervals

# Example 5a: Two-Sided 95
ci.prop(mtcars, vs, am, plot = "ci")

# Example 5b: Grouping variable
ci.prop(mtcars, vs, am, group = "am", plot = "ci")

# Example 5c: Split variable
ci.prop(mtcars, vs, am, split = "am", plot = "ci")

# Example 5d: Save plot as PDF file
ci.prop(mtcars, vs, am, plot = "ci", saveplot = "CI_Prop.pdf",
        width = 9, height = 6)

# Example 5e: Save plot as PNG file
ci.prop(mtcars, vs, am, plot = "ci", saveplot = "CI_Prop.png",
        width = 9, height = 6)

#----------------------------------------------------------------------------
# Plot Bootstrap Samples

# Example 6a: Two-Sided 95
ci.prop(mtcars, vs, am, boot = "bc", plot = "boot")

# Example 6b: Grouping variable
ci.prop(mtcars, vs, am, group = "am", boot = "bc", plot = "boot")

# Example 6c: Split variable
ci.prop(mtcars, vs, am, split = "am", boot = "bc", plot = "boot")

# Example 6d: Save plot as PDF file
ci.prop(mtcars, vs, am, boot = "bc", plot = "boot",
        saveplot = "CI_Prop_Boot.pdf", width = 9, height = 6)

# Example 6e: Save plot as PNG file
ci.prop(mtcars, vs, am, boot = "bc", plot = "boot",
        saveplot = "CI_Prop_Boot.png", width = 9, height = 6)

## End(Not run)

Confidence Interval for the Difference in Proportions

Description

This function computes a confidence interval for the difference in proportions in a two-sample and paired-sample design for one or more variables, optionally by a grouping and/or split variable.

Usage

ci.prop.diff(x, ...)

## Default S3 method:
ci.prop.diff(x, y, method = c("wald", "newcombe"), paired = FALSE,
             alternative = c("two.sided", "less", "greater"), conf.level = 0.95,
             group = NULL, split = NULL, sort.var = FALSE, digits = 2,
             as.na = NULL, write = NULL, append = TRUE,
             check = TRUE, output = TRUE, ...)

## S3 method for class 'formula'
ci.prop.diff(formula, data, method = c("wald", "newcombe"),
             alternative = c("two.sided", "less", "greater"), conf.level = 0.95,
             group = NULL, split = NULL, sort.var = FALSE, na.omit = FALSE,
             digits = 2, as.na = NULL, write = NULL, append = TRUE,
             check = TRUE, output = TRUE, ...)

Arguments

Details

The Wald confidence interval which is based on the normal approximation to the binomial distribution are computed by specifying method = "wald", while the Newcombe Hybrid Score interval (Newcombe, 1998a; Newcombe, 1998b) is requested by specifying method = "newcombe". By default, Newcombe Hybrid Score interval is computed which have been shown to be reliable in small samples (less than n = 30 in each sample) as well as moderate to larger samples(n > 30 in each sample) and with proportions close to 0 or 1, while the Wald confidence intervals does not perform well unless the sample size is large (Fagerland, Lydersen & Laake, 2011).

Value

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

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

Fagerland, M. W., Lydersen S., & Laake, P. (2011) Recommended confidence intervals for two independent binomial proportions. Statistical Methods in Medical Research, 24, 224-254.

Newcombe, R. G. (1998a). Interval estimation for the difference between independent proportions: Comparison of eleven methods. Statistics in Medicine, 17, 873-890.

Newcombe, R. G. (1998b). Improved confidence intervals for the difference between binomial proportions based on paired data. Statistics in Medicine, 17, 2635-2650.

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

See Also

ci.prop, ci.mean, ci.mean.diff, ci.median, ci.var, ci.sd, descript

Examples

#----------------------------------------------------------------------------
# Two-sample design

# Example 1a: Two-Sided 95% CI for 'vs' by 'am'
# Newcombes Hybrid Score interval
ci.prop.diff(vs ~ am, data = mtcars)

# Example 1b: Two-Sided 95% CI for 'vs' by 'am'
# Wald CI
ci.prop.diff(vs ~ am, data = mtcars, method = "wald")

# Example 1c: Two-Sided 95% CI for the difference in proportions
# Newcombes Hybrid Score interval
ci.prop.diff(c(0, 1, 1, 0, 0, 1, 0, 1), c(1, 1, 1, 0, 0))

#----------------------------------------------------------------------------
# Paired-sample design

dat.p <- data.frame(pre = c(0, 1, 1, 0, 1), post = c(1, 1, 0, 1, 1))

# Example 2a: Two-Sided 95% CI for the difference in proportions 'pre' and 'post'
# Newcombes Hybrid Score interval
ci.prop.diff(dat.p$pre, dat.p$post, paired = TRUE)

# Example 2b: Two-Sided 95% CI for the difference in proportions 'pre' and 'post'
# Wald CI
ci.prop.diff(dat.p$pre, dat.p$post, method = "wald", paired = TRUE)

(Bootstrap) Confidence Intervals for Variances and Standard Deviations

Description

The function ci.var computes and plots confidence intervals for variances, and the function ci.sd computes confidence intervals for the standard deviations, optionally by a grouping and/or split variable. These functions also supports three types of bootstrap confidence intervals (e.g., bias-corrected (BC) percentile bootstrap or bias-corrected and accelerated (BCa) bootstrap confidence intervals) and plots the bootstrap samples with histograms and density curves.

Usage

ci.var(data, ..., method = c("chisq", "bonett"),
       boot = c("none", "perc", "bc", "bca"), R = 1000, seed = NULL,
       alternative = c("two.sided", "less", "greater"),
       conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
       na.omit = FALSE, digits = 2, as.na = NULL,
       plot = c("none", "ci", "boot"), point.size = 2.5, point.shape = 19,
       errorbar.width = 0.3, dodge.width = 0.5, hist = TRUE,
       binwidth = NULL, bins = NULL, hist.alpha = 0.4, fill = "gray85", density = TRUE,
       density.col = "#0072B2", density.linewidth = 0.5, density.linetype = "solid",
       point = TRUE, point.col = "#CC79A7", point.linewidth = 0.6,
       point.linetype = "solid", ci = TRUE, ci.col = "black",
       ci.linewidth = 0.6, ci.linetype = "dashed", line = FALSE, intercept = 0,
       linetype = "solid", line.col = "gray65", xlab = NULL, ylab = NULL,
       xlim = NULL, ylim = NULL, xbreaks = ggplot2::waiver(), ybreaks = ggplot2::waiver(),
       axis.title.size = 11, axis.text.size = 10, strip.text.size = 11, title = NULL,
       subtitle = NULL, group.col = NULL, plot.margin = NA,  legend.title = "",
       legend.position = c("right", "top", "left", "bottom", "none"),
       legend.box.margin = c(-10, 0, 0, 0), facet.ncol = NULL, facet.nrow = NULL,
       facet.scales = "free", filename = NULL, width = NA, height = NA,
       units = c("in", "cm", "mm", "px"), dpi = 600, write = NULL, append = TRUE,
       check = TRUE, output = TRUE)

ci.sd(data, ..., method = c("chisq", "bonett"),
      boot = c("none", "perc", "bc", "bca"), R = 1000, seed = NULL,
      alternative = c("two.sided", "less", "greater"),
      conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
      na.omit = FALSE, digits = 2, as.na = NULL,
      plot = c("none", "ci", "boot"), point.size = 2.5, point.shape = 19,
      errorbar.width = 0.3, dodge.width = 0.5, hist = TRUE,
      binwidth = NULL, bins = NULL, hist.alpha = 0.4, fill = "gray85", density = TRUE,
      density.col = "#0072B2", density.linewidth = 0.5, density.linetype = "solid",
      point = TRUE, point.col = "#CC79A7", point.linewidth = 0.6,
      point.linetype = "solid", ci = TRUE, ci.col = "black",
      ci.linewidth = 0.6, ci.linetype = "dashed", line = FALSE, intercept = 0,
      linetype = "solid", line.col = "gray65", xlab = NULL, ylab = NULL,
      xlim = NULL, ylim = NULL, xbreaks = ggplot2::waiver(), ybreaks = ggplot2::waiver(),
      axis.title.size = 11, axis.text.size = 10, strip.text.size = 11, title = NULL,
      subtitle = NULL, group.col = NULL, plot.margin = NA,  legend.title = "",
      legend.position = c("right", "top", "left", "bottom", "none"),
      legend.box.margin = c(-10, 0, 0, 0), facet.ncol = NULL, facet.nrow = NULL,
      facet.scales = "free", filename = NULL, width = NA, height = NA,
      units = c("in", "cm", "mm", "px"), dpi = 600, write = NULL,
      append = TRUE, check = TRUE, output = TRUE)

Arguments

Details

The confidence interval based on the chi-square distribution is computed by specifying method = "chisq", while the Bonett (2006) confidence interval is requested by specifying method = "bonett". By default, the Bonett confidence interval interval is computed which performs well under moderate departure from normality, while the confidence interval based on the chi-square distribution is highly sensitive to minor violations of the normality assumption and its performance does not improve with increasing sample size. Note that at least four valid observations are needed to compute the Bonett confidence interval.

Value

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

Note

Bootstrap confidence intervals are computed using the R package boot by Angelo Canty and Brain Ripley (2024).

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

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

Canty, A., & Ripley, B. (2024). boot: Bootstrap R (S-Plus) Functions. R package version 1.3-31.

Bonett, D. G. (2006). Approximate confidence interval for standard deviation of nonnormal distributions. Computational Statistics and Data Analysis, 50, 775-782. https://doi.org/10.1016/j.csda.2004.10.003

See Also

ci.mean, ci.mean.diff, ci.median, ci.prop, ci.prop.diff, ci.cor, descript

Examples

#----------------------------------------------------------------------------
# Confidence Interval (CI) for the Variance

# Example 1a: Two-Sided 95% CI
ci.var(mtcars)

# Example 1b: One-Sided 99% CI based on the chi-square distributio
ci.var(mtcars, alternative = "less", method = "chisq")

#----------------------------------------------------------------------------
# Confidence Interval (CI) for the Standard Deviation

# Example 2a: Two-Sided 95% CI
ci.sd(mtcars)

# Example 2b: One-Sided 99% CI based on the chi-square distributio
ci.sd(mtcars, alternative = "less", method = "chisq")

## Not run: 
#----------------------------------------------------------------------------
# Bootstrap Confidence Interval (CI)

# Example 3a: Bias-corrected (BC) percentile bootstrap CI
ci.var(mtcars, boot = "bc")

# Example 3b: Bias-corrected and accelerated (BCa) bootstrap CI,
# 5000 bootstrap replications, set seed of the pseudo-random number generator
ci.var(mtcars, boot = "bca", R = 5000, seed = 123)

#----------------------------------------------------------------------------
# Grouping and Split Variable

# Example 4a: Grouping variable
ci.var(mtcars, mpg, cyl, disp, group = "vs")

# Alternative specification without using the '...' argument
ci.var(mtcars[, c("mpg", "cyl", "disp")], group = mtcars$vs)

# Example 4b: Split variable
ci.var(mtcars, mpg, cyl, disp, split = "am")

# Alternative specification without using the '...' argument
ci.var(mtcars[, c("mpg", "cyl", "disp")], split = mtcars$am)

# Example 4c: Grouping and split variable
ci.var(mtcars, mpg, cyl, disp, group = "vs", split = "am")

# Alternative specification without using the '...' argument
ci.var(mtcars[, c("mpg", "cyl", "disp")], group = mtcars$vs, split = mtcars$am)

#----------------------------------------------------------------------------
# Write Output

# Example 5a: Text file
ci.var(mtcars, write = "CI_Var_Text.txt")

# Example 5b: Excel file
ci.var(mtcars, write = "CI_Var_Excel.xlsx")

#----------------------------------------------------------------------------
# Plot Confidence Intervals

# Example 6a: Two-Sided 95
ci.var(mtcars, plot = "ci")

# Example 6b: Grouping variable
ci.var(mtcars, disp, hp, group = "vs", plot = "ci")

# Example 6c: Split variable
ci.var(mtcars, disp, hp, split = "am", plot = "ci")

# Example 6d: Save plot as PDF file
ci.var(mtcars, disp, hp, plot = "ci", saveplot = "CI_Var.pdf",
       width = 9, height = 6)

# Example 6e: Save plot as PNG file
ci.var(mtcars, disp, hp, plot = "ci", saveplot = "CI_Var.png",
       width = 9, height = 6)

#----------------------------------------------------------------------------
# Plot Bootstrap Samples

# Example 7a: Two-Sided 95
ci.var(mtcars, disp, hp, boot = "bc", plot = "boot")

# Example 7b: Grouping variable
ci.var(mtcars, disp, hp, group = "vs", boot = "bc", plot = "boot")

# Example 7c: Split variable
ci.var(mtcars, disp, hp, split = "am", boot = "bc", plot = "boot")

# Example 7d: Save plot as PDF file
ci.var(mtcars, disp, hp, boot = "bc", plot = "boot",
       saveplot = "CI_Var_Boot.pdf", width = 12, height = 7)

# Example 7e: Save plot as PNG file
ci.var(mtcars, disp, hp, boot = "bc", plot = "boot",
       saveplot = "CI_Var_Boot.png", width = 12, height = 7)

## End(Not run)

Clear Console in RStudio

Description

This function clears the console equivalent to Ctrl + L in RStudio on Windows, Mac, UNIX, or Linux operating system.

Usage

clear()

Author(s)

Takuya Yanagida

See Also

restart, setsource

Examples

## Not run: 

# Clear console
clear()

## End(Not run)

Cluster Scores

Description

This function computes group means by default.

Usage

cluster.scores(data, ..., cluster,
               fun = c("mean", "sum", "median", "var", "sd", "min", "max"),
               expand = TRUE, append = TRUE, name = ".a", as.na = NULL,
               check = TRUE)

Arguments

Value

Returns a numeric vector or data frame containing cluster scores with the same length or same number of rows as data if expand = TRUE or with the length or number of rows as length(unique(cluster)) if expand = FALSE.

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

References

Hox, J., Moerbeek, M., & van de Schoot, R. (2018). Multilevel analysis: Techniques and applications (3rd. ed.). Routledge.

Snijders, T. A. B., & Bosker, R. J. (2012). Multilevel analysis: An introduction to basic and advanced multilevel modeling (2nd ed.). Sage Publishers.

See Also

item.scores, multilevel.descript, multilevel.icc

Examples

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

# Example 1: Compute cluster means for 'y1' and expand to match the input 'y1'
cluster.scores(Demo.twolevel, y1, cluster = "cluster", append = FALSE)

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

# Example 2: Compute standard deviation for each cluster
# and expand to match the input x
cluster.scores(Demo.twolevel, cluster = "cluster", fun = "sd")

# Example 3: Compute cluster means without expanding the vector
cluster.scores(Demo.twolevel, cluster = "cluster", expand = FALSE)

# Example 4: Compute cluster means for 'y1' and 'y2' and append to 'Demo.twolevel'
cluster.scores(Demo.twolevel, y1, y2, cluster = "cluster")

# Alternative specification without using the '...' argument
cbind(Demo.twolevel,
      cluster.scores(Demo.twolevel[, c("y1", "y2")], cluster = Demo.twolevel$cluster))

Coding Categorical Variables

Description

This function creates k - 1 variables for a categorical variable with k distinct levels. The coding system available in this function are dummy coding, simple coding, unweighted effect coding, weighted effect coding, repeated coding, forward Helmert coding, reverse Helmert coding, and orthogonal polynomial coding.

Usage

coding(data, ...,
       type = c("dummy", "simple", "effect", "weffect", "repeat",
                "fhelm", "rhelm", "poly"), base = NULL,
       name = c("dum.", "sim.", "eff.", "weff.", "rep.", "fhelm.", "rhelm.", "poly."),
       append = TRUE, as.na = NULL, check = TRUE)

Arguments

Details

Dummy Coding

Dummy or treatment coding compares the mean of each level of the categorical variable to the mean of a baseline group. By default, the first group or factor level is selected as baseline group. The intercept in the regression model represents the mean of the baseline group. For example, dummy coding based on a categorical variable with four groups A, B, C, D makes following comparisons: B vs A, C vs A, and D vs A with A being the baseline group.

Simple Coding

Simple coding compares each level of the categorical variable to the mean of a baseline level. By default, the first group or factor level is selected as baseline group. The intercept in the regression model represents the unweighted grand mean, i.e., mean of group means. For example, simple coding based on a categorical variable with four groups A, B, C, D makes following comparisons: B vs A, C vs A, and D vs A with A being the baseline group.

Unweighted Effect Coding

Unweighted effect or sum coding compares the mean of a given level to the unweighted grand mean, i.e., mean of group means. By default, the first group or factor level is selected as omitted group. For example, effect coding based on a categorical variable with four groups A, B, C, D makes following comparisons: B vs (A, B, C, D), C vs (A, B, C, D), and D vs (A, B, C, D) with A being the omitted group.

Weighted Effect Coding

Weighted effect or sum coding compares the mean of a given level to the weighed grand mean, i.e., sample mean. By default, the first group or factor level is selected as omitted group. For example, effect coding based on a categorical variable with four groups A, B, C, D makes following comparisons: B vs (A, B, C, D), C vs (A, B, C, D), and D vs (A, B, C, D) with A being the omitted group.

Repeated Coding

Repeated or difference coding compares the mean of each level of the categorical variable to the mean of the previous adjacent level. For example, repeated coding based on a categorical variable with four groups A, B, C, D makes following comparisons: B vs A, C vs B, and D vs C.

Foward Helmert Coding

Forward Helmert coding compares the mean of each level of the categorical variable to the unweighted mean of all subsequent level(s) of the categorical variable. For example, forward Helmert coding based on a categorical variable with four groups A, B, C, D makes following comparisons: (B, C, D) vs A, (C, D) vs B, and D vs C.

Reverse Helmert Coding

Reverse Helmert coding compares the mean of each level of the categorical variable to the unweighted mean of all prior level(s) of the categorical variable. For example, reverse Helmert coding based on a categorical variable with four groups A, B, C, D makes following comparisons: B vs A, C vs (A, B), and D vs (A, B, C).

Orthogonal Polynomial Coding

Orthogonal polynomial coding is a form of trend analysis based on polynomials of order k - 1, where k is the number of levels of the categorical variable. This coding scheme assumes an ordered-categorical variable with equally spaced levels. For example, orthogonal polynomial coding based on a categorical variable with four groups A, B, C, D investigates a linear, quadratic, and cubic trends in the categorical variable.

Value

Returns a data frame with k - 1 coded variables or a data frame with the same length or same number of rows as ... containing the coded variables.

Note

This function uses the contr.treatment function from the stats package for dummy coding and simple coding, a modified copy of the contr.sum function from the stats package for effect coding, a modified copy of the contr.wec function from the wec package for weighted effect coding, a modified copy of the contr.sdif function from the MASS package for repeated coding, a modified copy of the code_helmert_forward function from the codingMatrices for forward Helmert coding, a modified copy of the contr_code_helmert function from the faux package for reverse Helmert coding, and the contr.poly function from the stats package for orthogonal polynomial coding.

Author(s)

Takuya Yanagida takuya.yanagida@univie.ac.at

See Also

rec, item.reverse

Examples

# Example 1: Dummy coding for 'gear', baseline group = 3
coding(mtcars, gear)

# Alternative specification without using the '...' argument
coding(mtcars$gear)

# Example 2: Dummy coding for 'gear', baseline group = 4
coding(mtcars, gear, base = 4)

# Example 3: Effect coding for 'gear', omitted group = 3
coding(mtcars, gear, type = "effect")

# Example 3: Effect coding for 'gear', omitted group = 4
coding(mtcars, gear, type = "effect", base = 4)

# Example 4a: Dummy-coded variable names with prefix "gear3."
coding(mtcars, gear, name = "gear3.")

# Example 4b: Dummy-coded variables named "gear_4vs3" and "gear_5vs3"
coding(mtcars, gear, name = c("gear_4vs3", "gear_5vs3"))

Cohen's d

Description

This function computes Cohen's d for one-sample, two-sample (i.e., between-subject design), and paired-sample designs (i.e., within-subject design) for one or more variables, optionally by a grouping and/or split variable. In a two-sample design, the function computes the standardized mean difference by dividing the difference between means of the two groups of observations by the weighted pooled standard deviation (i.e., Cohen's d_s according to Lakens, 2013) by default. In a paired-sample design, the function computes the standardized mean difference by dividing the mean of the difference scores by the standard deviation of the difference scores (i.e., Cohen's d_z according to Lakens, 2013) by default. Note that by default Cohen's d is computed without applying the correction factor for removing the small sample bias (i.e., Hedges' g).

Usage

cohens.d(x, ...)

## Default S3 method:
cohens.d(x, y = NULL, mu = 0, paired = FALSE, weighted = TRUE, cor = TRUE,
         ref = NULL, correct = FALSE, alternative = c("two.sided", "less", "greater"),
         conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
         digits = 2, as.na = NULL, write = NULL, append = TRUE,
         check = TRUE, output = TRUE, ...)

## S3 method for class 'formula'
cohens.d(formula, data, weighted = TRUE, cor = TRUE, ref = NULL,
         correct = FALSE, alternative = c("two.sided", "less", "greater"),
         conf.level = 0.95, group = NULL, split = NULL, sort.var = FALSE,
         na.omit = FALSE, digits = 2, as.na = NULL, write = NULL, append = TRUE,
         check = TRUE, output = TRUE, ...)

Arguments

Details

Cohen (1988, p.67) proposed to compute the standardized mean difference in a two-sample design by dividing the mean difference by the unweighted pooled standard deviation (i.e., weighted = FALSE).

Glass et al. (1981, p. 29) suggested to use the standard deviation of the control group (e.g., ref = 0 if the control group is coded with 0) to compute the standardized mean difference in a two-sample design (i.e., Glass's \Delta) since the standard deviation of the control group is unaffected by the treatment and will therefore more closely reflect the population standard deviation.

Hedges (1981, p. 110) recommended to weight each group's standard deviation by its sample size resulting in a weighted and pooled standard deviation (i.e., weighted = TRUE, default). According to Hedges and Olkin (1985, p. 81), the standardized mean difference based on the weighted and pooled standard deviation has a positive small sample bias, i.e., standardized mean difference is overestimated in small samples (i.e., sample size less than 20 or less than 10 in each group). However, a correction factor can be applied to remove the small sample bias (i.e., correct = TRUE). Note that the function uses a gamma function for computing the correction factor, while a approximation method is used if computation based on the gamma function fails.

Note that the terminology is inconsistent because the standardized mean difference based on the weighted and pooled standard deviation is usually called Cohen's d, but sometimes called Hedges' g. Oftentimes, Cohen's d is called Hedges' d as soon as the small sample correction factor is applied. Cumming and Calin-Jageman (2017, p.171) recommended to avoid the term Hedges' g , but to report which standard deviation was used to standardized the mean difference (e.g., unweighted/weighted pooled standard deviation, or the standard deviation of the control group) and whether a small sample correction factor was applied.

As for the terminology according to Lakens (2013), in a two-sample design (i.e., paired = FALSE) Cohen's d_s is computed when using weighted = TRUE (default) and Hedges's g_s is computed when using correct = TRUE in addition. In a paired-sample design (i.e., paired = TRUE), Cohen's d_z is computed when using weighted = TRUE, default, while Cohen's d_{rm} is computed when using weighted = FALSE and cor = TRUE, default and Cohen's d_{av} is computed when using weighted = FALSE and cor = FALSE. Corresponding Hedges' g_z, g_{rm}, and g_{av} are computed when using correct = TRUE in addition.

Value

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