Knowledge check!

What type of test would you use for each of the following scenarios?

  1. They say that there is a 50/50 chance of getting heads (or tails) when you flip a coin. You want to put it to the test by flipping a coin 100 times to see if it matches those chances. What test would you use?

  2. The meditation app, Headspace, has released a new version of their app and they want to see if there has been an increase in user satisfaction. They survey the same users before and after the update asking them if they’re satisfied with the app (yes or no). How would you test this difference?

  3. A small veterinary clinic wants to test the association between dog breed (pug vs french bulldog) and breathing problems (yes or no). The problem is, they have a small sample size. Which test can they use?

Load packages

p_needed <- c("tidyr",
          "tidyverse",
          "dplyr",
          "haven",
          "ggplot2",
          "readxl",
          "exact2x2",
          "gmodels",
          "WebPower")

packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]

if(length(p_to_install) > 0){
  install.packages(p_to_install)
}

lapply(p_needed, library, character.only = TRUE)

Load dataframe

# read excel files with the 'readxl' package
df <- read_xlsx("/Users/kareenadelrosario/Desktop/local r code/regression_code/data6.xlsx")

# what variables do we have
colnames(df)
##  [1] "...1"                    "PRE_male"               
##  [3] "PRE_female"              "POST_male"              
##  [5] "POST_female"             "POST_male_salary"       
##  [7] "POST_female_salary"      "POST_male_friendly"     
##  [9] "POST_female_friendly"    "POST_male_intelligent"  
## [11] "POST_female_intelligent" "Hire_m"                 
## [13] "Hire_f"                  "Choice_m"               
## [15] "Choice_f"                "itemnum"                
## [17] "partnum"                 "Gender"                 
## [19] "Ide"                     "Political"              
## [21] "Edu"

Descriptives - CHALLENGE!

Calculate the mean and standard deviation of the PRE_female variable, which indicates percent of women chosen

# Calculate mean

# Calculate sd

Let’s visualize the data in R using a bar graph and histogram like we did in Python

We want to plot the mean and sd of PRE_female. ggplot’s bar plots does not like it if you don’t specify both an x and y-axis. For that reason, we have to get creative

bar_plot <- df %>%
  dplyr::select(PRE_female) %>% # we only need this variable
  pivot_longer(PRE_female) %>% # pivot_longer separates out the variable name "PRE_female" and the value, giving us two new cols: name, value
  ggplot(aes(x = name, y = value, fill = name)) + # have to add "fill" in aes() if you want to change the color manually later
  stat_summary(fun=mean, geom="bar", width = .3) + # plot mean as barplot
  stat_summary(fun.data = mean_se, geom = "errorbar", width = .1) # mean_se is an argument that calculates SD (can calculate SE with additional arguments)

print(bar_plot)

Let’s make prettier to match python

If you’re interested in more color options, check out RColorBrewer: https://r-graph-gallery.com/38-rcolorbrewers-palettes.html

# library(RColorBrewer) # more color options

bar_plot + 
  theme_classic() +
  scale_fill_manual(values = "deepskyblue4") +
  theme(axis.title.x=element_blank(), # get rid of x-axis title
                 axis.title.y=element_blank())  # get rid of y-axis title

Density plot

density_plot <- df %>% 
  ggplot (aes(PRE_female)) +
  geom_density() + 
  scale_x_continuous(limits=c(-0.5,1.5)) + # specify range of x-axis so it doesn't get cut off
  geom_vline(aes(xintercept=mean(PRE_female)), # add a line that represents the mean
            color="deepskyblue4", linetype="dashed", size=1) +
  theme_classic()

print(density_plot)

We can layer the density plot on top of a histogram, if that’s more informative

density_plot + geom_histogram(aes(y=..density..), fill = "blue", alpha = 0.4, bins = 9) # bins is the width of the histogram bars

Inspect the data

