---
title: "Regression Analysis"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Regression Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

```{r}
library(mariposa)
library(dplyr)
data(survey_data)
```

## Overview

Regression analysis predicts an outcome from one or more predictors. mariposa provides two regression functions with SPSS-compatible output:

| Function | Use when |
|----------|----------|
| `linear_regression()` | Outcome is continuous (e.g., income, satisfaction score) |
| `logistic_regression()` | Outcome is binary (e.g., yes/no, high/low) |

Both functions support two interface styles:

- **Formula**: `linear_regression(data, y ~ x1 + x2)` --- standard R syntax
- **SPSS-style**: `linear_regression(data, dependent = y, predictors = c(x1, x2))`

## Linear Regression

### Simple Regression

```{r}
linear_regression(survey_data, life_satisfaction ~ age)
```

### Detailed Output

```{r}
result <- linear_regression(survey_data, life_satisfaction ~ age)
summary(result)
```

The detailed output includes four sections matching SPSS REGRESSION:

- **Model Summary**: R, R-squared, Adjusted R-squared
- **ANOVA Table**: Overall model significance
- **Coefficients**: B (unstandardized), Beta (standardized), t, p, confidence intervals
- **Descriptives**: Mean and SD for all variables

### Understanding Coefficients

- **B (unstandardized)**: For each 1-unit increase in the predictor, the outcome changes by B units
- **Beta (standardized)**: Allows comparison across predictors on different scales. A Beta of 0.30 means a 1-SD increase in the predictor is associated with a 0.30-SD change in the outcome
- **p-value**: Below 0.05 indicates statistical significance

### Multiple Regression

```{r}
linear_regression(survey_data,
                  life_satisfaction ~ age + income + trust_government)
```

Compare Beta values to identify the strongest predictor.

### SPSS-Style Interface

```{r}
linear_regression(survey_data,
                  dependent = life_satisfaction,
                  predictors = c(age, income, trust_government))
```

### With Survey Weights

```{r}
linear_regression(survey_data,
                  life_satisfaction ~ age + income,
                  weights = sampling_weight)
```

Weights are treated as frequency weights, matching SPSS `WEIGHT BY` behavior.

### Grouped Analysis

Run separate regressions for each subgroup:

```{r}
survey_data %>%
  group_by(region) %>%
  linear_regression(life_satisfaction ~ age + income)
```

### Interpreting R-squared

R-squared tells you how much variance the predictors explain:

- 0.01 -- 0.05: Small effect
- 0.06 -- 0.13: Medium effect
- 0.14+: Large effect

These benchmarks follow Cohen (1988). Always check the ANOVA table to confirm overall model significance.

### Using Transformed Predictors

Combine with data transformation functions for better models:

```{r}
# Standardize predictors for comparable coefficients
survey_data_z <- survey_data %>%
  std(age, income, suffix = "_z")

linear_regression(survey_data_z,
                  life_satisfaction ~ age_z + income_z + trust_government,
                  weights = sampling_weight)
```

## Logistic Regression

### When to Use

Use `logistic_regression()` when your outcome is binary. First, create a binary variable:

```{r}
survey_data <- survey_data %>%
  mutate(high_satisfaction = ifelse(life_satisfaction >= 4, 1, 0))
```

### Basic Logistic Regression

```{r}
logistic_regression(survey_data, high_satisfaction ~ age + income)
```

### Detailed Output

```{r}
log_result <- logistic_regression(survey_data, high_satisfaction ~ age + income)
summary(log_result)
```

The detailed output includes five sections matching SPSS LOGISTIC REGRESSION:

- **Omnibus Test**: Overall model significance
- **Model Summary**: -2 Log Likelihood and pseudo R-squared values
- **Hosmer-Lemeshow Test**: Model fit assessment
- **Classification Table**: Prediction accuracy
- **Coefficients**: B, Wald, Exp(B) (odds ratios), confidence intervals

### Understanding Odds Ratios

Exp(B) is the odds ratio --- the key statistic in logistic regression:

- **Exp(B) > 1**: Each unit increase raises the odds (e.g., 1.50 = 50% higher odds)
- **Exp(B) < 1**: Each unit increase lowers the odds (e.g., 0.80 = 20% lower odds)
- **Exp(B) = 1**: No effect

