POLS 1140

Measure	Mean
Unconditional
Overall	56.88
Conditional
No college degree	53.33
College degree	61.03
Under 30	61.95
Over Thirty	55.68

Measure	Mean
Unconditional
Overall	56.88
Conditional
No college degree	53.33
College degree	61.03
Under 30	61.95
Over Thirty	55.68

term	estimate	std.error	statistic	p.value	conf.low	conf.high	df	outcome
(Intercept)	61.95	1.33	46.50	0	59.34	64.57	1691	ft_professors
age_catOver Thirty	-6.27	1.54	-4.07	0	-9.30	-3.25	1691	ft_professors

Statistical models
	Model 1	Model 2	Model 3	Model 4
(Intercept)	61.95^***	67.26^***	53.33^***	64.35^***
	(1.33)	(1.81)	(0.94)	(1.83)
age_catOver Thirty	-6.27^***
	(1.54)
age		-0.21^***		-0.23^***
		(0.04)		(0.04)
has_degreeCollege degree			7.70^***	8.36^***
			(1.34)	(1.34)
R²	0.01	0.02	0.02	0.04
Adj. R²	0.01	0.02	0.02	0.04
Num. obs.	1693	1693	1693	1693
RMSE	27.80	27.65	27.64	27.35
^*p < 0.001; ^p < 0.01; ^*p < 0.05

Inference and Uncertainty

Statistical inference involves quantifying uncertainty about what could have happened.

Today, we’ll introduce the concepts of:

Samples and Populations
Confidence Intervals and Hypotheses Tests

There is more content here than we’ll discuss in class.

You don’t need to know how to conduct a hypothesis test or construct a confidence interval
You do need a functional understanding about how to use these tools to understand claims about statistical significance

Sampling distributions

When we conduct a survey we are trying to learn about a population by generalizing from a specific sample
What would have happened if we had a different sample? How different would our result be?
Let’s treat the 2024 NES pilot as the population
Take repeated samples of size N = 10, 30, 300
For each sample of size N, calculate the sample mean of feelings torward professors
Plot the distribution of sample means (i.e. the sampling distribution)

# Load Data
# load(url("https://pols1600.paultesta.org/files/data/nes24.rda"))

# ---- Population ----

# Population average
mu_prof <- mean(df$ft_professors, na.rm=T)
# Population standard deviation
sd_prof <- sd(df$ft_professors, na.rm = T)

# ---- Function to Take Repeated Samples From Data ----

sample_data_fn <- function(
    dat=df, var=ft_professors, samps=1000, sample_size=10,
    resample = F){
  if(resample == F){
  df <- tibble(
  sim = 1:samps,
  distribution = "Sampling",
  size = sample_size,
  sample_from = "Population",
  pop_mean = dat %>% pull(!!enquo(var)) %>% mean(., na.rm=T),
  pop_sd = dat %>% pull(!!enquo(var)) %>% sd(., na.rm=T),
  se_asymp = pop_sd / sqrt(size),
  ll_asymp = pop_mean - 1.96*se_asymp,
  ul_asymp = pop_mean + 1.96*se_asymp,
) %>% 
  mutate(
    sample = purrr::map(sim, ~ slice_sample(dat %>% select(!!enquo(var)), n = sample_size, replace = F)),
    sample_mean = purrr::map_dbl(sample, \(x) x %>% pull(!!enquo(var)) %>% mean(.,na.rm=T)),
    ll = sample_mean - 1.96*sd(sample_mean),
    ul = sample_mean + 1.96*sd(sample_mean)
  )
  }
  if(resample == T){
    df <- tibble(
  sim = 1:samps,
  distribution = "Resampling",
  size = sample_size,
  sample_from = "Sample",
  pop_mean = dat %>% pull(!!enquo(var)) %>% mean(., na.rm=T),
  pop_sd = dat %>% pull(!!enquo(var)) %>% sd(., na.rm=T),
  se_asymp = pop_sd / sqrt(size),
  ll_asymp = pop_mean - 1.96*se_asymp,
  ul_asymp = pop_mean + 1.96*se_asymp,
) %>% 
  mutate(
    sample = purrr::map(sim, ~ slice_sample(dat %>% select(!!enquo(var)), n = sample_size, replace = T)),
    sample_mean = purrr::map_dbl(sample, \(x) x %>% pull(!!enquo(var)) %>% mean(.,na.rm=T))
  )
  }
  return(df)
}

