---
title: "Data Analysis 2"
author: "Ryan Womack"
date: "2025-09-15"
toc:toc: true true
number-sections: true
highlight-style: pygments
output: html_document
format:
  html: 
    toc: true
    code-fold: true
    html-math-method: katex
  pdf:
    geometry: 
      - top=30mm
      - left=30mm
    fontfamily: libertinus
    colorlinks: true
    papersize: A4
  docx: default
theme: litera
---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(root.dir = "/home/ryan/R/data_topics/data_analysis_2/")

```

Copyright Ryan Womack, 2025. This work is licensed under [CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/)

**Data Analysis 2**

_statistical tests, distributions, regression, bootstrap, and more + comparison with Python_

# Overview

This workshop reviews the implementation of basic statistical tests and methods in R, with discussion contextualizing the appropriate use of these tests. For comparison purposes, an alternative Python approach to performing similar statistical analysis is illustrated in some cases.  This workshop takes inspiration from the **Introduction to Modern Statistics** available via [github](https://github.com/OpenIntroStat/ims) and [online text](https://openintro-ims2.netlify.app/), part of the [OpenIntro](https://openintro.org) project.

Note that "packages" in R are the equivalent of "modules" in Python.

# Setup and preparing data

## R packages required

We will use the [_pak_](https://pak.r-lib.org/) package for installation as a more complete approach to package management. Replace the pkg commands with _install.packages()_ versions if you prefer.

This session relies on the [_tidyverse_](https://tidyverse.org) suite of packages for data manipulation as a preference, although the same tasks could be accomplished in base R. The statistical functions used are those from base R. Install _pak_ and _tidyverse_ if you don't already have them on your system. We will also need _reticulate_ to run Python code chunks.

```{r install packages, eval=FALSE}

install.packages("pak", dependencies=TRUE)
library(pak)
pkg_install("tidyverse")
pkg_install("reticulate")
pkg_install("infer")
pkg_install("TOSTER")
devtools::session_info()

```

Now let's load the tidyverse, infer, TOSTER, and reticulate (for Python support only)
```{r tidyverse}

library(tidyverse)
library(infer)
library(TOSTER)
Sys.setenv(RETICULATE_PYTHON = "/usr/bin/python3")
library(reticulate)

```

## Data import and preparation

Now let's grab some data. We will use a realistic data example, the [World Bank's Gender Statistics database](https://genderdata.worldbank.org/), whose [raw data](https://databank.worldbank.org/data/download/Gender_Stats_CSV.zip) is directly downloadable. Other [World Bank Open Data](https://data.worldbank.org/) is available as well. See [genderdata.worldbank.org](https://genderdata.worldbank.org) for more background on the Gender Data portal. Note that we have to inspect the data and understand the variables first before manipulating in R -- this is not an automatic process.

```{r download and import data}

getOption("timeout")
options(timeout=6000)
download.file("https://databank.worldbank.org/data/download/Gender_Stats_CSV.zip", "gender.zip")
unzip("gender.zip")
gender_data <- read_csv("Gender_StatsCSV.csv")

```

Now we'll perform a few steps to clean the data, focusing on generating a useable file for a few countries (Central Asia and Mongolia plus selected high population or high income countries), from the latest available data year with complete data, typically the year before the last in the data set. For this session, we'll just run these steps without explaining them. Data cleaning and wrangling is covered in more detail in the Data Analysis 1 workshop (forthcoming). The final filtered output is the _gender_data_final_ file, which we _attach_ so that a copy of the data is made the default dataset for this session.

```{r data wrangling, results='hide'}

# clean the data to remove superfluous columns
names(gender_data)
gender_data <- gender_data[,c(-2,-4)]
names(gender_data)

# select countries of interest
country_list <- c("China", "Germany", "India", "Japan", "Kazakhstan", "Kyrgyz Republic", "Mongolia", "Russian Federation", "Tajikistan", "Turkmenistan", "United States", "Uzbekistan")
gender_data2 <-
  gender_data %>%
  filter(`Country Name` %in% country_list)

# clean the data to focus on a recent more complete time period
gender_data3 <-
   gender_data2 %>%
   pivot_longer(3:67, names_to = "Year", values_to = "Value")

#filter by year
gender_data2023 <-
  gender_data3 %>%
  filter(Year=="2023")

gender_data2023 <- gender_data2023[,-3]

gender_data2023wide <-
  gender_data2023 %>%
  pivot_wider(names_from = "Indicator Name", values_from = "Value")

