# Load in Data github_url <- "https://raw.githubusercontent.com/jorgeandr3s/heal/main/cscs/public_data/CSCS2025_full_cleaned_deidentified_data_and_metadata.RData" load(url(github_url)) data <- data[complete.cases( data[, c( "LONELY_dejong_emotional_social_loneliness_scale_score", "DEMO_age", "DEMO_gender" )] ), ] ########################################################################################### ########################################################################################### ################################ DATA EXPLORATION ######################################## ########################################################################################### ########################################################################################### #Explore Outcome hist(data$LONELY_dejong_emotional_social_loneliness_scale_score) summary(data$LONELY_dejong_emotional_social_loneliness_scale_score) sd(data$LONELY_dejong_emotional_social_loneliness_scale_score, na.rm = TRUE) #Explore Exposure summary(data$DEMO_age) sd(data$DEMO_age, na.rm = TRUE) table(data$DEMO_gender) #Explore Association plot(LONELY_dejong_emotional_social_loneliness_scale_score ~ DEMO_age, data = data) abline(lm(formula = LONELY_dejong_emotional_social_loneliness_scale_score ~ DEMO_age, data = data)) boxplot(LONELY_dejong_emotional_social_loneliness_scale_score ~ DEMO_gender, data = data) ########################################################################################### ########################################################################################### ################################ LINEAR REGRESSION ####################################### ########################################################################################### ########################################################################################### #Set Up and Run Linear Regression model1 <- lm(formula = LONELY_dejong_emotional_social_loneliness_scale_score ~ DEMO_age, data = data) summary(model1) #Gives you Beta Coefficients, Standard errors, and p-values. #R squared tells you how much variation in outcome the model explains, adjusted R squared penalizes unnecessary variables #F-test tells you whether the model explains something model2 <- lm(formula = LONELY_dejong_emotional_social_loneliness_scale_score ~ DEMO_age + DEMO_gender, data = data) anova(model1, model2) #Compares two models; if p-value is signficiant for nested model, it indicates that model fit is significantly improved #Including an interaction term between age and gender model1_interaction <- lm(formula = LONELY_dejong_emotional_social_loneliness_scale_score ~ DEMO_age*DEMO_gender, data = data) summary(model1_interaction) anova(model2, model1_interaction) confint(model1) #Gives you 95% confidence intervals for the coefficient. ########################################################################################### ########################################################################################### ################################ LINEAR MODEL CHECKING #################################### ########################################################################################### ########################################################################################### # Model checking combines: # 1. Graphical diagnostics, which show the magnitude and pattern of deviations. # 2. Formal statistical tests or quantitative diagnostics, which provide # statistical evidence of deviations. # # These should be interpreted together. # # IMPORTANT: # Formal tests are affected by sample size. In large samples, even very small # departures from model assumptions can produce statistically significant # p-values. Therefore, a statistically significant test does not necessarily # indicate a practically important violation. # # Conversely, graphical diagnostics can help assess whether an identified # deviation is substantial enough to matter for the analysis. # # General principle: # Statistical test = Is there evidence of a deviation? # Plot/diagnostic = What does the deviation look like, and how large is it? # # Independence is different: it is primarily assessed from the study design # and data structure rather than from the standard diagnostic plots. ########################################################################################### # 1. LINEARITY ########################################################################################### # Graphical assessment: # Residuals vs. Fitted plot # # Checks whether the relationship between the predictors and the outcome # is appropriately represented by the linear model. # # Look for: # - Residuals randomly scattered around zero # - No systematic curvature or other pattern # # A curved pattern may suggest that a linear functional form is inadequate. plot(model2, which = 1) # Formal/quantitative assessment: # # There is no single universally required formal test of linearity for # ordinary linear regression. Linearity is primarily assessed graphically # and through the specification of alternative models. # # For example, a quadratic term can be added to assess whether a curved # relationship substantially improves model fit: # # model_quadratic <- lm( # outcome ~ age + I(age^2) + other_predictors, # data = data # ) # # The linear and quadratic models can then be compared using an F-test: # # anova(model2, model_quadratic) # # p > 0.05: # No statistical evidence that the quadratic term improves model fit. # A linear specification is reasonable. # # p < 0.05: # Evidence that the quadratic term improves model fit, suggesting # that the relationship may be nonlinear. # # IMPORTANT: # A large sample can detect very small departures from linearity. # Therefore, statistical evidence of nonlinearity should be considered # alongside the residual plot and the practical magnitude of the deviation. ########################################################################################### # 2. CONSTANT VARIANCE (HOMOSCEDASTICITY) ########################################################################################### # Graphical assessment: # Scale-Location plot # # Checks whether the variability of the residuals is approximately constant # across the range of fitted values. # # Look for: # - Approximately equal vertical spread across fitted values # - No systematic increase or decrease in spread # # A funnel-shaped pattern suggests heteroscedasticity. plot(model2, which = 3) # Formal test: # Breusch-Pagan Test # # Tests whether the variance of the residuals is related to the fitted values. install.packages("lmtest") # Only need to run this once library(lmtest) bptest(model2) # Null hypothesis (H0): # The residual variance is constant. # # Alternative hypothesis (H1): # The residual variance is not constant. # # p > 0.05: # No statistical evidence of heteroscedasticity. # # p < 0.05: # Statistical evidence of heteroscedasticity. # # Interpretation: # Use the p-value to determine whether there is statistical evidence # of non-constant variance, and use the Scale-Location plot to determine # whether the deviation appears substantial or systematic. # # IMPORTANT: # In large samples, the Breusch-Pagan test can detect very small departures # from constant variance. A significant p-value should therefore be interpreted # alongside the diagnostic plot and the practical magnitude of the problem. # Alternative formal test: # Non-Constant Variance Test # # ncvTest() provides another test of non-constant error variance. install.packages("car") # Only need to run this once library(car) ncvTest(model2) # Null hypothesis (H0): # The error variance is constant. # # Alternative hypothesis (H1): # The error variance is not constant. # # p > 0.05: # No statistical evidence of non-constant variance. # # p < 0.05: # Statistical evidence of non-constant variance. # # ncvTest() is an alternative to the Breusch-Pagan test. # You do not generally need to run both. ########################################################################################### # 3. NORMALITY OF RESIDUALS ########################################################################################### # Graphical assessment: # Normal Q-Q plot # # Checks whether the distribution of residuals is approximately normal. # # Look for: # - Points approximately following the reference line # - No substantial systematic departures from the line # # Deviations at the extreme ends of the plot indicate departures in # the tails of the residual distribution. plot(model2, which = 2) # Formal test: # Shapiro-Wilk Test # # Tests whether the residuals are consistent with a normal distribution. summary(model2) set.seed(123) resid_sample <- sample( residuals(model2), 5000 ) shapiro.test(resid_sample) # Null hypothesis (H0): # The residuals are normally distributed. # # Alternative hypothesis (H1): # The residuals are not normally distributed. # # p > 0.05: # No statistical evidence of non-normality. # # p < 0.05: # Statistical evidence that the residuals are not normally distributed. # # Interpretation: # Use the Shapiro-Wilk test to assess statistical evidence of non-normality, # while using the Q-Q plot to assess the magnitude and pattern of the deviation. # # IMPORTANT: # The Shapiro-Wilk test is highly sensitive in large samples. With sufficiently # large samples, very small departures from normality can produce a significant # p-value. Therefore, a significant result does not necessarily indicate that # the model is seriously affected. # # The Q-Q plot is particularly important for determining whether deviations # are substantial, especially in the tails of the distribution. ########################################################################################### # 4. INFLUENTIAL OBSERVATIONS ########################################################################################### # Graphical assessment: # Residuals vs. Leverage plot # # Identifies observations that have unusual combinations of predictor values # and residuals and therefore may have a disproportionate effect on the model. # # Look for: # - Observations with high leverage # - Observations with large residuals # - Observations near or beyond Cook's distance contours plot(model2, which = 5) # Quantitative assessment: # Cook's Distance # # Cook's distance measures how much the estimated regression model would # change if an observation were removed. # # There is no p-value associated with Cook's distance. # Larger values indicate observations that may have greater influence # on the estimated regression coefficients. cooks_d <- cooks.distance(model2) # Plot Cook's distance for every observation plot( cooks_d, type = "h", ylab = "Cook's Distance", xlab = "Observation" ) # Add the 4/n screening threshold abline( h = 4 / nrow(data), lty = 2 ) # Identify observations above the 4/n threshold which(cooks_d > 4 / nrow(data)) # See the largest Cook's distances sort(cooks_d, decreasing = TRUE)[1:10] # Interpretation: # Cook's distance is used as a screening diagnostic rather than a formal # hypothesis test. # # The Residuals vs. Leverage plot helps identify observations that may be # influential, while Cook's distance quantifies their potential influence. # # A common rule of thumb is to investigate observations with: # # Cook's distance > 4 / n # # This is NOT a definitive cutoff for removing observations. # An observation should not be removed solely because it exceeds this threshold. # Investigate the observation, determine why it is influential, and consider # whether it represents a data error, an unusual but valid observation, # or an important part of the population being studied. ########################################################################################### # 5. INDEPENDENCE ########################################################################################### # Independence is primarily assessed from the study design and data structure, # rather than from the standard four regression diagnostic plots. # # Consider whether observations are independent of one another. # # Examples of potential violations include: # - Repeated measurements from the same participant # - Participants clustered within schools, workplaces, hospitals, etc. # - Longitudinal observations # - Spatial or temporal clustering # # If observations are not independent, a different modeling approach may be # required, such as a mixed-effects model or generalized estimating equations. # # Independence cannot generally be established simply because a diagnostic # plot looks acceptable. ########################################################################################### # 6. COLLINEARITY ########################################################################################### # Quantitative assessment: # Variance Inflation Factor (VIF) # # Collinearity occurs when predictors are strongly correlated with one another. # # Why it matters: # High collinearity can make it difficult to distinguish the individual # associations of predictors. It can increase standard errors, produce # unstable coefficients, and make statistically significant associations # become non-significant. install.packages("regclass") # Only need to run this once regclass::VIF(model2) # Interpretation: # # VIF = 1: # No collinearity with the other predictors. # # VIF between 1 and 5: # Generally considered acceptable. # # VIF > 5: # Potentially problematic collinearity. # # VIF > 10: # Often considered serious collinearity. # # These are rules of thumb rather than strict thresholds. # Consider the magnitude of the VIF alongside the variables, # study design, and research context. # # Unlike linearity, constant variance, and normality, collinearity is not # primarily assessed using the standard residual diagnostic plots. # VIF provides the main quantitative assessment. ########################################################################################### # OVERALL INTERPRETATION ########################################################################################### # Model assumptions should not be evaluated using a single diagnostic. # # Graphical diagnostics show the pattern and practical magnitude of potential # problems, while formal tests provide statistical evidence that a deviation # exists. # # In large samples, formal tests have high statistical power and can identify # very small departures from assumptions. Therefore: # # Statistical significance does not necessarily mean practical importance. # # Similarly, a non-significant test does not prove that an assumption is # perfectly satisfied. The appropriate conclusion considers the statistical # test, diagnostic plot, sample size, and substantive context together. # # Independence is primarily evaluated from study design and data structure, # while influential observations and collinearity are evaluated using # quantitative diagnostics. ########################################################################################### ########################################################################################### # OVERALL CONCLUSION: LINEAR REGRESSION ########################################################################################### # The linearity assumption appears reasonable based on the residual diagnostics, # although this should be confirmed by examining the Residuals vs. Fitted plot. # # There is strong statistical evidence of heteroscedasticity from both the # Breusch-Pagan and non-constant variance tests (p < 0.001). Given the large # sample size, these tests may detect relatively small departures from constant # variance; the Scale-Location plot should therefore be used to assess the # practical magnitude of the deviation. This is somewhat subjective, but # the Scale-Location plot shows moderate heteroscedasticity # # There is also statistical evidence of non-normal residuals based on the # Shapiro-Wilk test (p < 0.001). The Q-Q plot should be used to determine whether # the deviation is substantial. With a large sample, minor departures from # normality are expected to produce significant results. # The Q-Q plot shows short-tailed, negatively skewed residuals. # # Cook's distance identifies many observations above the 4/n screening threshold, # but the largest Cook's distance is approximately 0.025. These observations # should be investigated but should not be removed solely because they exceed # the screening threshold. # # Both violations are somewhat typical with a bounded discrete score, so you can # choose to maybe run a logistic regression as a sensitivity test to see how different # your results would be # # VIF values are close to 1, indicating no evidence of problematic # multicollinearity between the predictors. # # Independence cannot be established from these diagnostics and should instead # be evaluated based on the study design and data structure. # # Overall, the linear model shows statistical evidence of heteroscedasticity and # non-normality, but no clear evidence of problematic linearity, multicollinearity, # or severe influential observations. If the linear model is retained, robust # standard errors could be considered to account for non-constant variance. ########################################################################################### ########################################################################################### ############################## LOGISTIC Regression ######################################## ########################################################################################### ########################################################################################### logistic_regression_model1 <- glm( formula = LONELY_dejong_emotional_social_loneliness_scale_score_y_n ~ DEMO_age + DEMO_gender, family = "binomial", data = data ) summary(logistic_regression_model1) ########################################################################################### # Odds Ratios and 95% Confidence Intervals ########################################################################################### # Logistic regression coefficients are expressed as log-odds. # Exponentiating the coefficients gives odds ratios. exp(coef(logistic_regression_model1)) # 95% confidence intervals for odds ratios exp(confint(logistic_regression_model1)) # OR > 1: Higher odds of the outcome # OR < 1: Lower odds of the outcome # OR = 1: No association ########################################################################################### ########################################################################################### ############################## LOGISTIC MODEL CHECKING #################################### ########################################################################################### ########################################################################################### # Logistic regression does NOT assume: # # - Normally distributed residuals # - Constant variance (homoscedasticity) # # Therefore, Shapiro-Wilk, Breusch-Pagan, and NCV tests used for # linear regression are not standard assumption tests for logistic regression. # # The main issues to assess are: # # 1. Linearity of continuous predictors on the logit scale # 2. Multicollinearity # 3. Influential observations # 4. Overall model fit # 5. Model discrimination ########################################################################################### # 1. Linearity of the Logit ########################################################################################### # Logistic regression assumes that continuous predictors have a # linear relationship with the LOGIT (log-odds) of the outcome. # # This does NOT mean that age must have a linear relationship with # the probability of the outcome. # Compare the linear model with a model containing a quadratic age term. logistic_regression_model_quadratic <- glm( formula = LONELY_dejong_emotional_social_loneliness_scale_score_y_n ~ DEMO_age + I(DEMO_age^2) + DEMO_gender, family = "binomial", data = data ) summary(logistic_regression_model_quadratic) # Compare model fit anova( logistic_regression_model1, logistic_regression_model_quadratic, test = "Chisq" ) # p > 0.05: # No statistical evidence that the quadratic term improves model fit. # # p < 0.05: # Evidence that the relationship between age and the logit may be nonlinear. # # A statistically significant nonlinear term should be interpreted # alongside the research question and graphical evidence. ########################################################################################### # 2. Multicollinearity ########################################################################################### # Checks whether predictors are strongly correlated with one another. regclass::VIF(logistic_regression_model1) # VIF = 1: No collinearity # VIF between 1 and 5: Generally acceptable # VIF > 5: Potentially problematic # VIF > 10: Often considered serious # # These are rules of thumb rather than strict cutoffs. ########################################################################################### # 3. Influential Observations: Cook's Distance ########################################################################################### # Cook's distance assesses how much the estimated model changes # when an observation is removed. cooks_d <- cooks.distance(logistic_regression_model1) plot( cooks_d, type = "h", ylab = "Cook's Distance", xlab = "Observation" ) # Common screening rule: # Cook's distance > 4/n abline( h = 4 / nrow(data), lty = 2 ) # Identify observations above the screening threshold which(cooks_d > 4 / nrow(data)) # Examine the largest Cook's distances sort(cooks_d, decreasing = TRUE)[1:10] # IMPORTANT: # A large Cook's distance identifies an observation for investigation. # It does NOT automatically mean the observation should be removed. ########################################################################################### # 4. Overall Model Fit: Likelihood-Ratio Test ########################################################################################### # Compare the fitted model with an intercept-only model. null_model <- glm( LONELY_dejong_emotional_social_loneliness_scale_score_y_n ~ 1, family = "binomial", data = data ) anova( null_model, logistic_regression_model1, test = "Chisq" ) # H0: The predictors do not improve model fit compared with the null model. # H1: The predictors improve model fit. # # p < 0.05: # The model provides statistically significant improvement in fit # compared with the intercept-only model. ########################################################################################### # 5. Hosmer-Lemeshow Goodness-of-Fit Test ########################################################################################### # The Hosmer-Lemeshow test assesses how well a logistic regression model's predicted probabilities correspond to the outcomes that are actually observed. install.packages("ResourceSelection") # Only need to run this once library(ResourceSelection) hoslem.test( data$LONELY_dejong_emotional_social_loneliness_scale_score_y_n, fitted(logistic_regression_model1), g = 10 ) # H0: Observed and expected outcomes do not differ significantly. # H1: Observed and expected outcomes differ significantly. # # p > 0.05: # No statistical evidence of poor model fit. # # p < 0.05: # Evidence of lack of fit. # # IMPORTANT: # The Hosmer-Lemeshow test is sensitive to sample size and to the # number of groups used. Therefore, it should not be the sole # assessment of model fit. ########################################################################################### # 6. Pseudo R-Squared: McFadden's R2 ########################################################################################### # McFadden's R2 provides a measure of model fit based on the # likelihood of the fitted model relative to the null model. # install.packages("pscl") # Only need to run this once library(pscl) pR2(logistic_regression_model1) # IMPORTANT: # McFadden's R2 is NOT interpreted in the same way as R2 in # ordinary linear regression. # # It should not be interpreted as the percentage of variance # explained by the model. ########################################################################################### # 7. ROC Curve and AUC ########################################################################################### # AUC describes the model's ability to discriminate between participants # with and without the outcome. # # AUC = 0.50: # No better than chance discrimination. # # AUC = 1.00: # Perfect discrimination. # # Higher AUC indicates better discrimination. # # AUC should be interpreted as a measure of discrimination, # NOT calibration or causal validity. install.packages("pROC") # Only need to run this once library(pROC) # Generate predicted probabilities predicted_probability <- predict( logistic_regression_model1, type = "response" ) # Generate ROC curve roc_model <- roc( data$LONELY_dejong_emotional_social_loneliness_scale_score_y_n, predicted_probability ) plot( roc_model, xlab = "1 - Specificity", ylab = "Sensitivity" ) # Area Under the Curve (AUC) auc(roc_model) ########################################################################################### # OVERALL CONCLUSION: LOGISTIC REGRESSION ########################################################################################### # The linearity assumption for continuous predictors on the logit scale appears # reasonable. Adding a quadratic age term did not significantly improve model # fit (likelihood-ratio test p = 0.568), providing no statistical evidence that # the association between age and the log-odds of the outcome is nonlinear. # # VIF values are close to 1, indicating no evidence of problematic # multicollinearity between the predictors. # # Cook's distance identifies observations above the 4/n screening threshold, # but the largest Cook's distance is approximately 0.031. These observations # should be investigated but should not be removed solely because they exceed # the screening threshold. # # Independence cannot be established from these diagnostics and should instead # be evaluated based on the study design and data structure. # # The likelihood-ratio test indicates that the model provides a statistically # significant improvement in fit compared with the intercept-only model # (p < 0.001). This is a model-fit assessment rather than an assumption test. # # The Hosmer-Lemeshow test returned NA because the outcome was coded as a factor # in a way that was not compatible with the test. Therefore, the Hosmer-Lemeshow # result cannot be interpreted and should not be used as evidence of good or # poor model fit. # # The AUC of 0.655 indicates that the model has some ability to discriminate # between outcome categories. AUC is a measure of discrimination rather than # an assumption test. # # Overall, the logistic regression diagnostics provide no clear evidence of # violations of the principal model assumptions assessed here. The linearity # of age on the logit scale and multicollinearity appear acceptable, and no # obviously dominant influential observations were identified. Independence # should be justified from the study design.