# ---- Plot Single Distribution -----

plot_distribution <- function(the_pop,the_samp, the_var, ...){
  mu_pop <- the_pop %>% pull(!!enquo(the_var)) %>% mean(., na.rm=T)
  mu_samp <- the_samp %>% pull(!!enquo(the_var)) %>% mean(., na.rm=T)
  ll <- the_pop %>% pull(!!enquo(the_var)) %>% as.numeric() %>%  min(., na.rm=T)
  ul <- the_pop %>% pull(!!enquo(the_var)) %>% as.numeric() %>% max(., na.rm=T)
  p<- the_samp %>% 
    ggplot(aes(!!enquo(the_var)))+
    geom_density()+
    geom_rug()+
    theme_void()+
    geom_vline(xintercept = mu_samp, col = "red")+
    geom_vline(xintercept = mu_pop, col = "grey40",linetype = "dashed")+
    xlim(ll,ul)
  return(p)
}

# ---- Plot multiple distributions ----

plot_samples <- function(pop, x, variable,n_rows = 4, ...){
  sample_plots <- x$sample[1:(4*n_rows)] %>% 
  purrr::map( \(x) plot_distribution(the_pop=pop, the_samp = x, 
                                     the_var = !!enquo(variable)))
  p <- wrap_elements(wrap_plots(sample_plots[1:(4*n_rows)], ncol=4))
  return(p)
  
}

# ---- Plot Combined Figure ----

plot_figure_fn <- function(
    d=df, 
    v=age, 
    sim=1000, 
    size=10,
    rows = 4){
  # Population average
  mu <- d %>% pull(!!enquo(v)) %>% mean(., na.rm=T)
  sd <- d %>% pull(!!enquo(v)) %>% sd(., na.rm=T)
  se <- sd/sqrt(size)
  # Range
  ll <- d %>% pull(!!enquo(v)) %>% as.numeric() %>%  min(., na.rm=T)
  ul <- d %>% pull(!!enquo(v)) %>% as.numeric() %>% max(., na.rm=T)
  # Population standard deviation
  # Sample data
  samp_df <- sample_data_fn(dat=d, var = !!enquo(v), samps = sim, sample_size = size)
  # Plot Population
  p_pop <- d %>%
    ggplot(aes(!!enquo(v)))+
      geom_density(col ="grey60")+
      geom_rug(col = "grey60", )+
      geom_vline(xintercept = mu, col="grey40", linetype="dashed")+
      theme_void()+
      labs(title ="Population")+
      xlim(ll,ul)+
      theme(plot.title = element_text(hjust = 0))

  
  p_samps <- plot_samples(pop=d, x= samp_df,variable = !!enquo(v),
                          n_rows = rows)
  p_samps <- p_samps + 
    ggtitle(paste("Repeated samples of size N =",size,"from the population"))+
    theme(plot.title = element_text(hjust = 0.5), 
          plot.background = element_rect(
            fill = NA, colour = 'black', linewidth = 2)
          )
  
  
  p_dist <- samp_df %>% 
  ggplot(aes(sample_mean))+
  geom_density(col="red",aes(y= after_stat(ndensity)))+
  geom_rug(col="red")+
  geom_density(data = df, aes(!!enquo(v), y= after_stat(ndensity)),
               col="grey60")+
  geom_vline(xintercept = mu, col="grey40", linetype="dashed")+
  xlim(ll,ul)+
  theme_void()+
    labs(
      title = "Sampling Distribution"
    )+  theme(plot.title = element_text(hjust = 0))
  
  range_upper_df <- tibble(
  x = seq( ((ll+ul)/2 -5), ((ll+ul)/2 +5), length.out = 20),
  xend = seq(ll-5, ul+5, length.out = 20),
  y = rep(9, 20),
  yend = rep(1, 20)
)
p_upper <- range_upper_df %>% 
  ggplot(aes(x=x, xend = xend, y=y,yend=yend))+
  geom_segment(
    arrow = arrow(length = unit(0.05, "npc"))
  )+
  theme_void()+
  coord_fixed(ylim=c(0,10),
              xlim =c(ll-5,ul+5),clip="off")
  # Lower
  range_df <- samp_df %>% 
  summarise(
    min = min(sample_mean),
    max = max(sample_mean),
    mean = mean(sample_mean)
  )
  
  plot_df <- tibble(
  id = 1:50,
  # x = sort(rnorm(50, mu, sd)),
  x = sort(runif(50, ll, ul)),
  xend = sort(rnorm(50, mu, se)),
  y = 9,
  yend = 1
)

p_lower <- plot_df %>%
  ggplot(aes(x,y, group =id))+
  geom_segment(aes(xend=xend, yend=yend),
               col = "red",arrow = arrow(length = unit(0.05, "npc"))
               )+
  theme_void()+
  coord_fixed(ylim=c(0,10),xlim = c(ll,ul),clip="off")

  
  design <-"##AAAA##
            ##AAAA##
            ##AAAA##
            BBBBBBBB
            BBBBBBBB
            #CCCCCC#
            #CCCCCC#
            #CCCCCC#
            #CCCCCC#
            DDDDDDDD
            DDDDDDDD
            ##EEEE##
            ##EEEE##
            ##EEEE##"
  
  fig <- p_pop / p_upper / p_samps / p_lower / p_dist +
    plot_layout(design = design)
  return(fig)


  
  
  
}

