R/estimate_mdiff_contrast_bs_examples.R

In rcalinjageman/esci2: Effect Sizes and Confidence Intervals

# Code for testing estimate_mdiff_contrast_bs
estimate_mdiff_contrast_bs.testing <- function() {
  
  # From summary data ---------------------------------------------------
  
  # Test: Example from Bonett word file
  # Raw score effect sizes should be:
  # Estimate                                SE         t       df      p-value        LL        UL
  # Equal Variances Assumed:        -5.35 1.300136 -4.114955 36.00000 0.0002152581 -7.986797 -2.713203
  # Equal Variances Not Assumed:    -5.35 1.300136 -4.114955 33.52169 0.0002372436 -7.993583 -2.706417
  # And for smds:
  # Bonett script gives:
  #                               Estimate        SE        LL         UL
  # Equal Variances Not Assumed: -1.301263 0.3692800 -2.025039 -0.5774878
  # Equal Variances Assumed:     -1.301263 0.3514511 -1.990095 -0.6124317
  # esci won't match Bonett's equal variance exactly due to different
  #  implementation, but will match MBESS
  # without correction, MBESS gives:
  #   effect_size     lower      upper df        se       moe
  # 1   -1.301263 -1.982268 -0.6055378 36 0.3394148 0.6883652
  # 
  means <- c(33.5, 37.9, 38.0, 44.1)
  sds <- c(3.84, 3.84, 3.65, 4.98)
  ns <- c(10,10,10,10)
  contrast <- c(.5, .5, -.5, -.5)
  CI_mdiff_contrast_bs(means, sds, ns, contrast)
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = contrast,
    outcome_variable_name = "Time in maze",
    grouping_variable_name = "Condition",
    assume_equal_variance = FALSE
  )
  estimate
  my_test <- test_mdiff_contrast_bs(estimate)
  plot_esci_test(estimate)
  plot_mdiff_contrast_bs(estimate, y.axis.text = 10, y.axis.title = 16,
                         x.axis.text = 14, x.axis.title = 10
                         )
  

  
  
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = NULL,
    outcome_variable_name = "Time in maze",
    grouping_variable_name = "Condition",
    assume_equal_variance = FALSE
  )
  
  
  # Example from Kline, Table 3.4, worked out in Chapter 7, p. 202
  #  when worked out with non-central t estimated CI
  # Should give d = -0.8528029, [-2.1213, 0.4483] (non-corrected)
  # and MBESS gives:
  #   effect_size     lower     upper df       se      moe
  # 1  -0.8528029 -2.121155 0.4482577 12 0.589636 1.284707

  x1 <- c(9, 12, 13, 15, 16)
  x2 <- c(8, 12, 11, 10, 14)
  x3 <- c(13, 14, 16, 14, 18)
  means <- c(mean(x1), mean(x2), mean(x3))
  sds <- c(sd(x1), sd(x2), sd(x3))
  ns <- c(length(x1), length(x2), length(x3))
  contrast <- c(1, 0, -1)
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = contrast,
    assume_equal_variance = TRUE
  )
  estimate
  plot_mdiff_contrast_bs(estimate)
  
  
  
  # Test: Two-group example from esci ITNS for Bob's class
  # Bonett's script agrees with esci ITNS:
  #                              Estimate        SE        t       df      p-value       LL       UL
  # Equal Variances Assumed:         5.21 0.8992455 5.793746 149.0000 3.951374e-08 3.433079 6.986921
  # Equal Variances Not Assumed:     5.21 0.8166012 6.380104 117.4315 3.646704e-09 3.592826 6.827174
  means <- c(6.88, 12.09)
  sds <- c(4.22, 5.52)
  ns <- c(48, 103)
  contrast <- c(-1, 1)
  group_labels <- c("Small Group", "Big Group")
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    group_labels = group_labels,
    contrast = contrast,
    conf_level = 0.95,
    assume_equal_variance = FALSE
  )
  estimate
  
  # Test 3: Contrast example from esci ITNS
  # Bonett's script agrees with esci ITNS:
  #                              Estimate       SE         t        df   p-value        LL       UL
  # Equal Variances Assumed:          2.6 2.610673 0.9959117 108.00000 0.3215192 -2.574807 7.774807
  # Equal Variances Not Assumed:      2.6 2.812776 0.9243536  64.31133 0.3587576 -3.018643 8.218643
  # esci ITNS with equal variance assumed also reports contrast groups:
  # Comparison = 37.3 [33.641, 40.959]
  # Reference = 34.7 [31.041, 38.359]
  means <- c(37.5, 31.9, 41.2, 33.4, 29.9, 38.3)
  sds <- c(10, 13.5, 14.8, 10, 8.7, 10)
  ns <- c(19, 19, 19, 19, 19, 19)
  contrast <- c(
    "NFree10" = -1/2,
    "AFree10" = -1/2,
    "ADiet10" = 1/2,
    "NFree17" = 1/2
  )
  group_labels = c(
    "NFree10",
    "AFree10",
    "ADiet10",
    "NFree17",
    "AFree17",
    "ADiet17"
  )
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = contrast,
    group_labels = group_labels,
    conf_level = 0.95,
    assume_equal_variance = TRUE
  )
  estimate
  
  
  # From vectors ---------------------------------------------------
  # Test: From vectors with random data, 4 groups
  IV <- c(
    rep("Control", 50),
    rep("Group1", 50),
    rep("Group2", 50),
    rep("Group3", 50)
  )
  IV <- as.factor(IV)
  DV <- c(
    rnorm(n = 50, mean = 0, sd = 1),
    rnorm(n = 50, mean = 0.1, sd = 1),
    rnorm(n = 50, mean = 0.2, sd = 1),
    rnorm(n = 50, mean = 0.5, sd = 1)
  )
  contrast <- c("Control" = -1, "Group1" = 0.5, "Group2" = 0.5)

  estimate <- estimate_mdiff_contrast_bs(
    grouping_variable = IV,
    outcome_variable = DV,
    contrast = contrast,
    save_raw_data = TRUE
  )
  estimate
  
  
  # Test: vectors with NA values and unused group levels
  IV <- c(
    rep("Control", 50),
    rep("Group1", 50),
    rep("Group2", 50),
    rep("Group3", 50)
  )
  IV <- as.factor(IV)
  DV <- c(
    rnorm(n = 50, mean = 0, sd = 1),
    rnorm(n = 50, mean = 0.1, sd = 1),
    rnorm(n = 50, mean = 0.2, sd = 1),
    rnorm(n = 50, mean = 0.5, sd = 1)
  )
  contrast <- c("Control" = -1, "Group1" = 0.5, "Group2" = 0.5)
  IV[45:55] <- NA
  DV[90:110] <- NA
  levels(IV) <- c(levels(IV), "Unused_level")
  estimate <- estimate_mdiff_contrast_bs(
    grouping_variable = IV,
    outcome_variable = DV,
    contrast = contrast,
    save_raw_data = TRUE
  )
  estimate
  
  
  # From data frames ------------------------------------------------------
  # Test from data frame 
  IV <- c(
    rep("Control", 50),
    rep("Group1", 50),
    rep("Group2", 50),
    rep("Group3", 50),
    rep("Group_Bad", 1)
  )
  IV <- as.factor(IV)
  DV <- c(
    rnorm(n = 50, mean = 0, sd = 1),
    rnorm(n = 50, mean = 1, sd = 0.5),
    rnorm(n = 50, mean = 1, sd = 0.5),
    rnorm(n = 50, mean = 0.5, sd = 4),
    1
  )
  IV[45] <- NA
  DV[90:110] <- NA
  contrast <- c("Control" = -1, "Group1" = 0.5, "Group2" = 0.5)

  mydata <- data.frame(
    "my ind variable" = as.factor(IV),
    myoutcomevariable = DV,
    myextra = rnorm(n = 201, mean = 0, sd = 1)
  )
  
  estimate <- estimate_mdiff_contrast_bs(
    mydata, 
    "my ind variable", 
    "myoutcomevariable",
    contrast = contrast,
  )
  estimate
  
  my_test <- test_mdiff_contrast_bs(estimate)
  plot_esci_test(estimate)
  plot_mdiff_contrast_bs(estimate)
  
  estimate <- estimate_mdiff_contrast_bs(
    mydata, 
    "myindvariable", 
    "myoutcomevariable",
    contrast = NULL,
  )
  
  
  estimate <- estimate_mdiff_contrast_bs(
    mydata, myindvariable, myextra,
    contrast = contrast,
  )
  
  estimate <- estimate_mdiff_contrast_bs.jamovi(
    data = mydata, 
    "myindvariable", 
    c("myextra", "myoutcomevariable"),
    contrast = contrast,
  )
  
  
  
  # Errors from Summary Data ------------------------------------------------
  # Reject NAs
  means = c(33.5, 37.9, 38.0, 44.1)
  sds = c(3.84, 3.84, 3.65, 4.98)
  ns = c(10,10,10,10)
  contrast = c(.5, .5, -.5, -.5)
  estimate <- estimate_mdiff_contrast_bs(
    means = c(33.5, 37.9, 38.0, NA),
    sds = sds,
    ns = ns,
    contrast = contrast,
    outcome_variable_name = "Time in maze",
    grouping_variable_name = "Condition",
    assume_equal_variance = FALSE
  )
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = c(3.84, 3.84, NA, 4.98),
    ns = ns,
    contrast = contrast,
    outcome_variable_name = "Time in maze",
    grouping_variable_name = "Condition",
    assume_equal_variance = FALSE
  )
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = c(10,NA,10,10),
    contrast = contrast,
    outcome_variable_name = "Time in maze",
    grouping_variable_name = "Condition",
    assume_equal_variance = FALSE
  )
  
  # Reject ns < 2
  means = c(33.5, 37.9, 38.0, 44.1)
  sds = c(3.84, 3.84, 3.65, 4.98)
  ns = c(10,10,10,10)
  contrast = c(.5, .5, -.5, -.5)
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = c(1,10,10,10),
    contrast = contrast,
    outcome_variable_name = "Time in maze",
    grouping_variable_name = "Condition",
    assume_equal_variance = FALSE
  )

  # Reject bad conf_level or assume_equal_variance
  means = c(33.5, 37.9, 38.0, 44.1)
  sds = c(3.84, 3.84, 3.65, 4.98)
  ns = c(10,10,10,10)
  contrast = c(.5, .5, -.5, -.5)
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = contrast,
    assume_equal_variance = 22
  )
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = contrast,
    conf_level = 1
  )
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    contrast = contrast,
    conf_level = "1"
  )
  
  # Reject misformed contrasts
  means = c(33.5, 37.9, 38.0, 44.1)
  sds = c(3.84, 3.84, 3.65, 4.98)
  ns = c(10,10,10,10)
  group_labels = c("g1", "g2", "g3", "g4")
  contrast = c(.5, .5, -.5, -.5)
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    group_labels = group_labels,
    contrast = c(.5, .5, -.5, -.5, 0)
  )
  
  estimate <- estimate_mdiff_contrast_bs(
    means = means,
    sds = sds,
    ns = ns,
    group_labels = group_labels,
    contrast = c("g1" = -1, "g2" = 1, "bob" = 4)
  )
  
  # Errors from Vector Data ------------------------------------------------
  IV <- c(
    rep("Control", 50),
    rep("Group1", 50),
    rep("Group2", 50),
    rep("Group_Evil", 1),
    rep("Group3", 50),
    rep("Group_Bad", 1)
  )
  IV <- as.factor(IV)
  DV <- c(
    rnorm(n = 50, mean = 0, sd = 1),
    rnorm(n = 50, mean = 1, sd = 0.5),
    rnorm(n = 50, mean = 1, sd = 0.5),
    1,
    rnorm(n = 50, mean = 0.5, sd = 4),
    1
  )
  IV[45] <- NA
  DV[90:110] <- NA
  contrast <- c("Control" = -1, "Group1" = 0.5, "Group2" = 0.5)
  estimate <- estimate_mdiff_contrast_bs(
    grouping_variable = IV,
    outcome_variable = DV,
    contrast = contrast,
    save_raw_data = TRUE
  )
  contrast <- c("Control" = -1, "Group1" = 0.5, "Group2" = 0.5, "No Group" = -1)
  estimate <- estimate_mdiff_contrast_bs(
    grouping_variable = IV,
    outcome_variable = DV,
    contrast = contrast,
    save_raw_data = TRUE
  )
  contrast <- c("Control" = -1, "Group1" = 0.5, "Group2" = 0.5, "Group_Bad" = -1)
  estimate <- estimate_mdiff_contrast_bs(
    grouping_variable = IV,
    outcome_variable = DV,
    contrast = contrast,
    save_raw_data = TRUE
  )
  contrast <- c(-1, 0.5, 0.5, 0, 1)
  estimate <- estimate_mdiff_contrast_bs(
    grouping_variable = IV,
    outcome_variable = DV,
    contrast = contrast,
    save_raw_data = TRUE
  )
}