---
title: "08 Linear Regression. Crash Course in Statistics (Summer 2025)"
subtitle: "Neuroscience Center Zurich, University of Zurich"
author: "Zofia Baranczuk"
date: "2025-08-25"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 1. Example 1 
We consider the dataset teengamb (teenage gambling in Britain). Variables:

sex (0 = male, 1 = female)

status (socioeconomic status score)

income (pounds/week)

verbal (verbal score out of 12 correctly defined)

gamble (expenditure on gambling in pounds/year)

The data are available  in the package faraway.
### a) Load the data and perform exploratory data analysis
```{r}
library(faraway)
library(GGally)
data(teengamb)
str(teengamb)
summary(teengamb)
ggpairs(
  teengamb,
  columns = c("gamble", "income", "status", "verbal"),
  mapping = aes(color = factor(sex), alpha =0.4),   # or aes(colour = factor(sex))
  progress = FALSE
)
```

### b) Fit a multiple regression model with the expenditure on gambling as a the response and the sex, status, income, and verbal score as the predictors.
```{r}
lm_model <- lm(gamble ~ sex + status + income + verbal, data = teengamb)
summary(lm_model)
# Task (reminder from the theory part): Which predictors are significant, what is the meaning of p-values for predictors, what is the meaning of F-value. How to get the t-value when all other values in the summary are given? What do Multiple R-squared and	Adjusted R-squared mean?
```
### c) Assume, that the errors are independent. What are other assumptions of the linear model? Check these assumptions. Check for the influential points.
```{r}
#install.packages("ggfortify")
library(ggfortify)
autoplot(lm_model, which = 1:6, ncol = 3) # the ggplot version of diagnostic plots
plot(lm_model) # standard diagnostic plots
plot(lm_model, 4) # we can choose one of the standard plots, or also e.g. Cook's distance
# Discuss the diagnostic plots. 
library(car)
vif(lm_model)
```
### f) For all other predictors held constant, what would be the difference in predicted expenditure on gambling for a person with a status 5 compared to a person with status 4.
```{r}

coef(lm_model)["status"] 

```

### g) What is the 90% confidence interval of the income coefficient in this model?
```{r}

confint(lm_model, "income", level = 0.9)
#What is the meaning of the CI here? 

```

### h) Interpret the p-value for the intercept. Which hypothesis is being tested there?
The intercept has a p-value of 0.1968. This comes from a two-sided t-test with 42 degrees of freedom that tests the null hypothesis that the intercept equals 0 versus the alternative that it is not 0. The p-value is the probability, assuming the null is true, of getting an estimate at least as far from 0 as the one observed. Because 0.1968 is much larger than common alpha levels (0.05 or 0.01), we do not reject the null; there isn’t enough evidence that the intercept differs from 0.

### i) What is a predicted value of money spent on gambling for a female earning 8 pounds per week, having the status 50 and the verbal score 4?
```{r}
predict(lm_model, data.frame(sex = 1, status = 50, income = 8,
verbal = 4), interval = "confidence")

predict(lm_model, data.frame(sex = 1, status = 50, income = 8,
verbal = 4), interval = "prediction")

# Task: What is the difference between interval = "confidence" and interval = "prediction"?

```


### j) Use AIC criterion to check if you can use a model with less predictors. (hint: step())
```{r}

step_model <- step(lm_model)
summary(step_model)
#Discuss the outcomes, advantages and disadvantages of step. 


```
### k) Extra: Use bootstrapping to compute confidence interval of Income. Compare with the CI from the model. 
```{r}
set.seed(1) # what will happen if we don't set seed?
n <- nrow(teengamb) 
coefs <- numeric(n) 
B <- 1000 
 
for (i in 1:B) { 
  sample_indices <- sample(1:n, n, replace = TRUE) 
  sample_df <- teengamb[sample_indices, ] 
   
  model_boot <- lm(gamble ~ ., data = sample_df)
  coefs[i] <- coef(model_boot)["income"] 
} 
 
sol <- quantile(coefs, c(0.05, 0.95)) 
sol
# Task: look for a function in R that does bootstrapping and use the function to compute CI. 
# Compare the results. 


```

## 2. Is lm the right model? Discussion
Below we load three small datasets. After that, we discuss whether linear regression is appropriate and, if not, what to use instead.
```{r}

# Data we’ll use
library(faraway)   # orings, divusa (Divorce in the USA 1920–1996)
library(lme4)      # sleepstudy
library(ggplot2)

# a) O-rings: failures (counts out of 6) vs launch temperature
data(orings, package = "faraway")
str(orings); head(orings)
# quick look: proportion damaged vs temperature
ggplot(orings, aes(temp, damage/6)) +
  geom_point() + labs(y = "Proportion damaged", x = "Launch temperature (F)")

# b) Reaction times across days for multiple subjects
data(sleepstudy, package = "lme4")
str(sleepstudy); head(sleepstudy)
ggplot(sleepstudy, aes(Days, Reaction, group = Subject)) +
  geom_line(alpha = 0.4) + geom_point(alpha = 0.6)

# c) US divorce rates 
data(divusa, package = "faraway")
str(divusa)
pairs(divusa)

```

a) How does launch temperature relate to the probability of O-ring damage? Why or why not lm.
b) How does reaction time change with days of sleep restriction, allowing each subject their own baseline/slope?
c) Can we predict the divorce rate based on other sociodemographic factors?