### Multiple Predictors

```{r}
logistic_regression(survey_data,
                    high_satisfaction ~ age + income + trust_government + education)
```

### SPSS-Style Interface

```{r}
logistic_regression(survey_data,
                    dependent = high_satisfaction,
                    predictors = c(age, income, trust_government))
```

### With Survey Weights

```{r}
logistic_regression(survey_data,
                    high_satisfaction ~ age + income,
                    weights = sampling_weight)
```

### Grouped Analysis

```{r}
survey_data %>%
  group_by(region) %>%
  logistic_regression(high_satisfaction ~ age + income)
```

### Interpreting Model Fit

**Pseudo R-squared** values are not directly comparable to linear regression R-squared:

- **Nagelkerke R-squared**: Adjusted to reach 1.0, most commonly reported
- **Cox & Snell R-squared**: Cannot reach 1.0, always lower
- **McFadden R-squared**: Values above 0.20 indicate good fit

**Hosmer-Lemeshow Test**: Non-significant ($p > .05$) means the model fits well.

**Classification Table**: Compare correct predictions to the base rate --- your model should outperform guessing the most common category.

### Average Marginal Effects

Odds ratios are hard to communicate. `marginal_effects()` (Stata's
`margins, dydx(*)`) translates the model onto the probability scale:
the average change in predicted probability per unit of each predictor,
with delta-method standard errors. Factor predictors show the average
discrete change against their reference level:

```{r}
model <- survey_data %>%
  logistic_regression(high_satisfaction ~ age + income + education)

marginal_effects(model)
summary(marginal_effects(model))
```

An AME of 0.03 reads: "one more unit increases the probability of the
outcome by 3 percentage points, averaged over the sample" --- usually
the sentence your readers actually need.

## Complete Example

```{r}
# 1. Explore relationships first
survey_data %>%
  pearson_cor(life_satisfaction, age, income, trust_government)

# 2. Run linear regression
lm_result <- linear_regression(survey_data,
                               life_satisfaction ~ age + income + trust_government,
                               weights = sampling_weight)
lm_result
summary(lm_result)

# 3. Create binary outcome
survey_data <- survey_data %>%
  mutate(high_satisfaction = ifelse(life_satisfaction >= 4, 1, 0))

# 4. Run logistic regression
log_result <- logistic_regression(survey_data,
                                  high_satisfaction ~ age + income + trust_government,
                                  weights = sampling_weight)
log_result
summary(log_result)
```

## Practical Tips

1. **Check correlations first.** Use `pearson_cor()` to explore bivariate relationships before building a model.

2. **Center or standardize predictors.** Centering makes the intercept interpretable; standardizing (with `std()`) makes Beta coefficients comparable. Use `std(method = "2sd")` to make continuous predictors comparable to binary ones.

3. **Compare Beta values.** In multiple regression, standardized Beta coefficients reveal which predictor has the strongest effect, regardless of measurement scale.

4. **Watch for multicollinearity.** Highly correlated predictors produce unstable coefficients. Check bivariate correlations before interpreting results.

5. **Never use `linear_regression()` with a binary outcome.** Predicted values can fall outside 0--1, and the statistical tests are invalid. Use `logistic_regression()` instead.

6. **Report completely.** Include R-squared, F-test or Omnibus test, individual coefficients, and sample size.

## Summary

1. `linear_regression()` predicts continuous outcomes (B, Beta, ANOVA, R-squared)
2. `logistic_regression()` predicts binary outcomes (odds ratios, classification, pseudo R-squared)
3. Both support formula and SPSS-style interfaces, survey weights, and grouped analysis
4. Use `std()` and `center()` to prepare predictors for better interpretability
5. Always explore bivariate relationships before building regression models

## Next Steps

- Prepare variables with `std()` and `center()` --- see `vignette("data-transformation")`
- Construct reliable scale scores --- see `vignette("scale-analysis")`
- Compare groups directly --- see `vignette("hypothesis-testing")`
- Handle survey weights --- see `vignette("survey-weights")`