# now use a little sapply trick to select variables that don't have much missing data - here the proportion is 0.75 (the 0.25 in the function is 1-proportion desired)

gender_data_filtered <- gender_data2023wide[,!sapply(gender_data2023wide, function(x) mean(is.na(x)))>0.25]

# and lastly simplify the dataset by removing some of the topics we won't use

phrases <- c("Worried", "Made", "Received", "Saved", "Used", "Coming", "Borrowed")

gender_data_final <- 
  gender_data_filtered %>%
  select(!starts_with(phrases))

# we'll also generate a couple of variables for future use

gender_data_final$female_high_labor <- gender_data_final$`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`>70
gender_data_final$male_high_labor <- gender_data_final$`Labor force participation rate, male (% of male population ages 15-64) (modeled ILO estimate)`>70 

attach(gender_data_final)

```

Let's test for some basic relationships between labor participation, gender, fertility, and income, using the [World Bank's own definition of gender](https://genderdata.worldbank.org/en/about): 

>"Gender refers to the social, behavioral, and cultural attributes, expectations, and norms associated with being male or female."

Note that this is just for the purposes of demonstration, and not a serious investigation into these important research issues.

# Statistical Tests

## t-test

In statistics there are many **hypothesis tests**. Per the Wikipedia entry on [statistical hypothesis tests](https://en.wikipedia.org/wiki/Statistical_hypothesis_test), 

>"a statistical hypothesis test is a method of statistical inference used to decide whether the data sufficiently supports a particular hypothesis. A statistical hypothesis test typically involves a calculation of a test statistic. Then a decision is made, either by comparing the test statistic to a critical value or equivalently by evaluating a p-value computed from the test statistic. Roughly 100 specialized statistical tests have been defined."

The hypothesis is usually framed as a "Null hypothesis", describing the situation where there is no statistically significant difference, and the alternative hypothesis. Then a test is chosen appropriate to the situation, which often involves invoking a statistical distribution.

The *t-test* is one of the most common tests, used to determine if the difference between two groups, according to some numerical measure, is statistically significant. The Null hypothesis, usually denoted $H_{0}$, is that there is no difference (a difference of zero). The Alternative, usually denoted $H_{1}$ or $H_{A}$, is that the difference is not zero.  

The t-test was originally, and sometimes still is called ["Student's t-test"](https://en.wikipedia.org/wiki/Student's_t-test). The story of the student who invented this test, and his affiliation with a certain well known Irish stout, is interesting from both a statistical and human perspective. The $t$-test has a few variants: a one sample test, a two sample test with unpaired results, and two sample test with paired results.  The one sample test simply checks whether one measure is significantly different from the null. The two sample test with unpaired results compares whether two separate sets of observations are different, such as sampling the populations of two different cities. The two sample test with paired results implies that we have the same subjects in the dataset who are measured twice, perhaps at different intervals or via different measures. The $t$-distribution is the underlying statistical distribution determining the test statistic and critical value that we check for. These are conveniently summarized by the *p-value*, which expresses the level of significance. We typically look for a $p$<.05 to determine statistical significance, but this is a convention: other $p$-value cutoffs can be used.

## One Sample t-test

We start with a one sample $t$-test to check the labor force participation rates of females and males across the countries in our dataset.

```{r t-test}

