Exercise 06 Solution

• Solution

Preliminaries

Load in dataset and libraries of interest…

library(tidyverse)
library(mosaic)
library(infer)
f <- "https://raw.githubusercontent.com/difiore/ada-datasets/main/tbs-2006-2008-ranges.csv"
d <- read_csv(f, col_names = TRUE)
head(d)
## # A tibble: 6 × 9
##   id       sex   sex.code mcp50 mcp80 mcp95 kernel50 kernel80 kernel95
##   <chr>    <chr>    <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>    <dbl>
## 1 Ana      F            1 123.   182.  331.     75.1     185.     330.
## 2 Andreo   M            2 136.   238.  367.    101.      232.     416.
## 3 Buka     F            1  98.2  202.  234.     60.5     161.     296.
## 4 Eva      F            1 104.   193.  352.     82.2     202.     355.
## 5 Evita    F            1 104.   183.  365.     82.6     190.     367.
## 6 Geronimo M            2 116.   286.  402.    125.      285.     502.

Challenge

Step 1

  • Reduce the dataset to just the variables of interest.
d <- d |> select(sex, kernel95)
head(d)
## # A tibble: 6 × 2
##   sex   kernel95
##   <chr>    <dbl>
## 1 F         330.
## 2 M         416.
## 3 F         296.
## 4 F         355.
## 5 F         367.
## 6 M         502.

Step 2

  • Determine the mean, standard deviation, and standard error in “kernel95” home range size for each sex.
hr_summary <- d |>
  group_by(sex) |>
  summarize(
    mean = mean(kernel95),
    sd = sd(kernel95),
    n = n(),
    se = sd / sqrt(n)
  )
hr_summary
## # A tibble: 2 × 5
##   sex    mean    sd     n    se
##   <chr> <dbl> <dbl> <int> <dbl>
## 1 F      319.  65.7    11  19.8
## 2 M      430.  58.3     9  19.4

Step 3

  • Create boxplots comparing “kernel95” home range size by sex.
p <- ggplot(data = d, aes(x = sex, y = kernel95)) +
  geom_boxplot() +
  geom_jitter()
p

Step 4

  • For each sex, generate a bootstrap distribution for mean kernel95 home range size.

NOTE: The code below does this for males… simply replace “M” with “F” in the filter() statement to do the same for females.

# ci in mean HR size
n_boot <- 10000
s <- d |> filter(sex == "M")

# option 1 - using {base} R
boot <- vector()
# the size of each bootstrap sample should equivalent to the size our original sample
for (i in 1:n_boot) {
  boot[[i]] <-
    mean(sample(s$kernel95, length(s$kernel95), replace = TRUE))
}

# option 2 - using {mosaic}
# note that the entire do()* loop code needs to be wrapped in braces or parentheses before pull()
boot <-
  {
    do(n_boot) *
      mean(sample(s$kernel95, size = length(s$kernel95), replace = TRUE))
  } |>
  pull(mean) # pull out mean column as vector

# option 3 - using {infer}
boot <- s %>%
  rep_sample_n(
    replace = TRUE,
    size = nrow(.),
    reps = 10000
  ) |>
  group_by(replicate) |>
  summarize(mean = mean(kernel95)) |>
  pull(mean) # pull out mean column as vector

Step 5

  • Plot the resulting bootstrap sampling distribution and plot an appropriate normal distribution over the bootstrap sampling distribution.
se <- sd(boot)
histogram(boot)
plotDist("norm", mean(boot), se, add = TRUE)

Step 6

  • Calculate a 95% confidence interval around for the mean kernel95 home range size for each sex…

    • Using the quantile() method applied directly to your bootstrap sampling distribution, and

    • Using the theory-based “standard error” method, based on qnorm() and the standard deviation of your bootstrap sampling distribution.

(ci_boot <- c(
  quantile(boot, 0.025),
  quantile(boot, 0.975)
))
##     2.5%    97.5% 
## 393.5361 466.1644
(ci_theory <- mean(s$kernel95) + c(-1, 1) * qnorm(0.975) * se)
## [1] 393.6887 465.7927
# or

(ci_theory <- mean(s$kernel95) + qnorm(c(0.025, 0.975)) * se)
## [1] 393.6887 465.7927
ladd(
  panel.abline(v = ci_boot, col = "red", lty = 3, lwd = 2)
)

ladd(
  panel.abline(v = ci_theory, col = "blue", lty = 1, lwd = 2)
)

We can also do all of the above using the {infer} package’s specify() → generate() → calculate() → visualize() workflow…

# option 4 - using {infer}
boot <- s |>
  specify(response = kernel95) |>
  generate(reps = n_boot, type = "bootstrap") |>
  calculate(stat = "mean")

ci_boot <- boot |>
  get_confidence_interval(
    type = "percentile", level = 0.95
  )

ci_theory <- boot |>
  get_confidence_interval(
    type = "se", level = 0.95,
    point_estimate = mean(s$kernel95)
  )

visualize(boot) + shade_confidence_interval(endpoints = ci_theory, color = "blue", lty = 1, fill = "#c0c0c0") + shade_confidence_interval(endpoints = ci_boot, color = "red", lty = 3, fill = "#c0c0c0")

Step 7

  • Use simulation-based permutation to evaluate the difference in mean kernel95 home range size for males versus females. To do this, you will want to shuffle either the variable “sex” or “kernel95” home range size a total of 10,000 times, recalculating mean kernel95 size by sex for each permuted sample and then compare the difference in male and female kernel95 means from your original sample to the permutation distribution for the difference in means. We can do this with a for() loop…