# ---- Samples and Figures Varying Sample Size ----
## N = 10
set.seed(1234)
samp_n10 <- sample_data_fn(sample_size  = 10, samps = 1000)
set.seed(1234)
fig_n10 <- plot_figure_fn(v=ft_professors,size = 10)

## N = 30
set.seed(1234)
samp_n30 <- sample_data_fn(sample_size  = 30, samps = 1000)
set.seed(1234)
fig_n30 <- plot_figure_fn(v=ft_professors,size = 30,rows=4)

## N = 300
set.seed(1234)
samp_n300 <- sample_data_fn(sample_size  = 300, samps = 1000)
set.seed(1234)
fig_n300 <- plot_figure_fn(v=ft_professors,size = 300)

Standard errors

The standard error (SE) is simply the standard deviation of the sampling distribution.
The SE decreases as the sample size increases (by the LLN):
Approximately 95% of the sample means will be within 2 SEs of the population mean (CLT)

se_df <- tibble(
  `Sample Size` = factor(paste("N =",c(10,30, 300))),
  se = c(sd(samp_n10$sample_mean),
         sd(samp_n30$sample_mean),
         sd(samp_n300$sample_mean)),
  SE = paste("SE =", round(se,2)),
  ll = mu_prof,
  ul = mu_prof + se,
  y = c(.3,.3,.45),
  yend = y
)

ci_df <- tibble(
  `Sample Size` = factor(paste("N =",c(10,30, 300))),
  se = c(sd(samp_n10$sample_mean),
         sd(samp_n30$sample_mean),
         sd(samp_n300$sample_mean)),
  mu = mu_prof,
  ll = round(mu_prof - 1.96 *se,2),
  ul = round(mu_prof + 1.96 *se,2),
  ci = paste("95 % Coverage Interval [",ll,";",ul,"]",sep=""),
  y = c(.3,.3,.45),
  yend = y
)
sim_df <- samp_n10 %>% 
  bind_rows(samp_n30) %>% 
  bind_rows(samp_n300) %>% 
  mutate(
    `Sample Size` = factor(paste("N =",size))
    ) %>% 
  left_join(ci_df) %>% 
  mutate(
    Coverage = case_when(
      sample_mean > ll_asymp & sample_mean < ul_asymp  & size == 10~ "#F8766D",
      sample_mean > ll_asymp & sample_mean < ul_asymp  & size == 30~ "#00BA38",
      sample_mean > ll_asymp & sample_mean < ul_asymp  & size == 300~ "#619CFF",
      T ~ "grey"
    )
  )



fig_se <- sim_df %>% 
  ggplot(aes(sample_mean, col = `Sample Size`))+
  geom_density()+
  geom_rug()+
  geom_vline(xintercept = mu_prof, linetype = "dashed")+
  theme_minimal()+
  facet_wrap(~`Sample Size`, ncol=1)+
  ylim(0,.5)+
  guides(col="none")+
  geom_segment(
    data = se_df,
    aes(x= ll, xend =ul, y = y, yend = yend)
  )+
  geom_text(
    data = se_df,
    aes(x = ul, y =y, label = SE),
    hjust = -.25
  ) +
  labs(
    y = "",
    x = "Sampling Distributions of Sample Means",
    title = "Standard Errors decrease with Sample Size"
  )

