Learning objectives

By the end of the lecture, you will be able to …

  • Use R to calculate confidence intervals for means and proportions
  • Use R to conduct one and two sample t-tests

Mini-task

Install the Rmisc() package, available on CRAN.


Run the code chunk in your code-along.

OR: Copy and paste the code into your Console pane. Then hit “Enter”.


install.packages("Rmisc")

Mini-task

Load the standard packages and our new Rmisc() package.


library(tidyverse)
library(haven) # not core tidyverse
library(gssr)
library(gssrdoc) # load the GSS codebook
library(summarytools)
library(Rmisc) # load the new package

Mini-task

Load the 2022 GSS data.


# Get the data from the 2022 survey
gss22 <- gss_get_yr(2022)

Variables

  • wwwhr Not counting e-mail, about how many minutes or hours per week do you use the Web?
  • postlife Do you believe there is a life after death?
  • childs How many children have you ever had?
  • sex Respondent’s sex

Mini-task

Create a dataset named my_data with the variables wwwhr, postlife, childs, & sex with no missing data


# Create a new df
my_data <-gss22 |>
  select(wwwhr, postlife, childs, sex) |>
  drop_na()

Confidence Intervals

Review: summarise()

my_data |>
  summarise(
    n = n(),
    mean = mean(wwwhr),
    sd = sd(wwwhr)
    )
# A tibble: 1 × 3
      n  mean    sd
  <int> <dbl> <dbl>
1   970  17.7  22.3

summarise() with equation

wwwhr_sum <- my_data |>
  summarise(
    n = n(),
    mean_wwwhr = mean(wwwhr),
    sd_wwwhr = sd(wwwhr),
    se_wwwhr = sd_wwwhr / sqrt(n),
    ci_lower = mean_wwwhr - (1.96 * se_wwwhr),
    ci_upper = mean_wwwhr + (1.96 * se_wwwhr)
    )

summarise() + CI()

wwwhr_ci <- my_data |>
  summarise(
    n = n(),
    mean_wwwhr = mean(wwwhr),
    sd_wwwhr = sd(wwwhr),
    lowCI_wwwhr = Rmisc::CI(wwwhr)['lower'],
    hiCI_wwwhr  = Rmisc::CI(wwwhr)['upper']
    )

wwwhr_ci
# A tibble: 1 × 5
      n mean_wwwhr sd_wwwhr lowCI_wwwhr hiCI_wwwhr
  <int>      <dbl>    <dbl>       <dbl>      <dbl>
1   970       17.7     22.3        16.3       19.1

Neat! But how confident are we?

summarise() + CI()

wwwhr_ci99 <- my_data |>
  summarise(
    n = n(),
    mean_wwwhr = mean(wwwhr),
    sd_wwwhr = sd(wwwhr),
    lowCI_wwwhr = CI(wwwhr, ci = 0.99)['lower'],
    hiCI_wwwhr  = CI(wwwhr, ci = 0.99)['upper']
    )

wwwhr_ci99
# A tibble: 1 × 5
      n mean_wwwhr sd_wwwhr lowCI_wwwhr hiCI_wwwhr
  <int>      <dbl>    <dbl>       <dbl>      <dbl>
1   970       17.7     22.3        15.8       19.5

CIs & percision

95% confidence intervals

wwwhr_ci
# A tibble: 1 × 5
      n mean_wwwhr sd_wwwhr lowCI_wwwhr hiCI_wwwhr
  <int>      <dbl>    <dbl>       <dbl>      <dbl>
1   970       17.7     22.3        16.3       19.1

99% confidence intervals

wwwhr_ci99
# A tibble: 1 × 5
      n mean_wwwhr sd_wwwhr lowCI_wwwhr hiCI_wwwhr
  <int>      <dbl>    <dbl>       <dbl>      <dbl>
1   970       17.7     22.3        15.8       19.5


How does a 95% confidence interval differ from a 99% confidence interval in terms of width and certainty?

CIs (wrong) with proportions

my_data |>
  summarise(
    n = n(),
    mean_ghosts = mean(postlife),
    sd_ghosts = sd(postlife),
    lowCI_ghosts = CI(postlife, ci = 0.99)['lower'],
    hiCI_ghosts  = CI(postlife, ci = 0.99)['upper']
    )