# What is our range?
range(df$PRE_female)
## [1] 0 1
# What is the variable class?
class(df$PRE_female)
## [1] "numeric"
# What if you want to see the class of variables for each column in the df
glimpse(df)
## Rows: 260
## Columns: 21
## $ ...1                    <dbl> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, …
## $ PRE_male                <dbl> 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1…
## $ PRE_female              <dbl> 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0…
## $ POST_male               <dbl> 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1…
## $ POST_female             <dbl> 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 1…
## $ POST_male_salary        <dbl> 57, 63, 16, 27, 53, 49, 70, 67, 34, 47, 45, 37…
## $ POST_female_salary      <dbl> 88, 44, 21, 18, 44, 72, 50, 66, 50, 85, 64, 68…
## $ POST_male_friendly      <dbl> 95, 66, 34, 26, 76, 98, 89, 96, 87, 88, 65, 28…
## $ POST_female_friendly    <dbl> 73, 79, 32, 30, 90, 99, 90, 93, 88, 94, 78, 64…
## $ POST_male_intelligent   <dbl> 90, 74, 17, 34, 73, 90, 89, 92, 88, 88, 37, 45…
## $ POST_female_intelligent <dbl> 92, 96, 34, 36, 85, 92, 91, 91, 87, 88, 39, 63…
## $ Hire_m                  <dbl> 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1…
## $ Hire_f                  <dbl> 1, 0, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1…
## $ Choice_m                <dbl> 1, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0…
## $ Choice_f                <dbl> 1, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1…
## $ itemnum                 <dbl> 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1…
## $ partnum                 <dbl> 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7…
## $ Gender                  <chr> "Female", "Female", "Female", "Female", "Femal…
## $ Ide                     <dbl> 100, 100, 42, 42, 51, 51, 92, 92, 2, 2, 13, 13…
## $ Political               <chr> "Democrat", "Democrat", "Other", "Other", "Ind…
## $ Edu                     <chr> "College Degree", "College Degree", "High Scho…

One-sample chi-square test

CHALLENGE: Let’s first count how many 0, 1 in PRE_female. How can we do that?

hint: we can use dplyr’s group_by and count functions

chisq.test(x = table(df$PRE_female)) # for one-sample, need to include table function
## 
##  Chi-squared test for given probabilities
## 
## data:  table(df$PRE_female)
## X-squared = 12.938, df = 1, p-value = 0.0003219

Calculate the effect size (Cohen’s w)

\[w = \sqrt\frac{\chi^2} {N} \]

# Calculate the observed frequencies of 0,1
observed_freq <- table(df$PRE_female)

# Perform the Chi-square test
chi_square_result <- chisq.test(observed_freq)

# Total sample size
N <- sum(observed_freq) 

# Calculate effect size using the equation above
effect_size <- sqrt(chi_square_result$statistic / N)

# Print the effect size
print(effect_size)
## X-squared 
## 0.2230769

Independent proportions chi-square

Visualize the data

bar_plot_twogender <- df %>%
  dplyr::select(PRE_female, PRE_male) %>% # we only need these variable
  pivot_longer(cols = c(PRE_female, PRE_male), names_to = "category", values_to = "value") %>% # nesting male and female under one category
  ggplot(aes(x = category, y = value, fill = category)) + # have to add "fill" in aes() if you want to change the color manually later
  stat_summary(fun=mean, geom="bar", position = position_dodge(), width = .3) + # plot mean as barplot
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.9), width = .1) + # mean_se is an argument that calculates SD (can calculate SE with additional arguments)
  scale_fill_manual(values = c("deepskyblue4", "deepskyblue2")) + # how to add two colors. they're in alphabetical order of the variables (e.g., PRE_female comes before PRE_male)
  labs(x = "Gender", y = "Mean") + # add title to x and y-axis (labs = "labels")
  theme_minimal() + # preset theme
  ggtitle("Means of Female and Males in this Dataframe")

bar_plot_twogender