fig_coverage <- sim_df %>% 
  ggplot(aes(sample_mean,col=`Sample Size`))+
  geom_density()+
  geom_rug(col=sim_df$Coverage)+
  geom_vline(xintercept = mu_prof, linetype = "dashed")+
  theme_minimal()+
  facet_wrap(~`Sample Size`, ncol=1)+
  ylim(0,.55)+
  guides(col="none")+
  geom_segment(
    data = ci_df,
    aes(x= ll, xend =ul, y = y, yend = yend)
  )+
  geom_text(
    data = ci_df,
    aes(x = mu, y =y, label = ci),
    hjust = .5,
    nudge_y =.1
  ) +
  labs(
    y = "",
    x = "Sampling Distributions of Sample Means",
    title = "Approximately 95% of sample means are within 2 SE of the population mean"
  )

How do we calculate a standard error from a single sample?

Calculating standard errors

Simulation:
- Treat sample as population
- Sample with replacement (“bootstrapping”)
- Estimate SE from standard deviation of resampling distribution (“plug-in principle”)
Analytic
- Characterize sampling distribution from sample mean and variance via asymptotic theory (the LLT and CLT)
- For a sample mean, $\bar{x}$

$S E_{\bar{x}} = \frac{σ_{x}}{\sqrt{(} n)}$

plot_resampling_fn <- function(d=df, v=age, sim=1000, size=10,rows=3){
  # Population average
  mu <- d %>% pull(!!enquo(v)) %>% mean(., na.rm=T)
  # Population standard deviation and SE
  sd <- d %>% pull(!!enquo(v)) %>% sd(., na.rm=T)
  se <- sd/sqrt(size)
  # Range
  ll <- d %>% pull(!!enquo(v)) %>% as.numeric() %>%  min(., na.rm=T)
  ul <- d %>% pull(!!enquo(v)) %>% as.numeric() %>% max(., na.rm=T)
  # Resampling with replace
  # Draw 1 Sample
  sample <- sample_data_fn(dat=d, var = !!enquo(v), samps = 1, sample_size = size, resample = F)
  samp_df <- as.data.frame(sample$sample)
  # Resample from sample with replacement
  resamp_df <- sample_data_fn(dat=samp_df, var = !!enquo(v), samps = sim, sample_size = size, resample = T)
  # Plot Population
  p_pop <- d %>%
    ggplot(aes(!!enquo(v)))+
      geom_density(col ="grey60")+
      geom_rug(col = "grey60", )+
      geom_vline(xintercept = mu, col="grey40", linetype="dashed")+
      theme_void()+
      labs(title ="Population")+
      xlim(ll,ul)+
      theme(plot.title = element_text(hjust = 0))

  p_samp <- plot_distribution(the_pop = d,
                              the_samp = samp_df,
                              the_var = age)+
    labs(title ="Sample")+
      xlim(ll,ul)+
      theme(plot.title = element_text(hjust = 0))
  
  p_samps <- plot_samples(pop=d, x= resamp_df,variable = !!enquo(v), n_rows =rows)
  p_samps <- p_samps + 
    ggtitle(paste("Repeated samples with replacement\nof size N =",size,"from sample"))+
    theme(plot.title = element_text(hjust = 0.5), 
          plot.background = element_rect(
            fill = NA, colour = 'black', linewidth = 2)
          )
  
  # Resampling Distribution
  
  
  p_dist <- resamp_df %>% 
  ggplot(aes(sample_mean))+
  geom_density(col="red",aes(y= after_stat(ndensity)))+
  geom_rug(col="red")+
  geom_density(data = df, aes(!!enquo(v), y= after_stat(ndensity)),
               col="grey60")+
  geom_vline(xintercept = unique(resamp_df$pop_mean), col="red", linetype="solid")+
  geom_vline(xintercept = mu, col="grey40", linetype="dashed")+
  xlim(ll,ul)+
  theme_void()+
    labs(
      title = "Reampling Distribution"
    )+  theme(plot.title = element_text(hjust = 0))
  
   range_upper_df <- tibble(
  x = seq( ((ll+ul)/2 -5), ((ll+ul)/2 +5), length.out = 20),
  xend = seq(ll-5, ul+5, length.out = 20),
  y = rep(9, 20),
  yend = rep(1, 20)
)
p_upper <- range_upper_df %>% 
  ggplot(aes(x=x, xend = xend, y=y,yend=yend))+
  geom_segment(
    arrow = arrow(length = unit(0.05, "npc"))
  )+
  theme_void()+
  coord_fixed(ylim=c(0,10),
              xlim =c(ll-5,ul+5),clip="off")
  # Lower
  range_df <- resamp_df %>% 
  summarise(
    min = min(sample_mean),
    max = max(sample_mean),
    mean = mean(sample_mean)
  )
  
  plot_df <- tibble(
  id = 1:50,
  # x = sort(rnorm(50, mu, sd)),
  x = sort(runif(50, ll, ul)),
  xend = sort(rnorm(50, unique(resamp_df$pop_mean), se)),
  y = 9,
  yend = 1
)

p_lower <- plot_df %>%
  ggplot(aes(x,y, group =id))+
  geom_segment(aes(xend=xend, yend=yend),
               col = "red",arrow = arrow(length = unit(0.05, "npc"))
               )+
  theme_void()+
  coord_fixed(ylim=c(0,10),xlim = c(ll,ul),clip="off")

  
  design <-"##AAAA##
            ##AAAA##
            ##AAAA##
            ##BBBB##
            ##BBBB##
            ##BBBB##            
            CCCCCCCC
            CCCCCCCC
            #DDDDDD#
            #DDDDDD#
            #DDDDDD#
            #DDDDDD#
            EEEEEEEE
            EEEEEEEE
            ##FFFF##
            ##FFFF##
            ##FFFF##"
  
  fig <- p_pop / p_samp /p_upper / p_samps / p_lower / p_dist +
    plot_layout(design = design)
  return(fig)


  
  
  
}
set.seed(123)
resamp_n10 <- sample_data_fn(
  dat = sample_data_fn(samps = 1, sample_size = 10, resample = T)$sample %>%  as.data.frame(),
  sample_size = 10, 
  resample = T)