# A tibble: 1 × 5
      n mean_ghosts sd_ghosts lowCI_ghosts hiCI_ghosts
  <int>       <dbl>     <dbl>        <dbl>       <dbl>
1   970        1.19     0.389         1.15        1.22

This doesn’t work because the values aren’t coded 0 and 1.

CIs (correct) with proportions

# Use logic to recode (yes: 1 = 1; no: 2 = 0)
my_data$ghosts <- ifelse(my_data$postlife == 1, 1, 0)
  
my_data |>
  summarise(
    n = n(),
    mean_ghosts = mean(ghosts),
    sd_ghosts = sd(ghosts),
    lowCI_ghosts = CI(ghosts, ci = 0.99)['lower'],
    hiCI_ghosts  = CI(ghosts, ci = 0.99)['upper']
    )
# A tibble: 1 × 5
      n mean_ghosts sd_ghosts lowCI_ghosts hiCI_ghosts
  <int>       <dbl>     <dbl>        <dbl>       <dbl>
1   970       0.814     0.389        0.782       0.847

Comparing CIs

my_data$sex <- zap_missing(my_data$sex) 
my_data$sex <- as_factor(my_data$sex)   
my_data$sex <- droplevels(my_data$sex)  

df_ghosts <- my_data |>
  group_by(sex) |>
  summarise(
    n = n(),
    mean = mean(ghosts),
    sd = sd(ghosts),
    lower = CI(ghosts, ci = 0.95)['lower'],
    upper = CI(ghosts, ci = 0.95)['upper']
    )
df_ghosts
# A tibble: 2 × 6
  sex        n  mean    sd lower upper
  <fct>  <int> <dbl> <dbl> <dbl> <dbl>
1 male     472 0.773 0.419 0.735 0.811
2 female   498 0.853 0.354 0.822 0.885

Comparing CIs

df_wwwhr <- my_data |>
  group_by(sex) |>
  summarise(
    n = n(),
    mean = mean(wwwhr),
    sd = sd(wwwhr),
    lower = CI(wwwhr, ci = 0.95)['lower'],
    upper  = CI(wwwhr, ci = 0.95)['upper']
    )
df_wwwhr
# A tibble: 2 × 6
  sex        n  mean    sd lower upper
  <fct>  <int> <dbl> <dbl> <dbl> <dbl>
1 male     472  16.2  19.8  14.4  17.9
2 female   498  19.1  24.4  16.9  21.2

Hypothesis Testing

t.test()

Performs one and two sample t-tests. To use this function, we’ll need to specify a few things:

  • variable of interest
  • mu (H0): a number indicating the true value of the mean
  • alternative hypothesis (H1): “two.sided”, “greater”, or “less”

t.test() in action

In your family sociology course you learned that fertility rates have been dropping.

You used to think that people had 2 children, on average. But, now you suspect that people are having fewer than 2 children.


How do you test your hypothesis?

Review: Summary Statistics

my_data |>
  summarise(
    n = n(),
    mean = mean(childs),
    median = median(childs),
    sd = sd(childs))
# A tibble: 1 × 4
      n  mean median    sd
  <int> <dbl>  <dbl> <dbl>
1   970  1.79      2  1.64

Is the difference between our sample mean and 2 statistically significant or could we have gotten this mean simply by chance (sample variation)?

Review: CIs

my_data |>
  summarise(
    n = n(),
    mean = mean(childs),
    median = median(childs),
    sd = sd(childs),
    se = sd / sqrt(n),
    lower = CI(childs, ci = 0.99)['lower'],
    upper  = CI(childs, ci = 0.99)['upper']
    )
# A tibble: 1 × 7
      n  mean median    sd     se lower upper
  <int> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>
1   970  1.79      2  1.64 0.0526  1.66  1.93

Does the confidence interval overlap with 2? What does that tell you?

One sample t.test() for means

t.test(my_data$childs)

    One Sample t-test

data:  my_data$childs
t = 34.102, df = 969, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 1.690587 1.897041
sample estimates:
mean of x 
 1.793814

