--- title: "Prediction" author: Maureen Ritchey format: html: df-print: kable theme: cosmo toc: true toc-depth: 3 embed-resources: true execute: echo: true editor: source --- # Introduction Multiple regression is a statistical technique used to understand how multiple predictor variables influence an outcome variable. It is widely used in **psychology and neuroscience** to examine relationships between brain activity, cognitive performance, and behavioral variables. # Load packages ```{r, message=FALSE} # Install packages if not already installed # install.packages(c("tidyverse", "caret", "broom")) # Load packages library(tidyverse) library(caret) library(broom) ``` # Understand the data We’ll use the built-in `mtcars` dataset, predicting **miles per gallon (`mpg`) based on horsepower (`hp`) and weight (`wt`)**. ```{r} # Load dataset data(mtcars) # View first few rows head(mtcars) ``` # Fit a multiple regression model The regression equation is: $$ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon $$ Where: - $Y$ = Outcome variable (`mpg`) - $X_1, X_2$ = Predictor variables (`hp` and `wt`) - $\beta_0$ = Intercept (baseline when predictors are 0) - $\beta_1, \beta_2$ = Coefficients (effects of predictors) - $\varepsilon$ = Error term ```{r} # Fit the linear model model <- lm(mpg ~ hp + wt, data = mtcars) # View model summary summary(model) ``` ## Interpreting regression output When we run `summary(model)`, we get key information about our regression model: ### Residuals Residuals represent the differences between the actual and predicted values. Ideally, they should be symmetrically distributed. ### Coefficients The table shows: - **Estimate**: The effect of each predictor. - **Std. Error**: The uncertainty in the estimates. - **t-value & p-value**: Statistical significance of predictors. ### Model Fit - **R-squared**: Measures how well the model explains variance. - **F-statistic**: Tests whether the model is statistically significant. For our model: - `hp` and `wt` are **significant predictors** of `mpg`. - The model explains **82.68% of the variance**. - The **F-test (p \< 0.0001)** confirms overall significance. # Generate predicted values ```{r} # Generate predictions on the original dataset mtcars$predicted_mpg <- predict(model, newdata = mtcars) # View actual vs. predicted values head(mtcars[, c("mpg", "predicted_mpg")]) ``` # Perform cross validation ## What is Cross-Validation and Why Use It? Cross-validation helps **assess how well a model generalizes to new data**. Instead of evaluating a model on the same data it was trained on (which can lead to overfitting), we split the data into multiple folds, train the model on some folds, and test it on others. This is especially important in neuroscience and psychology, where **sample sizes are often small**, and we want to ensure our model isn't just capturing random noise. ## 5-Fold Cross-Validation ```{r} # Set up cross-validation using the caret package set.seed(123) # For reproducibility train_control <- trainControl(method = "cv", number = 5) # Train model with cross-validation - enter ?train for more information about the settings cv_model <- train(mpg ~ hp + wt, data = mtcars, method = "lm", trControl = train_control) # View cross-validation results print(cv_model) ``` ## Understanding the `cv_model` output Let’s break down each part: ### Model Details - **Linear Regression**: Confirms the type of model being used. - **32 samples**: The dataset has 32 observations. - **2 predictor variables**: The model includes `hp` (horsepower) and `wt` (weight) as predictors. ### Cross-Validation Setup - **5-fold cross-validation**: The dataset is split into **5 different folds**, where 4 folds are used for training and 1 fold is used for testing in each iteration. - **Sample sizes**: The number of observations used in each fold. ### Performance Metrics - **RMSE (Root Mean Squared Error)**: Measures how far off predictions are from actual values. - Lower values indicate better predictive accuracy. - Here, RMSE = **2.70**, meaning the model’s predictions are on average **2.70 mpg off** from actual values. - **R-squared (0.84)**: Indicates how much variance in `mpg` is explained by `hp` and `wt`. - Higher values (closer to 1) indicate better model fit. - Here, **84% of the variance** in `mpg` is explained by the model. - **MAE (Mean Absolute Error)**: Measures the average absolute difference between predicted and actual values. - Like RMSE, lower values indicate better accuracy. - Here, MAE = **2.14**, meaning the model’s average absolute prediction error is **2.14 mpg**. ### How to Interpret These Metrics? | **Metric** | **Interpretation** | |-----------------------------|-------------------------------------------| | **RMSE** | Lower is better; measures average prediction error. | | **R-squared** | Higher is better; shows proportion of variance explained. | | **MAE** | Lower is better; gives an intuitive measure of error. | If RMSE and MAE are **low**, and R² is **high**, the model is performing well.\ If RMSE is **high**, the model has poor predictive accuracy.\ If R² is **low**, the predictors do not explain much of the variance in the outcome. ### What to Do If Performance is Poor? If the model has **high RMSE / low R²**, consider: - **Adding more relevant predictors** to improve explanatory power. - **Transforming variables** (e.g., log, square root) to improve linearity. - **Using a different modeling approach** (e.g., polynomial regression, ridge regression). By checking cross-validation performance, we ensure the model is robust and generalizes well to new data. # Visualize predictions ```{r} ggplot(mtcars, aes(x = mpg, y = predicted_mpg)) + geom_point(color = "blue") + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") + labs(title = "Actual vs. Predicted MPG", x = "Actual MPG", y = "Predicted MPG") + theme_minimal() ``` # Summary -️ **Multiple regression** helps us understand relationships between variables.\ -️ **Predictions** allow us to estimate outcomes for new data.\ -️ **Cross-validation** ensures our model generalizes beyond the dataset.\ -️ **Applications in psychology/neuroscience** include predicting cognitive performance and brain activity. # Deep dive: Checking regression assumptions Before interpreting a regression model, it’s essential to check if it meets the standard **assumptions** of linear regression. If these assumptions are violated, the model’s estimates might be biased or misleading. We are not going to go into detail on this section, but it is included here for your future reference. ## 1. Linearity Assumption The relationship between the predictors (`hp`, `wt`) and the outcome (`mpg`) should be **linear**. ### How to Check: Scatterplots ```{r} # Scatterplots for each predictor vs. mpg par(mfrow = c(1,2)) # Arrange plots side by side plot(mtcars$hp, mtcars$mpg, main = "MPG vs. HP", xlab = "Horsepower", ylab = "MPG", pch = 19) abline(lm(mpg ~ hp, data = mtcars), col = "red") plot(mtcars$wt, mtcars$mpg, main = "MPG vs. Weight", xlab = "Weight", ylab = "MPG", pch = 19) abline(lm(mpg ~ wt, data = mtcars), col = "red") ``` - The red regression line should fit the scatterplot reasonably well.\ - If the relationship is **curved**, consider **polynomial regression** or transformations (e.g., log, square root). ## 2. Normality of Residuals The residuals should be **normally distributed** to ensure valid hypothesis testing. ### How to Check: Histogram & Q-Q Plot ```{r} # Histogram of residuals hist(resid(model), main = "Histogram of Residuals", col = "lightblue", breaks = 10) # Q-Q Plot qqnorm(resid(model)) qqline(resid(model), col = "red") ``` - The histogram should be **bell-shaped**.\ - The Q-Q plot should show points **falling along the red line**.\ - If residuals are **skewed**, consider transformations or robust regression. ## 3. Homoscedasticity (Constant Variance) Residuals should have **constant variance** across predicted values. ### How to Check: Residuals vs. Fitted Plot ```{r} plot(model$fitted.values, resid(model), main = "Residuals vs. Fitted Values", xlab = "Fitted Values", ylab = "Residuals", pch = 19) abline(h = 0, col = "red", lty = 2) ``` - Points should be **randomly scattered** with no clear pattern.\ - If you see a **funnel shape (increasing or decreasing spread)**, consider using **log transformations** or weighted regression. ## 4. No Multicollinearity Predictors should not be **highly correlated** with each other. ### How to Check: Variance Inflation Factor (VIF) ```{r, message=FALSE} # Install if necessary # install.packages("car") library(car) # Compute VIF values vif(model) ``` - **VIF \> 5** indicates **moderate multicollinearity**.\ - **VIF \> 10** suggests **severe multicollinearity**, and one of the predictors should be removed or combined. ## 5. No Autocorrelation of Residuals Residuals should not show patterns or correlations across observations. ### How to Check: Durbin-Watson Test ```{r, message=FALSE} # Install if necessary # install.packages("lmtest") library(lmtest) # Durbin-Watson Test dwtest(model) ``` - **p \> 0.05** suggests no autocorrelation.\ - **p \< 0.05** suggests correlation, meaning observations might not be independent. ## Conclusion: What to Do if Assumptions Are Violated? | Issue | Solution | |--------------------------|----------------------------------------------| | Non-linearity | Try log or polynomial transformations. | | Non-normal residuals | Apply transformations (log, sqrt) or use robust regression. | | Heteroscedasticity | Use weighted regression or log transformations. | | Multicollinearity | Remove highly correlated predictors or combine them. | | Autocorrelation | Consider adding lag variables or using time-series models. | If most assumptions hold, you can confidently interpret your regression results. ------------------------------------------------------------------------