set.seed(123)
fig_n10_bs <- plot_resampling_fn(size=10)

set.seed(12345)
resamp_n30 <- sample_data_fn(
  dat = sample_data_fn(samps = 1, sample_size = 30, resample = T)$sample %>%  as.data.frame(),
  samps = 1000, sample_size = 30, resample = T)

set.seed(12345)
fig_n30_bs <- plot_resampling_fn(size=30)

set.seed(1234)
resamp_n300 <- sample_data_fn(
  dat = sample_data_fn(samps = 1, sample_size = 300, resample = T)$sample %>%  as.data.frame(),
  samps = 1000, sample_size = 300, resample = T)
set.seed(1234)
fig_n300_bs <- plot_resampling_fn(size=300)

Bootstrap SE	Analytic SE
9.85	8.82
5.79	5.09
1.85	1.61

Confidence intervals

Confidence intervals:

provide a way of quantifying uncertainty about estimates
describe a range of plausible values for an estimate
are a function of the standard error of the estimate, and the a critical value determined by $α$ , which describes the degree of confidence we want

Calculating a confidence interval

Choose level of confidence $(1 - α) \times 100$
- $α = 0.05$ , corresponds to a 95% confidence level.
Derive the sampling distribution of the estimator
- Simulation: bootstrap re-sampling
- Analytically: computing its mean and variance.
Compute the standard error
Compute the critical value $z_{α / 2}$
- as the $1.96 = Φ (z_{0.5 / 2})$ for a 95% CI
Compute the lower and upper confidence limits
- lower limit = $\hat{θ} - z_{α / 2} \times S E$
- upper limit = $\hat{θ} + z_{α / 2} \times S E$

resamp_df <- 
  resamp_n10 %>% 
  bind_rows(resamp_n30) %>% 
  bind_rows(resamp_n300) %>% 
  mutate(
    `Sample Size` = factor(paste("N =",size))
    )

resamp_ci_df <- tibble(
  `Sample Size` = factor(paste("N =",c(10,30,300))),
  mu = unique(resamp_df$pop_mean),
  ll = unique(resamp_df$ll_asymp),
  ul = unique(resamp_df$ul_asymp),
  y = c(.3, .3,.5)
)