# Option 1: Classic chi-square output
chisq.test(df$PRE_female, df$PRE_male)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  df$PRE_female and df$PRE_male
## X-squared = 0.13286, df = 1, p-value = 0.7155
# Option 2: The kitchen sink
CrossTable(df$PRE_female, df$PRE_male,
           fisher = TRUE,
           chisq = TRUE,
           expected = TRUE,
           prop.c = FALSE, # This stops the column proportions from being displayed
           prop.t = FALSE, # This stops the total proportions from being displayed
           prop.chisq = FALSE, # This stops the chisquared proportions from being displayed
           sresid = TRUE, # Produces standardized residuals
           format = "SPSS") # To see residuals, must use SPSS format
## 
##    Cell Contents
## |-------------------------|
## |                   Count |
## |         Expected Values |
## |             Row Percent |
## |            Std Residual |
## |-------------------------|
## 
## Total Observations in Table:  260 
## 
##               | df$PRE_male 
## df$PRE_female |        0  |        1  | Row Total | 
## --------------|-----------|-----------|-----------|
##             0 |      104  |       55  |      159  | 
##               |  102.127  |   56.873  |           | 
##               |   65.409% |   34.591% |   61.154% | 
##               |    0.185  |   -0.248  |           | 
## --------------|-----------|-----------|-----------|
##             1 |       63  |       38  |      101  | 
##               |   64.873  |   36.127  |           | 
##               |   62.376% |   37.624% |   38.846% | 
##               |   -0.233  |    0.312  |           | 
## --------------|-----------|-----------|-----------|
##  Column Total |      167  |       93  |      260  | 
## --------------|-----------|-----------|-----------|
## 
##  
## Statistics for All Table Factors
## 
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 =  0.2472369     d.f. =  1     p =  0.6190274 
## 
## Pearson's Chi-squared test with Yates' continuity correction 
## ------------------------------------------------------------
## Chi^2 =  0.1328593     d.f. =  1     p =  0.7154857 
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio:  1.139965 
## 
## Alternative hypothesis: true odds ratio is not equal to 1
## p =  0.6907062 
## 95% confidence interval:  0.6553751 1.97583 
## 
## Alternative hypothesis: true odds ratio is less than 1
## p =  0.7361355 
## 95% confidence interval:  0 1.818849 
## 
## Alternative hypothesis: true odds ratio is greater than 1
## p =  0.3568877 
## 95% confidence interval:  0.7128049 Inf 
## 
## 
##  
##        Minimum expected frequency: 36.12692

That’s a lot of information. What do we actually care about?

  • Pearson’s chi-square has been criticized for being susceptible to Type I error. This test also gives you the Yates continuity correction which is a more conservative version of Pearson’s chi-square.

  • Copying what we did in Python, let’s say we had a small sample. Which test statistic would you use from this output?

# alternative for fisher's exact test
fisher.test(df$PRE_female, df$PRE_male)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  df$PRE_female and df$PRE_male
## p-value = 0.6907
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.6553751 1.9758303
## sample estimates:
## odds ratio 
##   1.139965

To report this output, we’ll need the means and sds for each group.

df %>% 
  summarize(mean_female = mean(PRE_female, na.rm = TRUE),
            mean_male = mean(PRE_male, na.rm = TRUE),
            sd_female = sd(PRE_female, na.rm = TRUE),
            sd_male = sd(PRE_male, na.rm = TRUE))
## # A tibble: 1 × 4
##   mean_female mean_male sd_female sd_male
##         <dbl>     <dbl>     <dbl>   <dbl>
## 1       0.388     0.358     0.488   0.480

McNemar Test: Repeated measures chi-square

Visualize the data

Let’s copy the code we used to visually compare PRE_female and PRE_male. How would you change it to show PRE_female and POST_female?