t.test(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`)

```
This first example returns a statistically significant, but not very interesting result. The default is to test whether the variable is different than zero. It is not surprising that female labor force participation is higher than zero.

```{r t-test mu}

t.test(`Labor force participation rate, male (% of male population ages 15-64) (modeled ILO estimate)`, mu=70)

```
We refine the test by passing the option $mu$=70, to test whether the male labor participation rate is significantly different than 70. With a $p$>.05 and a 95% confidence interval of (67,83), it is *not* significantly different than 70 in a statistical sense. 

Note the use of options is a common R syntax technique. Base execution of the command without options gives sensible results, but one can pass many options to tweak the function's behavior. Most of the things demanded by statisticians are available through these option tweaks. Typing _?t.test_ (the question mark followed by the function name) will pull up the help that displays the possibilities.

```{r help}

?t.test

```

## Paired t-test

Perhaps a more interesting question is to see if there is a significantly significant difference in male and female labor participation rates. We can compare them, country-by-country, using the paired=TRUE option, in the paired two-sample test below. For one explanation of this see [STHDA](http://www.sthda.com/english/wiki/paired-samples-t-test-in-r).

```{r t-test paired}

t.test(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`,`Labor force participation rate, male (% of male population ages 15-64) (modeled ILO estimate)`, paired = TRUE, alternative = "two.sided")

```
There is a statistically significant difference, with the 95% confidence interval of female participation being between 7 and 25 percentage points lower than male rates, and a mean difference of -16 points.

## Python t-test

We will be giving a flavor of the Python approach to these problems, without going in depth. Python is covered in other workshops in greater detail.

Running Python in RStudio requires the _reticulate_ package. Your Python installation should also have _numpy_ and _scipy_ installed. Python will not have direct access to the R data structures, so we use a simplified example here from [Builtin.com](https://builtin.com/data-science/t-test-python)

How to pass data from R to Python and vice-versa is described in this [blog post by Dima Diachkov](https://medium.com/data-and-beyond/how-to-seamlessly-integrate-python-into-r-rmarkdown-codes-2fe09cfdd0ee)

Note that our code chunk is labeled Python. For more _reticulate_ package, see [this post](https://www.r-bloggers.com/2022/04/getting-started-with-python-using-r-and-reticulate/)

Also note that it is also possible to invoke R from a Python environment using [rpy2](https://pypi.org/project/rpy2/)

```{python t-test python}

import numpy as np
import pandas as pd
from scipy import stats

# import from R
gender_python = r.gender_data_final

# print(gender_python)

labor = gender_python.loc[:,"Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)"]

# Hypothesized population mean
mu = 70

# Perform one-sample t-test
# t_stat, p_value = stats.ttest_1samp(labor, mu)
t_stat, p_value = stats.ttest_1samp(labor, mu)
print("T statistic:", t_stat)
print("P-value:", p_value)

# Setting significance level
alpha = 0.05

# Interpret the results
if p_value < alpha:
    print("Reject the null hypothesis; there is a significant difference between the sample mean and the hypothesized population mean.")
else:
    print("Fail to reject the null hypothesis; there is no significant difference between the sample mean and the hypothesized population mean.")

```
## Chi-squared test (in R and Python)

Many operations in R are a question of finding the appropriate function, either using the built-in R help via the ? (for help on a specific command) or ??  (to search across functions by keywords), or by using the built-in help system. Using a search engine like DuckDuckGo for a topic plus "in R" or a help site like stackoverflow.com will also turn up useful leads. Or one can consult books such as the many [R-related book series](https://link-springer-com.proxy.libraries.rutgers.edu/search?new-search=true&query=Use+R&dateFrom=&dateTo=&sortBy=relevance) in [SpringerLink](https://login.proxy.libraries.rutgers.edu/login?url=https://link.springer.com) for more comprehensive how-to's.  I would not recommend using AI tools for anything important, since one must understand the results well enough to check them for errors. 

So, for example, if we wanted to run a [chi-squared test](https://en.wikipedia.org/wiki/Chi-squared_test) instead of a t-test, some noodling about will turn up the _chisq.test_ function. We won't discuss the distribution or mathematical background of the chi-squared test here.

Take a look at this table of the number of countries with "high" labor participation rates of over 70%, for males and females.

```{r labor}

table(female_high_labor,male_high_labor)

```

We can test whether this is a statistically significant pattern via the chi-squared test, implemented with:

```{r chisq}

chisq.test(table(female_high_labor,male_high_labor))

```

Which is *not* statistically significant, probably due to the small size of the matrix.

## How to choose the appropriate statistical test?

This is a question that depends on many factors, that we won't go into detail about in this introductory workshop, but you can consult the following resources:

- [UCLA's quick guide](https://stats.oarc.ucla.edu/other/mult-pkg/whatstat/)
- [a nice brief summary article](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8327789/)
- You can also consult [Sage Research Methods' Which Stats Test?](https://methods-sagepub-com.proxy.libraries.rutgers.edu/which-stats-test) for more information on choosing the appropriate statistical test. The other guides in [Sage Research Methods](https://methods-sagepub-com.proxy.libraries.rutgers.edu/) are good too for more in-depth exploration of specific tests and topics! - _Rutgers-restricted_

## Distributions

One great aspect of R is the availability of tools to work with almost any statistical distribution. The density, distribution function, quantile function and random generation of a number are easily available with R functions, via the prefixes _d, p, q, and r_. For example for the normal distribution, we have:

```{r normal distribution}

# the density of a standard normal distribution at the value 2 (2 above mean of zero)

dnorm(2)

# the cumulative percentage of a standard normal distribution below the value 2

pnorm(2)

# mean and standard deviation of the distribution can also be specified
# in this example the cumulative percentage of a normal distribution with mean 100 and s.d. 20, below the value 90

pnorm(90, mean=100, sd=20)

# provides the numeric value of a particular quantile of the distribution (again standard normal in this example)

qnorm(.9)

# a single random draw from the standard normal distribution

rnorm(1)

# five random draws from a normal distribution with mean 100 and s.d. 20

rnorm(5, mean=100, sd=20)

# getting help

?rnorm

```

As usual the question mark will provide details about implementation, variables, and options.

Other distributions and their functions can be found on the [TaskViews on Distributions](https://cran.r-project.org/web/views/Distributions.html) page.

# Correlation

We can also perform simple correlations to evaluate the strength of a relationship between two variables, but there are some cautions.

```{r correlation incorrect, eval=FALSE}

cor(`GDP per capita (constant 2015 US$)`,`Fertility rate, total (births per woman)`, na.rm=TRUE)

```

Why won't this work? The defaults and options of R commands are not very standardized.  Instead try this option:

```{r correlation correct}

cor(`GDP per capita (constant 2015 US$)`,`Fertility rate, total (births per woman)`, use="complete.obs")

```

The _cor.test_ function gives more complete output 

```{r cor.test}

cor.test(`GDP per capita (constant 2015 US$)`,`Fertility rate, total (births per woman)`)

```

# Regression

*Regression* is a common and fundamental technique to model a relationship between a *response* variable (also called dependent or outcome variable) and one or more *explanatory* variables (also called independent or predictor variables). Proper implementation of regression requires careful attention to the data and examination of the model fit, variable selection, and more. The quick and dirty approach below is simply designed to show R syntax in action.

## Linear Regression

We start at the beginning, with linear regresion and the _lm_ command.  Let's see if there is a relationship between female labor participation and GDP:

```{r linear regression}

lm(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`GDP per capita (constant 2015 US$)`)

```

Which predicts about a 0.5 increase in the percentage labor participation of females for every $1000 rise in GDP. However, the presentation of the result is a bit odd. R only presents the coefficients by default. The _summary_ command must be used to tease out the additional information that we expect in a regression table to evaluate fit and significance.

```{r linear regression with summary}

summary(lm(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`GDP per capita (constant 2015 US$)`))

```

Which is significant at the 0.05 level. Note that we must decide on the appropriate $p$-value to test for in advance and resist the temptation to say that this is "almost significant".  Note that this is heavily influenced by the fact that we are using a small sample of 12 countries, which could be an argument to use a 0.10 signficance level.

Let's check some additional relationships:

```{r linear regression additional relationships}

summary(lm(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`Fertility rate, total (births per woman)`))

summary(lm(`Fertility rate, total (births per woman)`~`GDP per capita (constant 2015 US$)`))

```

Here our p-values are significant at the 10% level but not the 5% level. Also note that if want to run the regression without an intercept, we can just add a -1 to the equation.

```{r linear regression with no intercept}

summary(lm(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`GDP per capita (constant 2015 US$)`-1))

```

## Multiple Regression

And if we need to explore *multiple regression*, with multiple explanatory variables, this is as easy as using a + sign in the equation:

```{r multiple linear regression}

summary(lm(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`Fertility rate, total (births per woman)`+`GDP per capita (constant 2015 US$)`))

```

## Stored Regression Objects

A very important feature of R is the flexibility that derives from the ability to store (and use) our regression output as an R object.  The output becomes an R object, and we can access its variables/components using the familiar $ notation.  This can allow us to easily store and compare multiple models and use regression data such as residuals or predicted values as inputs into other equations and analyses.

```{r stored regression output}

regoutput<-lm(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`Fertility rate, total (births per woman)`+`GDP per capita (constant 2015 US$)`)
names(regoutput)
regoutput$residuals

```
Once we have store the regression output, numerous quick functions are available that build off of the regression results, such as _predict_ for predicted values and _anova_ for the Analysis of Variance.

```{r stored regression output - quick functions}

# predicted values
predict(regoutput)

# analysis of variance (anova)
anova(regoutput)

```

## Regression Diagnostics

_Regression diagnostics_ are used to evaluate the fit and appropriateness of the regression model as applied. Once again, we are providing some illustrative examples without going into depth here.

One thing that we may wish to test is the normality of the regression residuals (i.e., the difference between predicted and actual values of the observations is normally distributed). See this [discussion of Normality at U Wisconsin](https://sscc.wisc.edu/sscc/pubs/RegDiag-R/normality.html) and this [article available in PubMedCentral](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3693611/) for a brief explanation.

We can apply a typically used test for normality, the Shapiro-Wilk Normality Test, as follows. We get standardized residuals from the model with _rstandard_ (or more simply the residuals accessed with _regoutput$residuals_), and then use the _shapiro.test_ function to perform the test, or the (less powerful) Kolmogorov-Smirnov Test with _ks.test_, which also allows us to test against different distributions (here we just specific _'pnorm'_ for the normal distribution).

A $p$<.05 in both the Shapiro-Wilk and Kolmogorov-Smirnov test implies _rejection_ of normality. But note in the example below that we must be careful to standardize the residuals so that the test is applied correctly.

```{r tests for normality}

shapiro.test(rstandard(regoutput))

ks.test(regoutput$residuals, 'pnorm')

ks.test(rstandard(regoutput), 'pnorm', mean=0, sd=1)

```
Note that standardizing the residuals makes a big difference for the Kolmogorov-Smirnov test!

We can also easily pull up regression diagnostic plots:

```{r diagnostic plots}

plot(regoutput, pch=3)

```

## glm and Logistic Regression

R is developed by statisticians for statisticians, so you will find every variant of regression that you might eventually need. The _glm_ command implements generalized linear models (or alternatively, the [_glm2_](https://cran.r-project.org/package=glm2) package).

One commonly used alternative, necessary where the response variable is binary (categorical) is *logistic regression*. The _glm_ command is used, with "family=binomial" as an argument. The binomial distribution represents the fact that the response variable is binary (either a two-level factor or directly coded as 0/1 - referred to as "one-hot" coding in certain circles).

We can try this out with our categorized variable that labels female labor participation as high if over 70% (otherwise low):

```{r logistic regression}

logistic_output <- glm(female_high_labor ~ `Fertility rate, total (births per woman)`+`GDP per capita (constant 2015 US$)`, family=binomial)

summary(logistic_output)

```

We won't delve further into the interpretation of logistic regression or other regression alternatives here, but just know that nearly anything is possible with R!

## Python regression

Linear Regression is part of the [scikit-learn](https://scikit-learn.org/) module which forms the core of machine learning in Python.

Some basic information on linear regression in Python is at [RealPython](https://realpython.com/linear-regression-in-python/#simple-linear-regression)

We include this code as a quick illustration of the approach. As is typical, the output in Python must be individually extracted with commands.  The [stargazer](https://github.com/StatsReporting/stargazer) module is at least one approach to automating the process of generating output tables.  The [statsmodels](https://www.statsmodels.org/stable/index.html) module is less commonly used, but provides an alternative in Python that generates output tables more easily.

```{python regression in python}

import numpy as np
import pandas as pd
from scipy import stats
from sklearn.linear_model import LinearRegression

# import from R
gender_python = r.gender_data_final

# specify variables
labor = gender_python.loc[:,"Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)"]
gdp = gender_python.loc[:,"GDP per capita (constant 2015 US$)"]

# this step is important to put the data in Python-friendly form
labor2 = labor.array.reshape(-1, 1)
gdp2 = gdp.array.reshape(-1, 1)
model = LinearRegression().fit(labor2,gdp2)

print(model.intercept_, model.coef_)
r_sq = model.score(labor2, gdp2)
print(f"R-squared: {r_sq}")

```

# Bootstrap

This section is modeled off the approach used by [Introduction to Modern Statistics](https://openintro-ims.netlify.app/), Second Edition by Mine Çetinkaya-Rundel and Johanna Hardin, which goes into greater depth and provides many useful code exercises. We'll refer to this text as *IMS* below. As elsewhere, we are providing a quick and dirty demo, but please use these resources to go further!

In particular take a look at the section on [bootstrapping](https://openintro-ims.netlify.app/foundations-bootstrapping) and the [supplementary tutorials](https://openintrostat.github.io/ims-tutorials/).

Why the "bootstrap"? Many traditional statistical techniques rely on assumptions about the behavior or distribution of the data, assumptions that we cannot always be sure of satsifying. According to IMS, 
"some statistics do not have simple theory for how they vary, and bootstrapping provides a computational approach for providing interval estimates for almost any population parameter."

The bootstrap works by sampling and resampling the original collection of data.  The variability of the samples can be used to estimate the variability of the underlying population, without requiring any further assumptions about its characteristics.

Essentially, we resample our data to construct an empirical distribution, which we can then use to estimate standard errors and corresponding confidence intervals, as well as other parameters. We call it the "bootstrap" because we are "pulling ourselves up by our bootstraps" and creating the standard error ourselves via this replicate sampling process.

R makes this easy via the [_infer_](https://infer.netlify.app/) package from the tidyverse.

The bootstrap in R is constructed by a three-part sequence of commands: _specify_, _generate_, and _calculate_:
 - we _specify_ the variable we want to sample
 - then _generate_ the number of sample replicates we want to generate
 - then _calculate_ the statistic of interest on the replicates

Before any sampling we can simply generate the point estimate of our proportion, in this case one of our indicator variables for the gender data, whether a woman can work in dangerous jobs in the same way as a man.

Note that |> is the "pipe" in base R (relatively newly introduced, so you won't see it as often yet), versus the %>% "pipe" in tidyverse via _magrittr_ (traditional).

```{r point_estimate}

point_estimate <- gender_data_final %>%
  specify(response = `A woman can work in a job deemed dangerous in the same way as a man (1=yes; 0=no)`) %>%
  calculate(stat = "mean")

point_estimate

```

Then run replicate samples that we will use to estimate the confidence interval.

```{r bootstrap confidence interval}

boot_dist_dangerous <- gender_data_final %>%
  specify(response = `A woman can work in a job deemed dangerous in the same way as a man (1=yes; 0=no)`) %>%
  generate(reps = 50000, type = "bootstrap") %>%
  calculate(stat = "mean")
  
```

We can easily visualize the distribution of the replicate samples.

```{r bootstrap replicates}

boot_dist_dangerous %>%
  visualize()

boot_dist_dangerous %>%
  visualize() +
  shade_p_value(obs_stat = point_estimate, direction = "two-sided")
  
boot_dist_dangerous %>%
  get_confidence_interval(
    point_estimate = point_estimate,
    level = 0.95,
    type = "se"
  )

```

Note that we are resampling from a small population of 12, so even with large replications, our confidence interval is wide. Using this technique on the type of country data we have is atypical, and is just for demonstration purposes. A more typical usage would be on survey data with individual responses selected from a larger population. Again to quote *IMS*, "bootstrapping is best suited for modeling studies where the data have been generated through random sampling from a population." 

## t-test via simulation

*IMS* has some examples of using _infer_ to simulate t-test type results, however this is not as well-suited to the case of a paired t-test, for which the TOSTER package has a built in function that quickly implements a bootstrapped version of the t-test. This is just one example of how we can switch among the many ways of accomplishing a task in R.

First, instead of preserving the pairs, we get at the same idea by computing the difference in labor participation for males and females for each country.

```{r labor_diff}

labor_diff <- `Labor force participation rate, male (% of male population ages 15-64) (modeled ILO estimate)` - `Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`

labor_diff

```

Then apply the bootstrap t-test from the _TOSTER_ package, specifying 1000 replications.

```{r boot_t_test}

boot_t_test(labor_diff, R=1000)

```

## bootstrap to estimate variability of the slope of the regression 

We can also bootstrap our way to some of the results of a regression model. Here we also use the _hypothesize_ option of the _infer_ package. Note that for this small sample size, the bootstrap estimate fluctuates quite a bit if you repeat the experiment.


```{r boot slope}

var_slope <- gender_data_final |>
  specify(`Labor force participation rate, female (% of female population ages 15-64) (modeled ILO estimate)`~`Fertility rate, total (births per woman)`) |>
  hypothesize(null = "independence") |>
  generate(reps = 500, type = "permute") |>
  calculate(stat = "slope")

var_slope |> 
  # Ungroup the dataset
  ungroup() |> 
  # Calculate summary statistics
  summarize(
    # Mean of stat
    mean_stat = mean(stat), 
    # Std error of stat
    std_err_stat = sd(stat)
  )  

```

## Python bookstrap

Finally we'll briefy look at the [Python approach to bootstrapping](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bootstrap.html).

```{python bootstrap}

from scipy.stats import bootstrap
import numpy as np

labor = r.labor_diff

#convert array to sequence
data = (labor,)

#calculate 95% bootstrapped confidence interval for median
bootstrap_ci = bootstrap(data, np.mean, confidence_level=0.95,
                         random_state=1, method='percentile')

#view 95% boostrapped confidence interval
print(bootstrap_ci.confidence_interval)

```

There's [more that can be done](https://machinelearningmastery.com/calculate-bootstrap-confidence-intervals-machine-learning-results-python/) of course, but we'll stop here.

*Enjoy R!* _(and Python!)_