obs_diff <- d |>
  summarize(diff = mean(kernel95[sex == "M"]) -
    mean(kernel95[sex == "F"])) |>
  pull(diff)

n_perm <- 10000 # number of permutations
# create a dummy vector to hold results for each permutation
perm_diff <- vector()
for (i in 1:n_perm) {
  # scramble the sex vector and then calculate difference between means
  # `sample()` with a vector as an argument yields a random permutation of the vector
  d |>
    mutate(sex = sample(sex)) |>
    summarize(diff = mean(kernel95[sex == "M"]) -
      mean(kernel95[sex == "F"])) |>
    pull(diff) -> perm_diff[[i]]
}
histogram(perm_diff)

p <- (sum(perm_diff >= abs(obs_diff)) +
  sum(perm_diff <= -abs(obs_diff))) / n_perm
# or
p <- sum(abs(perm_diff) >= abs(obs_diff)) / n_perm
p
## [1] 2e-04
ladd(
  panel.abline(v = obs_diff, col = "red", lty = 3, lwd = 2)
)

We can also use the do()* or replicate() functions to do our permutations…

# using do()*...
perm_diff <-
  {
    do(n_perm) *
      d |>
        mutate(sex = sample(sex)) |>
        summarize(diff = mean(kernel95[sex == "M"]) -
          mean(kernel95[sex == "F"]))
  } |>
  pull(diff) # extract column of value as vector
histogram(perm_diff)

# using replicate()... note the braces `{}` wrapped around the code to perform in each replicate within the `replicate()` function
perm_diff <- replicate(n_perm, {
  d |>
    mutate(sex = sample(sex)) |>
    summarize(diff = mean(kernel95[sex == "M"]) -
      mean(kernel95[sex == "F"])) |>
    pull(diff) # extract column of value as vector
})
histogram(perm_diff)

And, again, we could do this whole process using the {infer} package workflow…

perm_diff <- d |>
  specify(formula = kernel95 ~ sex) |>
  hypothesize(null = "independence") |>
  generate(reps = n_perm, type = "permute") |>
  calculate(stat = "diff in means", order = c("M", "F"))

obs_diff <- d |>
  specify(formula = kernel95 ~ sex) |>
  calculate(stat = "diff in means", order = c("M", "F"))
obs_diff
## Response: kernel95 (numeric)
## Explanatory: sex (factor)
## # A tibble: 1 × 1
##    stat
##   <dbl>
## 1  111.
perm_diff |>
  get_p_value(obs_stat = obs_diff, direction = "both")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1  0.0006
visualize(perm_diff) +
  shade_p_value(obs_stat = obs_diff, size = 0.5, fill = "#c0c0c0", direction = "both")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the infer package.
##   Please report the issue at <https://github.com/tidymodels/infer/issues>.

  • Under this approach, what is the “null” hypothesis? What is the test statistic? Is the difference in mean kernel95 home range size “significant”?

The null hypothesis here is that the difference in mean kernel95 home range size for males versus females is zero. The test statistic is the actual difference in mean kernel95 home range size for males versus females, and it is evaluated relative to a permutation distribution for this statistic. Based on our very low p value (less than 0.001), under a null hypothesis significance testing framework, we would conclude that the difference is “significant” (assuming that p is less than the alpha level we specify).

Step 8

  • Finally, use a theory-based parametric test (e.g., a t-test) to also calculate an appropriate test statistic and associated p value for the comparison of male and female mean kernel95 home range size. Is the difference in mean kernel95 home range size “significant”?
f <- d |> filter(sex == "F")
m <- d |> filter(sex == "M")

# first, what is ratio of variances in our two samples?
var_f <- var(f$kernel95)
var_m <- var(m$kernel95)
var_f / var_m
## [1] 1.272656
# ratio is less than 2.0, so we can use equal variance version of t test...

num_f <- nrow(f)
num_m <- nrow(m)

mean_f <- mean(f$kernel95)
mean_m <- mean(m$kernel95)

# hand-calculate the test statistic and p value...
s2 <- ((num_m - 1) * var_m + (num_f - 1) * var_f) / (num_m + num_f - 2)
t_stat <- (mean_m - mean_f) / sqrt(s2 * (1 / num_m + 1 / num_f))
t_stat
## [1] 3.94816
df <- num_m + num_f - 2
p <- 2 * (1 - pt(t_stat, df))
p
## [1] 0.0009426349
# or use the t.test() function
t_test <- t.test(x = m$kernel95, y = f$kernel95, var.equal = TRUE)
t_test
## 
##  Two Sample t-test
## 
## data:  m$kernel95 and f$kernel95
## t = 3.9482, df = 18, p-value = 0.0009426
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   51.91339 169.99890
## sample estimates:
## mean of x mean of y 
##  429.7407  318.7846
  • Under this approach, what is the test statistic? Is the difference in mean kernel95 home range size “significant”?

Here, the null hypothesis again is that the difference in mean kernel95 home range size for males versus females is zero. The test statistic is the actual difference in means scaled by something equivalent to a standard error that takes into account the variance and size of our two samples. This test statistic is then evaluated relative to a t distribution that depends on the number of degrees of freedom, which depends on the sample size of our two samples (if the variances were not roughly equal, we would calculate degrees of freedom taking those unequal variances into account). Again, based on the very low p value, we would conclude that the difference is “significant” under a null hypothesis significance testing framework (assuming that p is less than the alpha level we specify).