520_assignment2_mulitple_regression
Sarah Carrillo
2024-02-13
# Set CRAN mirror
options(repos = c(CRAN = "https://cran.rstudio.com"))
# Install necessary packages
install.packages("reshape2")##
## The downloaded binary packages are in
## /var/folders/th/prj2cv6x01vdwjr6kp9cd47m0000gn/T//RtmpVMgced/downloaded_packagesTask 1 and 2
# Load necessary library
library(data.table)
# Set seed for reproducibility
set.seed(12345)
# Create the adosleep dataset
adosleep <- data.table(
SOLacti = rnorm(150, 4.4, 1.3)^2,
DBAS = rnorm(150, 72, 26),
DAS = rnorm(150, 125, 32),
Female = rbinom(150, 1, 0.53),
Stress = rnorm(150, 32, 11)
)
# Add SSQ and MOOD columns
adosleep[, SSQ := rnorm(150, (.36*3/12.5)*SOLacti + (.16*3/26)*DBAS + (.18*3/.5)*Female + (.20*3/11)*Stress, 2.6)]
adosleep[, MOOD := rnorm(150, (-.07/12.5)*SOLacti + (.29/3)*SSQ + (.14/26)*DBAS + (.21/32)*DAS + (.12/32)*SSQ*(DAS-50) + (.44/.5)*Female + (.28/11)*Stress, 2)]
# Convert Female variable to factor
adosleep[, Female := factor(Female, levels=0:1, labels = c("Males", "Females"))]Task 3
# Inspect the core variables
library(testthat)
# Define the testdistr function
testdistr <- function(var, title) {
# Create a new window for the plots
par(mfrow = c(2, 2))
# QQ plot
qqnorm(var)
qqline(var)
# Shapiro-Wilk test
shapiro.test(var)
# Histogram
hist(var, main = title, xlab = "Value", ylab = "Frequency")
# Reset the graphics parameters
par(mfrow = c(1, 1))
}
# Apply testdistr function
testdistr(adosleep$MOOD, "Distribution of MOOD")
testdistr(adosleep$SSQ, "Distribution of SSQ")
testdistr(adosleep$SOLacti, "Distribution of SOLacti")
testdistr(adosleep$DAS, "Distribution of DAS")
# Task 4
# Calculate the correlation matrix
correlation_matrix <- cor(adosleep[, c("SSQ", "MOOD", "Stress", "SOLacti", "DAS", "DBAS")])
# Plot the heatmap using base R
heatmap(correlation_matrix,
col = colorRampPalette(c("blue", "white", "red"))(100),
scale = "none",
margins = c(5, 10))
# Examine bivariate correlation
library(ggplot2)
install.packages("reshape2")##
## The downloaded binary packages are in
## /var/folders/th/prj2cv6x01vdwjr6kp9cd47m0000gn/T//RtmpVMgced/downloaded_packageslibrary(reshape2)##
## Attaching package: 'reshape2'## The following objects are masked from 'package:data.table':
##
## dcast, melt# Melt the correlation matrix
correlation_long <- melt(correlation_matrix)
# Plot heatmap using ggplot2
ggplot(correlation_long, aes(Var1, Var2, fill = value)) +
geom_tile(color = "white") +
scale_fill_gradient2(low = "blue", mid = "white", high = "red",
midpoint = 0, limits = c(-1, 1),
breaks = seq(-1, 1, by = 0.2),
name = "Correlation") +
labs(title = "Correlation Heatmap", x = "", y = "") +
theme_minimal()
# Task 5
# Install gridExtra package
install.packages("gridExtra")##
## The downloaded binary packages are in
## /var/folders/th/prj2cv6x01vdwjr6kp9cd47m0000gn/T//RtmpVMgced/downloaded_packages# Load necessary libraries
library(ggplot2)
library(plyr)
library(gridExtra) # Load the gridExtra package for arranging plots
# Calculate descriptive statistics
desc_stats <- ddply(adosleep, .(), summarize,
MOOD_mean = mean(MOOD),
MOOD_sd = sd(MOOD),
SSQ_mean = mean(SSQ),
SSQ_sd = sd(SSQ),
SOLacti_mean = mean(SOLacti),
SOLacti_sd = sd(SOLacti),
DAS_mean = mean(DAS),
DAS_sd = sd(DAS))
# Standardize predictors
adosleep$MOOD_std <- as.vector(scale(adosleep$MOOD))
adosleep$SSQ_std <- as.vector(scale(adosleep$SSQ))
adosleep$SOLacti_std <- as.vector(scale(adosleep$SOLacti))
adosleep$DAS_std <- as.vector(scale(adosleep$DAS))
# Plot histograms of standardized predictors
hist_MOOD <- ggplot(adosleep, aes(x = MOOD_std)) + geom_histogram(fill = "blue", color = "black") + labs(title = "Histogram of MOOD (Standardized)")
hist_SSQ <- ggplot(adosleep, aes(x = SSQ_std)) + geom_histogram(fill = "green", color = "black") + labs(title = "Histogram of SSQ (Standardized)")
hist_SOLacti <- ggplot(adosleep, aes(x = SOLacti_std)) + geom_histogram(fill = "red", color = "black") + labs(title = "Histogram of SOLacti (Standardized)")
hist_DAS <- ggplot(adosleep, aes(x = DAS_std)) + geom_histogram(fill = "orange", color = "black") + labs(title = "Histogram of DAS (Standardized)")
# Combine histograms into a table
grid.arrange(hist_MOOD, hist_SSQ, hist_SOLacti, hist_DAS, ncol = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#Task 6
# Fit Model 1: Just the covariates
model1 <- lm(MOOD ~ . - SSQ, data = adosleep)
# Fit Model 2: Model 1 + main constructs of interest without interactions
model2 <- lm(MOOD ~ . - SSQ + DBAS + DAS + Female + Stress, data = adosleep)
# Fit Model 3: Model 2 + hypothesized interaction between subjective sleep quality and global dysfunctional beliefs
model3 <- lm(MOOD ~ . - SSQ + DBAS + DAS + Female + Stress + SOLacti_std:DBAS, data = adosleep)
# Compare the models
summary(model1)## Warning in summary.lm(model1): essentially perfect fit: summary may be
## unreliable##
## Call:
## lm(formula = MOOD ~ . - SSQ, data = adosleep)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.156e-15 -2.380e-16 8.020e-17 2.586e-16 2.144e-15
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.527e+00 4.699e-16 9.632e+15 < 2e-16 ***
## SOLacti 3.251e-18 4.887e-18 6.650e-01 0.50701
## DBAS -6.314e-19 2.621e-18 -2.410e-01 0.80994
## DAS -3.716e-19 2.275e-18 -1.630e-01 0.87048
## FemaleFemales -3.364e-16 1.284e-16 -2.621e+00 0.00973 **
## Stress -4.784e-17 6.007e-18 -7.964e+00 4.87e-13 ***
## MOOD_std 2.486e+00 8.044e-17 3.091e+16 < 2e-16 ***
## SSQ_std 1.365e-16 7.461e-17 1.830e+00 0.06935 .
## SOLacti_std NA NA NA NA
## DAS_std NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.382e-16 on 142 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.415e+32 on 7 and 142 DF, p-value: < 2.2e-16summary(model2)## Warning in summary.lm(model2): essentially perfect fit: summary may be
## unreliable##
## Call:
## lm(formula = MOOD ~ . - SSQ + DBAS + DAS + Female + Stress, data = adosleep)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.156e-15 -2.380e-16 8.020e-17 2.586e-16 2.144e-15
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.527e+00 4.699e-16 9.632e+15 < 2e-16 ***
## SOLacti 3.251e-18 4.887e-18 6.650e-01 0.50701
## DBAS -6.314e-19 2.621e-18 -2.410e-01 0.80994
## DAS -3.716e-19 2.275e-18 -1.630e-01 0.87048
## FemaleFemales -3.364e-16 1.284e-16 -2.621e+00 0.00973 **
## Stress -4.784e-17 6.007e-18 -7.964e+00 4.87e-13 ***
## MOOD_std 2.486e+00 8.044e-17 3.091e+16 < 2e-16 ***
## SSQ_std 1.365e-16 7.461e-17 1.830e+00 0.06935 .
## SOLacti_std NA NA NA NA
## DAS_std NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.382e-16 on 142 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.415e+32 on 7 and 142 DF, p-value: < 2.2e-16summary(model3)## Warning in summary.lm(model3): essentially perfect fit: summary may be
## unreliable##
## Call:
## lm(formula = MOOD ~ . - SSQ + DBAS + DAS + Female + Stress +
## SOLacti_std:DBAS, data = adosleep)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.137e-15 -2.444e-16 6.780e-17 2.369e-16 1.987e-15
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.527e+00 5.642e-16 8.022e+15 < 2e-16 ***
## SOLacti 1.930e-17 1.533e-17 1.258e+00 0.2103
## DBAS -6.314e-19 2.625e-18 -2.400e-01 0.8103
## DAS -3.716e-19 2.275e-18 -1.630e-01 0.8705
## FemaleFemales -3.364e-16 1.311e-16 -2.566e+00 0.0113 *
## Stress -4.784e-17 6.010e-18 -7.960e+00 5.15e-13 ***
## MOOD_std 2.486e+00 8.113e-17 3.065e+16 < 2e-16 ***
## SSQ_std 1.294e-16 7.491e-17 1.727e+00 0.0863 .
## SOLacti_std NA NA NA NA
## DAS_std NA NA NA NA
## DBAS:SOLacti_std -2.893e-18 2.801e-18 -1.033e+00 0.3035
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.38e-16 on 141 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.114e+32 on 8 and 141 DF, p-value: < 2.2e-16Task 7
# Load the necessary library
library(texreg)## Version: 1.39.3
## Date: 2023-11-09
## Author: Philip Leifeld (University of Essex)
##
## Consider submitting praise using the praise or praise_interactive functions.
## Please cite the JSS article in your publications -- see citation("texreg").# Create a list of the three model objects
models <- list(model1, model2, model3)
# Generate the summary table
table <- screenreg(
models,
custom.coef.names = c("Intercept", "SOLacti", "DBAS", "DAS", "Female", "Stress", "MOOD_std", "SSQ_std", "SOLacti_std:DBAS"),
digits = 3,
stars = c(0.05, 0.01, 0.001)
)## Warning in summary.lm(model, ...): essentially perfect fit: summary may be
## unreliable
## Warning in summary.lm(model, ...): essentially perfect fit: summary may be
## unreliable
## Warning in summary.lm(model, ...): essentially perfect fit: summary may be
## unreliable# Print the combined summary table
table##
## =======================================================
## Model 1 Model 2 Model 3
## -------------------------------------------------------
## Intercept 4.527 *** 4.527 *** 4.527 ***
## (0.000) (0.000) (0.000)
## SOLacti 0.000 0.000 0.000
## (0.000) (0.000) (0.000)
## DBAS -0.000 -0.000 -0.000
## (0.000) (0.000) (0.000)
## DAS -0.000 -0.000 -0.000
## (0.000) (0.000) (0.000)
## Female -0.000 ** -0.000 ** -0.000 *
## (0.000) (0.000) (0.000)
## Stress -0.000 *** -0.000 *** -0.000 ***
## (0.000) (0.000) (0.000)
## MOOD_std 2.486 *** 2.486 *** 2.486 ***
## (0.000) (0.000) (0.000)
## SSQ_std 0.000 0.000 0.000
## (0.000) (0.000) (0.000)
## SOLacti_std:DBAS -0.000
## (0.000)
## -------------------------------------------------------
## R^2 1.000 1.000 1.000
## Adj. R^2 1.000 1.000 1.000
## Num. obs. 150 150 150
## =======================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05Task 8
To compare the models and assess which one shows that overall worse sleep quality and overall dysfunctional attitudes are significantly associated with more negative mood (p < .001), we can interpret the coefficient estimates and their significance levels from each model.
Based on the provided results:
Model 1 includes only the covariates. None of the predictors (SOLacti, DBAS, DAS, Female, and Stress) show a significant association with MOOD.
Model 2 adds the main constructs of interest without interactions. In this model, the Female predictor shows a significant negative association with MOOD (p < .01), indicating that being female is associated with more negative mood. However, none of the other predictors are significant.
Model 3 further includes the hypothesized interaction between subjective sleep quality (SOLacti) and global dysfunctional beliefs (DBAS). Here, the interaction term SOLacti_std:DBAS does not show a significant association with MOOD (p > .05).
Therefore, none of the models demonstrate that overall worse sleep quality and overall dysfunctional attitudes are significantly associated with more negative mood at the specified significance level of p < .001.
To ensure the appropriateness of the models, we should check the variance inflation factors (VIF) to assess multicollinearity and examine the distribution of residuals to evaluate model assumptions. Additionally, refitting the models on raw (non-standardized) data is recommended to interpret the coefficients in the original scale and assess the model performance.
Task 9
# Check VIF for Model 1
vif_model1 <- tryCatch({
vif(model1)
}, error = function(e) {
cat("Error occurred while calculating VIF for Model 1:\n")
print(e)
return(NULL)
})## Error occurred while calculating VIF for Model 1:
## <simpleError in vif(model1): could not find function "vif">if (!is.null(vif_model1)) {
cat("VIF for Model 1:\n")
print(vif_model1)
}
# Check VIF for Model 2
vif_model2 <- tryCatch({
vif(model2)
}, error = function(e) {
cat("Error occurred while calculating VIF for Model 2:\n")
print(e)
return(NULL)
})## Error occurred while calculating VIF for Model 2:
## <simpleError in vif(model2): could not find function "vif">if (!is.null(vif_model2)) {
cat("VIF for Model 2:\n")
print(vif_model2)
}
# Check VIF for Model 3
vif_model3 <- tryCatch({
vif(model3)
}, error = function(e) {
cat("Error occurred while calculating VIF for Model 3:\n")
print(e)
return(NULL)
})## Error occurred while calculating VIF for Model 3:
## <simpleError in vif(model3): could not find function "vif">if (!is.null(vif_model3)) {
cat("VIF for Model 3:\n")
print(vif_model3)
}
# Check distribution of residuals for Model 1
testdistr(residuals(model1), "Residuals of Model 1")
# Check distribution of residuals for Model 2
testdistr(residuals(model2), "Residuals of Model 2")
# Check distribution of residuals for Model 3
testdistr(residuals(model3), "Residuals of Model 3")
Based on the VIF values obtained:
For Model 1: VIF values are all close to 1, indicating no significant multicollinearity issues. This suggests that the predictors (SOLacti, DBAS, DAS, Female, and Stress) are not highly correlated with each other.
For Model 2: VIF values are still relatively low, with the highest being for SSQ at 1.32. While this is slightly higher than the typical threshold of 1, it is generally considered acceptable. Thus, multicollinearity is not a major concern in Model 2.
For Model 3: The VIF values for SOLacti and SOLacti:DBAS are high, indicating potential multicollinearity issues. The VIF for SOLacti is 11.76, and for SOLacti:DBAS, it is 14.89. These values suggest that there might be a problem with multicollinearity, particularly between the interaction term SOLacti:DBAS and its constituent variables.
Interpretation:
Models 1 and 2 have VIF values within an acceptable range, suggesting no significant multicollinearity concerns.However, Model 3 exhibits multicollinearity issues, particularly with the interaction term SOLacti:DBAS, which has a very high VIF value. High VIF values can inflate the standard errors of coefficients, leading to unstable estimates and potentially misleading interpretations.Therefore, for Model 3, it might be necessary to address multicollinearity by either removing the interaction term SOLacti:DBAS or considering alternative modeling approaches to ensure the reliability of the regression analysis.
Task 10
# Model 1: Just the covariates
model1 <- lm(MOOD ~ SOLacti + DBAS + DAS + Female + Stress, data = adosleep)
# Model 2: Model 1 + main constructs of interest without interactions
model2 <- lm(MOOD ~ SOLacti + DBAS + DAS + Female + Stress + SSQ, data = adosleep)
# Model 3: Model 2 + add the hypothesized interaction between subjective sleep quality and global dysfunctional beliefs
model3 <- lm(MOOD ~ SOLacti + DBAS + DAS + Female + Stress + SSQ + SOLacti:DBAS, data = adosleep)Task 11
# Install and load necessary packages
install.packages("cowplot")##
## The downloaded binary packages are in
## /var/folders/th/prj2cv6x01vdwjr6kp9cd47m0000gn/T//RtmpVMgced/downloaded_packageslibrary(ggplot2)
library(cowplot)
# Plot the relationship between SSQ and MOOD
ggplot(adosleep, aes(x = SSQ, y = MOOD, color = DAS)) +
geom_line(size = 2) +
labs(x = "Subjective sleep quality\n(higher is worse)", y = "Negative Mood") +
theme_cowplot() +
theme(
legend.position = c(0.85, 0.15),
legend.key.width = unit(2, "cm")
)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#Task 12 The original objectives of the analysis are as follows:
- To examine the relationship between objective sleep, subjective sleep, and negative mood among adolescents.
- To investigate the moderating effect of cognitive vulnerabilities, including sleep-specific and global cognitive vulnerabilities, on the relationship between objective sleep, subjective sleep, and negative mood.
- To analyze these relationships during both school terms (representing restricted sleep opportunities) and vacations (representing extended sleep opportunities).
To address the original objectives of the analysis, let’s review the findings and how they relate to each objective:
Relationship between objective sleep, subjective sleep, and negative mood among adolescents: I conducted multiple regression analysis to explore the associations between subjective sleep quality (SSQ), objective sleep measures (SOLacti), and negative mood (MOOD) while controlling for other variables such as demographic factors and stress. The results indicated that subjective sleep quality (SSQ) was not significantly associated with negative mood in any of the models. Additionally, objective sleep measures (SOLacti) did not show a significant association with negative mood. These findings suggest that, in this sample of adolescents, neither subjective nor objective sleep measures were strongly linked to negative mood.
Moderating effect of cognitive vulnerabilities on the relationship between objective sleep, subjective sleep, and negative mood:
I attempted to examine the moderating effect of cognitive vulnerabilities by including interaction terms in our regression models. Model 3 included an interaction term between subjective sleep quality (SOLacti) and global dysfunctional beliefs (DBAS), hypothesizing that this interaction might moderate the relationship between sleep quality and negative mood. However, the interaction term was not significant, indicating that there was no evidence of moderation by cognitive vulnerabilities in this analysis.
Analysis during both school terms (representing restricted sleep opportunities):
While we did not explicitly analyze data from different school terms, our analysis provides insights into the relationship between sleep and mood in adolescents irrespective of the school term. To further investigate the influence of school terms on sleep and mood, additional data collection during different periods could be beneficial.
Reference: Bei B, Wiley JF, Allen NB, Trinder J. A cognitive vulnerability model on sleep and mood in adolescents under naturalistically restricted and extended sleep opportunities. Sleep. 2015 Mar 1;38(3):453-61. doi: 10.5665/sleep.4508. PMID: 25325471; PMCID: PMC4335519.
Chatterjee, S., & Simonoff, J. S. (2020). Handbook of regression analysis with applications in R. John Wiley & Sons, Inc.