bar_plot_prepost <- df %>%
  dplyr::select(PRE_female, PRE_male) %>% # we only need these variable
  pivot_longer(cols = c(PRE_female, PRE_male), names_to = "category", values_to = "value") %>% 
  mutate(category = factor(category, levels = c("PRE_female", "PRE_male"))) %>% # Set the order of factors
  ggplot(aes(x = category, y = value, fill = category)) + # have to add "fill" in aes() if you want to change the color manually later
  stat_summary(fun=mean, geom="bar", position = position_dodge(), width = .3) + # plot mean as barplot
  stat_summary(fun.data = mean_se, geom = "errorbar", position = position_dodge(0.9), width = .1) + # mean_se is an argument that calculates SD (can calculate SE with additional arguments)
  scale_fill_manual(values = c("cornflowerblue", "deepskyblue4")) + # how to add two colors. they're in the order of the df
  labs(x = "Time", y = "Mean") + # add title to x and y-axis (labs = "labels")
  theme_minimal() + # preset theme
  ggtitle("Means of Pre- and Post-Measures for Women")

print(bar_plot_prepost)
# table (var1 = column, var2 = row). No names if two variables 
table(df$PRE_female, df$POST_female)
##    
##      0  1
##   0 85 74
##   1 24 77
# alternative that shows names
xtabs( ~ PRE_female + POST_female, df)
##           POST_female
## PRE_female  0  1
##          0 85 74
##          1 24 77

How many women are only in the PRE?

__

How many women are only in the POST?

__

How many women do we have in both?

__

mcnemar.test(df$PRE_female, df$POST_female, correct = FALSE)
## 
##  McNemar's Chi-squared test
## 
## data:  df$PRE_female and df$POST_female
## McNemar's chi-squared = 25.51, df = 1, p-value = 4.4e-07
# the continuity correction (default) is a conservative test, best for super small samples and uneven distributions
df %>% 
  summarize(mean_pre = mean(PRE_female, na.rm = TRUE),
            mean_post = mean(POST_female, na.rm = TRUE),
            sd_pre = sd(PRE_female, na.rm = TRUE),
            sd_post = sd(POST_female, na.rm = TRUE))
## # A tibble: 1 × 4
##   mean_pre mean_post sd_pre sd_post
##      <dbl>     <dbl>  <dbl>   <dbl>
## 1    0.388     0.581  0.488   0.494

What if I want an effect size for McNemar or Fisher’s exact test?

We don’t typically use Cohen’s W for these tests, which captures the strength of an association. Instead, we use an odds ratio. Odds ratios work for any categorical test and they’re pretty easy to interpret. Know that it can only compare two variables.

Odds ratio : odds of numerator vs denominator

  • odds < 1 = numerator more likely
  • odds > 1 = denomninator more likely
# Option 2
library(exact2x2)

mcnemar.exact(df$PRE_female, df$POST_female)
## 
##  Exact McNemar test (with central confidence intervals)
## 
## data:  df$PRE_female and df$POST_female
## b = 74, c = 24, p-value = 4.216e-07
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  1.922863 5.110849
## sample estimates:
## odds ratio 
##   3.083333
fisher.test(df$PRE_female, df$POST_female)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  df$PRE_female and df$POST_female
## p-value = 2.815e-06
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.052735 6.713070
## sample estimates:
## odds ratio 
##   3.666142

Power Analysis

Finally, we’re going to do a power analysis using the WebPower r package. To do this for a chi-square, we’re going to run a power analysis for “tests of proportions.” This is also what it’s called in G*Power.

# note that alpha = 0.5 is the default argument. 
# if you're running this post-hoc, include n1 (and n2) 

wp.prop(h = .2, # cohen's h: effect size
        type = "1p", # 1p = 1 sample, 2p = 2 sample, 2p2n = 2 sample, unequal sample size
        alpha = .05, 
        power = .95)
## Power for one-sample proportion test
## 
##       h        n alpha power
##     0.2 324.8677  0.05  0.95
## 
## URL: http://psychstat.org/prop