H0 = 0 by default

One sample t.test() for means

t.test(my_data$childs, mu=2, alternative="less") # Ho: mu=2

    One Sample t-test

data:  my_data$childs
t = -3.9197, df = 969, p-value = 4.744e-05
alternative hypothesis: true mean is less than 2
95 percent confidence interval:
    -Inf 1.88042
sample estimates:
mean of x 
 1.793814

H0 = 2; one-tail test

two sample t.test() in action

Are the proportions of men and women who believe in ghosts statistically significantly different?

t.test(my_data$ghosts ~ my_data$sex, alternative = "two.sided")

    Welch Two Sample t-test

data:  my_data$ghosts by my_data$sex
t = -3.2072, df = 923.29, p-value = 0.001387
alternative hypothesis: true difference in means between group male and group female is not equal to 0
95 percent confidence interval:
 -0.1291289 -0.0310882
sample estimates:
  mean in group male mean in group female 
           0.7733051            0.8534137

Aside: paired t.test()

Independent group t-test
Do two subsamples have the same mean?

t.test(my_data$ghosts ~ my_data$sex, alternative = "two.sided")


Paired t-test
Do two variables have the same mean?

t.test(Pair(hrs1_2018, hrs1_2020) ~ 1, data = gss_panel)

Joining Dataframes (if time allows…)

appending v.s. merging

APPEND

add new observations (rows) to existing variables


MERGE

add new variables (columns) to existing observations (many merge types)

appending v.s. merging

Example datasets

dataframe 1

  coupleid  name age
1        2  John  42
2        1 Megan  36
3        3   Bin  38

dataframe 2

  coupleid    name age
1        1     Sue  40
2        3  Ye-jin  39
3        2 Chrissy  35


dataframe 3

  coupleid marstat numchild country
1        3       1        1 S.Korea
2        1       0        0      US
3        4       1        2  Canada

append data with bind_rows()

df_all <- bind_rows(df_partner1, df_partner2)

tibble(df_all)
# A tibble: 6 × 3
  coupleid name      age
     <dbl> <chr>   <dbl>
1        2 John       42
2        1 Megan      36
3        3 Bin        38
4        1 Sue        40
5        3 Ye-jin     39
6        2 Chrissy    35

merge data with joins

add columns from df1 to df2, matching observations based on the keys

  • left_join() keeps all observations in df1.
  • right_join() keeps all observations in df2.
  • full_join() keeps all observations in df1 and df2.
  • inner_join() only keeps observations from df1 that have a matching key in df2

left_join()

df_couples <- left_join(df_partner1, df_family, by = "coupleid")

tibble(df_couples)
# A tibble: 3 × 6
  coupleid name    age marstat numchild country
     <dbl> <chr> <dbl>   <dbl>    <dbl> <chr>  
1        2 John     42      NA       NA <NA>   
2        1 Megan    36       0        0 US     
3        3 Bin      38       1        1 S.Korea

right_join()

df_couples <- right_join(df_partner1, df_family, by = "coupleid")

tibble(df_couples)
# A tibble: 3 × 6
  coupleid name    age marstat numchild country
     <dbl> <chr> <dbl>   <dbl>    <dbl> <chr>  
1        1 Megan    36       0        0 US     
2        3 Bin      38       1        1 S.Korea
3        4 <NA>     NA       1        2 Canada

full_join()

df_couples <- full_join(df_partner1, df_family, by = "coupleid")

tibble(df_couples)
# A tibble: 4 × 6
  coupleid name    age marstat numchild country
     <dbl> <chr> <dbl>   <dbl>    <dbl> <chr>  
1        2 John     42      NA       NA <NA>   
2        1 Megan    36       0        0 US     
3        3 Bin      38       1        1 S.Korea
4        4 <NA>     NA       1        2 Canada

inner_join()

df_couples <- inner_join(df_partner1, df_family, by = "coupleid")

tibble(df_couples)
# A tibble: 2 × 6
  coupleid name    age marstat numchild country
     <dbl> <chr> <dbl>   <dbl>    <dbl> <chr>  
1        1 Megan    36       0        0 US     
2        3 Bin      38       1        1 S.Korea