fig_ci1 <- resamp_df %>% 
  ggplot(aes(sample_mean,
             col = `Sample Size`))+
  geom_density()+
  geom_rug()+
  geom_vline(xintercept = mu_prof, linetype = "dashed")+
  geom_vline(data = resamp_ci_df,
             aes(xintercept = mu,
                 col = `Sample Size`))+
  geom_segment(data = resamp_ci_df,
               aes(x = ll, xend =ul, y = y, yend =y,
                   col = `Sample Size`))+
  facet_wrap(~`Sample Size`, ncol=1)+
  theme_minimal()+
  labs(
    y = "",
    x = "Resampling Distribution",
    title = "95% Confidence Intervals"
  )
  

samp_ci_df <- samp_n10 %>% 
  bind_rows(samp_n30) %>% 
  bind_rows(samp_n300) %>% 
  mutate(
    `Sample Size` = factor(paste("N =",size))
    ) %>% 
  mutate(
    Coverage = case_when(
      pop_mean > ll & pop_mean < ul ~ "red",
      T ~ "black"
    )
  )

fig_ci2 <- samp_ci_df %>% 
  filter(sim %in% 1:100) %>% 
  filter(size == 10) %>% 
  ggplot(aes(y = sample_mean, x= sim))+
  geom_pointrange(aes(ymin = ll, ymax =ul, col=Coverage))+
  geom_hline(yintercept = mu_prof, linetype = "dashed")+
  coord_flip()+
  theme_minimal()+
  guides(col = "none")+
  facet_wrap(~`Sample Size`)

fig_ci3 <- samp_ci_df %>% 
  filter(sim %in% 1:100) %>% 
  ggplot(aes(y = sample_mean, x= sim))+
  geom_pointrange(aes(ymin = ll, ymax =ul, col=Coverage))+
  geom_hline(yintercept = mu_prof, linetype = "dashed")+
  coord_flip()+
  theme_minimal()+
  guides(col = "none")+
  facet_wrap(~`Sample Size`)

Interpreting confidence intervals

Confidence intervals give a range of values that are likely to include the true value of the parameter $θ$ with probability $(1 - α) \times 100 %$
- $α = 0.05$ corresponds to a “95-percent confidence interval”
Our “confidence” is about the interval
In repeated sampling, we expect that $(1 - α) \times 100 %$ of the intervals we construct would contain the truth.
For any one interval, the truth, $θ$ , either falls within in the lower and upper bounds of the interval or it does not.

Hypothesis testing

What is a hypothesis test

A formal way of assessing statistical evidence. Combines
- Deductive reasoning distribution of a test statistic, if the a null hypothesis were true
- Inductive reasoning based on the test statistic we observed, how likely is it that we would observe it if the null were true?

What is a test statistic?

A way of summarizing data
- difference of means
- coefficients from a linear model
- coefficients from a linear model divided by their standard errors
- R^2
- Sums of ranks

Note

Different test statistics may be more or less appropriate depending on your data and questions.

What is a null hypothesis?

A statement about the world
- Only interesting if we reject it
- Would yield a distribution of test statistics under the null
- Typically something like “X has no effect on Y” (Null = no effect)
- Never accept the null can only reject

What is a p-value?

A p-value is a conditional probability summarizing the likelihood of observing a test statistic as far from our hypothesis or farther, if our hypothesis were true.

How do we do hypothesis testing?

Posit a hypothesis (e.g. $β = 0$ )
Calculate the test statistic (e.g. $(\hat{β} - β) / s e_{β}$ )
Derive the distribution of the test statistic under the null via simulation or asymptotic theory
Compare the test statistic to the distribution under the null
Calculate p-value (Two Sided vs One sided tests)
Reject or fail to reject/retain our hypothesis based on some threshold of statistical significance (e.g. p < 0.05)

Outcomes of hypothesis tests

Two conclusions from of a hypothesis test: we can reject or fail to reject a hypothesis test.
We never “accept” a hypothesis, since there are, in theory, an infinite number of other hypotheses we could have tested.

Our decision can produce four outcomes and two types of error:

	Reject $H_{0}$	Fail to Reject $H_{0}$
$H_{0}$ is true	False Positive	Correct!
$H_{0}$ is false	Correct!	False Negative

Type 1 Errors: False Positive Rate (p < 0.05)
Type 2 Errors: False negative rate (1 - Power of test)

Quantifying uncertainty in regression

How education condition the relationship between age and feelings toward professors?

Let’s fit the following “interaction” model to assess whether education moderates the relationship between age and evaluations

$y = β_{0} + β_{1} age + β_{2} degree + β_{3} age \times degree + ϵ$

m1 <- lm_robust(ft_professors ~ age*has_degree, df)

And unpack the output

tidy(m1) %>% 
  mutate_if(is.numeric, \(x) round(x, 3)) -> m1_sum
m1_sum

                          term estimate std.error statistic p.value conf.low
1                  (Intercept)   62.812     2.326    27.001   0.000   58.249
2                          age   -0.197     0.048    -4.069   0.000   -0.292
3     has_degreeCollege degree   12.168     3.602     3.378   0.001    5.103
4 age:has_degreeCollege degree   -0.076     0.072    -1.064   0.288   -0.217
  conf.high   df       outcome
1    67.375 1689 ft_professors
2    -0.102 1689 ft_professors
3    19.234 1689 ft_professors
4     0.064 1689 ft_professors

htmlreg(m1,include.ci=F)

Statistical models
	Model 1
(Intercept)	62.81^***
	(2.33)
age	-0.20^***
	(0.05)
has_degreeCollege degree	12.17^***
	(3.60)
age:has_degreeCollege degree	-0.08
	(0.07)
R²	0.04
Adj. R²	0.04
Num. obs.	1693
RMSE	27.35
^*p < 0.001; ^p < 0.01; ^*p < 0.05

htmlreg(m1,include.ci=T)

Statistical models
	Model 1
(Intercept)	62.81^*
	[58.25; 67.37]
age	-0.20^*
	[-0.29; -0.10]
has_degreeCollege degree	12.17^*
	[ 5.10; 19.23]
age:has_degreeCollege degree	-0.08
	[-0.22; 0.06]
R²	0.04
Adj. R²	0.04
Num. obs.	1693
RMSE	27.35
^* 0 outside the confidence interval.

pred_df <- expand_grid(
  has_degree = c("College degree", "No college degree"),
  age = 18:80
  )

pred_df <- cbind(pred_df,
    predict(m1, newdata = pred_df, interval = "confidence")$fit
  )

m1_plot <- pred_df %>% 
  ggplot(aes(age, fit, fill = has_degree))+
  geom_ribbon(aes(ymin=lwr, ymax=upr),alpha=0.5)+
  geom_line() +
  labs(
    fill = "Education",
    x = "Age",
    y = "Predicted Feelings toward Professors"
  )+
  theme_minimal()

$x_{i}$	$d_{i}$	$y_{i}$
chocolate	1	7
chocolate	1	8
chocolate	1	5
chocolate	1	4
fruit	0	10
fruit	0	9
chocolate	1	5
fruit	0	8
chocolate	1	4
chocolate	1	6

$x_{i}$	$d_{i}$	$y_{i}$
chocolate	1	7
chocolate	1	8
chocolate	1	5
chocolate	1	4
fruit	0	10
fruit	0	9
chocolate	1	5
fruit	0	8
chocolate	1	4
chocolate	1	6

$x_{i}$	$d_{i}$	$y_{i}$
chocolate	1	7
chocolate	1	8
chocolate	0	4
chocolate	1	4
fruit	0	10
fruit	1	8
chocolate	0	4
fruit	0	8
chocolate	1	4
chocolate	0	0

$x_{i}$	$d_{i}$	$y_{i}$
chocolate	1	7
chocolate	1	8
chocolate	0	4
chocolate	1	4
fruit	0	10
fruit	1	8
chocolate	0	4
fruit	0	8
chocolate	1	4
chocolate	0	0

$x_{i}$	$d_{i}$	$y_{i}$
chocolate	1	7
chocolate	1	8
chocolate	0	4
chocolate	1	4
fruit	0	10
fruit	1	8
chocolate	0	4
fruit	0	8
chocolate	1	4
chocolate	0	0

$Y_{i} (1)$	$Y_{i} (0)$	$τ_{i}$
7	3	4
8	6	2
5	4	1
4	3	1
6	10	-4
8	9	-1
5	4	1
7	8	-1
4	3	1
6	0	6

$Y_{i} (1)$	$Y_{i} (0)$	$τ_{i}$
7	3	4
8	6	2
5	4	1
4	3	1
6	10	-4
8	9	-1
5	4	1
7	8	-1
4	3	1
6	0	6

$Y_{i} (1)$	$Y_{i} (0)$	$τ_{i}$
7	3	4
8	6	2
5	4	1
4	3	1
6	10	-4
8	9	-1
5	4	1
7	8	-1
4	3	1
6	0	6