Correlation and simple Regression final
Mackie Dullen
2023-07-15
1 Loading Libraries
library(psych) # for the describe() command and the corr.test() command
library(apaTables) # to create our correlation table
library(kableExtra) # to create our correlation table
library(broom) # for the augment() command
library(ggplot2) # to visualize our results##
## Attaching package: 'ggplot2'## The following objects are masked from 'package:psych':
##
## %+%, alpha2 Importing Data
# import the dataset you cleaned previously
# this will be the dataset you'll use throughout the rest of the semester
# use ARC data downloaded previous for lab
d <- read.csv(file="Data/arc_data_final.csv", header=T)3 Correlation: State Your Hypothesis
We predict that stress (measured by the PSS-4), depression symptoms (measured by the PHQ-9), well-being (measured by swemws), and social support measure by percieved Social Support scale) will all be correlated with each other. Furthermore, we predict that depression symptoms will be higher in participants who have low social support.
4 Check Your Variables
# you only need to check the variables you're using in the current analysis
# although you checked them previously, it's always a good idea to look them over again and be sure that everything is correct
str(d)## 'data.frame': 2073 obs. of 41 variables:
## $ X : int 1 20 30 31 32 33 48 49 57 58 ...
## $ gender : chr "female" "male" "female" "female" ...
## $ trans : chr "no" "no" "no" "no" ...
## $ sexual_orientation : chr "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" ...
## $ ethnicity : chr "White - British, Irish, other" "White - British, Irish, other" "White - British, Irish, other" "White - British, Irish, other" ...
## $ relationship_status : chr "In a relationship/married and cohabiting" "Prefer not to say" "Prefer not to say" "In a relationship/married and cohabiting" ...
## $ age : chr NA "1 under 18" "1 under 18" "4 between 36 and 45" ...
## $ urban_rural : chr "city" "city" "city" "town" ...
## $ income : chr "3 high" NA NA "2 middle" ...
## $ education : chr "6 graduate degree or higher" "prefer not to say" "2 equivalent to high school completion" "5 undergraduate degree" ...
## $ employment : chr "3 employed" "1 high school equivalent" "1 high school equivalent" "3 employed" ...
## $ treatment : chr "no psychological disorders" "in treatment" "not in treatment" "no psychological disorders" ...
## $ health : chr "something else or not applicable" "something else or not applicable" "something else or not applicable" "two conditions" ...
## $ mhealth : chr "none or NA" "anxiety disorder" "none or NA" "none or NA" ...
## $ sleep_hours : chr "3 7-8 hours" "2 5-6 hours" "3 7-8 hours" "2 5-6 hours" ...
## $ exercise : num 0 2 3 1.5 NA 1 NA 2 2 1.7 ...
## $ pet : chr "cat" "cat" "dog" "no pets" ...
## $ covid_pos : int 0 0 0 0 0 0 0 0 0 0 ...
## $ covid_neg : int 0 0 0 0 0 0 0 0 0 0 ...
## $ big5_open : num 5.33 5.33 5 6 NA ...
## $ big5_con : num 6 3.33 5.33 5.67 NA ...
## $ big5_agr : num 4.33 4.33 6.67 4.67 NA ...
## $ big5_neu : num 6 6.67 4 4 NA ...
## $ big5_ext : num 2 1.67 6 5 NA ...
## $ pswq : num 4.94 3.36 1.86 3.94 NA ...
## $ iou : num 3.19 4 1.59 3.37 NA ...
## $ mfq_26 : num 4.2 3.35 4.65 4.65 NA 4.5 NA 4.3 5.25 4.45 ...
## $ mfq_state : num 3.62 3 5.88 4 NA ...
## $ rse : num 2.3 1.6 3.9 1.7 NA 3.9 NA 2.4 1.8 NA ...
## $ school_covid_support: num NA NA NA NA NA NA NA NA NA NA ...
## $ school_att : num NA NA NA NA NA NA NA NA NA NA ...
## $ pas_covid : num 3.22 4.56 3.33 4.22 NA ...
## $ pss : num 3.25 3.75 1 3.25 NA 2 NA 2 4 1.25 ...
## $ phq : num 1.33 3.33 1 2.33 NA ...
## $ gad : num 1.86 3.86 1.14 2 NA ...
## $ edeq12 : num 1.58 1.83 1 1.67 NA ...
## $ brs : num NA NA NA NA NA NA NA NA NA NA ...
## $ swemws : num 2.86 2.29 4.29 3.29 NA ...
## $ isolation_a : num 2.25 NA NA 2.5 NA 1.75 NA 2 1.25 NA ...
## $ isolation_c : num NA 3.5 1 NA NA NA NA NA NA 1 ...
## $ support : num 2.5 2.17 5 2.5 NA ...# we're going to create a fake variable for this lab, so that it has four variables and mirrors the homework assignment. SKIP THIS STEP FOR THE HOMEWORK
d$fake <- (d$pss*d$phq)/d$rse
# since we're focusing on our continuous variables, we're going to subset them into their own dataframe. this will make some stuff we're doing later easier.
cont <- subset(d, select=c(pss, phq, swemws , support))
# you can use the describe() command on an entire dataframe (d) or just on a single variable (d$pss)
describe(cont)## vars n mean sd median trimmed mad min max range skew kurtosis
## pss 1 1356 2.94 0.94 3.00 2.94 1.11 1 5 4 0.07 -0.73
## phq 2 1312 2.07 0.86 1.89 1.99 0.99 1 4 3 0.63 -0.69
## swemws 3 1235 3.16 0.84 3.14 3.18 0.85 1 5 4 -0.28 -0.26
## support 4 1298 3.57 0.95 3.67 3.62 0.99 1 5 4 -0.42 -0.56
## se
## pss 0.03
## phq 0.02
## swemws 0.02
## support 0.03# our fake variable has high kurtosis, which I'll ignore. you don't need to discuss univariate normality in the results write-ups for the labs/homework, but you will need to discuss it in your final manuscript
# also use histograms to examine your continuous variables
hist(d$pss)
hist(d$phq)
hist(d$support)
hist(d$swemws)
# last, use scatterplots to examine your continuous variables together
plot(d$pss, d$phq)
plot(d$pss, d$support)
plot(d$pss, d$swemws)
plot(d$phq, d$support)
plot(d$phq, d$swemws)
plot(d$swemws,d$support)
5 Correlation: Check Your Assumptions
5.1 Pearson’s Correlation Coefficient Assumptions
- Should have two measurements for each participant
- Variables should be continuous and normally distributed
- Outliers should be identified and removed
- Relationship between the variables should be linear
5.1.1 Checking for Outliers
Note: You are not required to screen out outliers or take any action based on what you see here. This is something you will check and then discuss in your write-up.
d$pss_std <- scale(d$pss, center=T, scale=T)
hist(d$pss_std)
sum(d$pss_std < -3 | d$pss_std > 3)## [1] NAd$phq_std <- scale(d$phq, center=T, scale=T)
hist(d$phq_std)
sum(d$phq_std < -3 | d$phq_std > 3)## [1] NAd$support_std <- scale(d$support, center=T, scale=T)
hist(d$support_std)
sum(d$support_std < -3 | d$support_std > 3)## [1] NAd$swemws_std <- scale(d$swemws, center=T, scale=T)
hist(d$swemws_std)
sum(d$swemws_std < -3 | d$swemws_std > 3)## [1] NA5.2 Issues with My Data
All of my variables meet all of the assumptions of Pearson’s correlation coefficient. The kurtosis for each variable are mostly normally distributed and continous with 0 presenting outlier, however, the PSS adn PHQ variables a kurtosis that is lower than the normal disrtribution.
6 Correlation: Create a Correlation Matrix
corr_output_m <- corr.test(cont)7 Correlation: View Test Output
corr_output_m## Call:corr.test(x = cont)
## Correlation matrix
## pss phq swemws support
## pss 1.00 0.75 -0.76 -0.47
## phq 0.75 1.00 -0.76 -0.51
## swemws -0.76 -0.76 1.00 0.59
## support -0.47 -0.51 0.59 1.00
## Sample Size
## pss phq swemws support
## pss 1356 1306 1230 1292
## phq 1306 1312 1205 1269
## swemws 1230 1205 1235 1218
## support 1292 1269 1218 1298
## Probability values (Entries above the diagonal are adjusted for multiple tests.)
## pss phq swemws support
## pss 0 0 0 0
## phq 0 0 0 0
## swemws 0 0 0 0
## support 0 0 0 0
##
## To see confidence intervals of the correlations, print with the short=FALSE option8 Correlation: Write Up Results
To test our hypothesis that stress (measured by the PSS-4), depression symptoms (measured by the PHQ-9), well-being (measured by swemws), and social support measure by percieved Social Support scale) will all be correlated with each other. Furthermore, we predict that depression symptoms will be higher in participants who have low social support. Most of our data met the assumptions of the test, with all variables meeting the standards of normality and no outliers.
As predicted, we found that all four variables were not significantly correlated (all ps < .001). The effect sizes of all correlations were large (rs > .5; Cohen, 1988). This test however did not support our second hypothesis,that depression symptoms will be higher in participants who have low social support, as can be seen by the correlation coefficients reported in Table 1.
The relationship between stress and depression symptoms was positive.The relationship between stress and well-being was negative.The relationship between stress and social support was negative. The relationship between depression symptoms and well-being was negative and was also negtaive with social support.The relationship between well-being and social support was positive.
| Variable | M | SD | 1 | 2 | 3 |
|---|---|---|---|---|---|
| Perceived stress (PSS-4) | 2.94 | 0.94 | |||
| Depression symptoms (PHQ-9) | 2.07 | 0.86 | .75** | ||
| [.72, .77] | |||||
| Well-being(swemws) | 3.16 | 0.84 | -.76** | -.76** | |
| [-.78, -.74] | [-.78, -.73] | ||||
| Social Support (SSS) | 3.57 | 0.95 | -.47** | -.51** | .59** |
| [-.51, -.43] | [-.55, -.47] | [.55, .63] | |||
| Note: | |||||
| M and SD are used to represent mean and standard deviation, respectively. Values in square brackets indicate the 95% confidence interval. The confidence interval is a plausible range of population correlations that could have caused the sample correlation. | |||||
| * indicates p < .05 | |||||
| ** indicates p < .01. |
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.
9 Regression: State Your Hypothesis
We hypothesize that stree (measured by the PSS) will significantly predict depressive symptoms(measured by the PHQ), and that the relationship will be positive.
10 Regression: Check Your Variables
# you only need to check the variables you're using in the current analysis
# although you checked them previously, it's always a good idea to look them over again and be sure that everything is correct
str(d)## 'data.frame': 2073 obs. of 46 variables:
## $ X : int 1 20 30 31 32 33 48 49 57 58 ...
## $ gender : chr "female" "male" "female" "female" ...
## $ trans : chr "no" "no" "no" "no" ...
## $ sexual_orientation : chr "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" ...
## $ ethnicity : chr "White - British, Irish, other" "White - British, Irish, other" "White - British, Irish, other" "White - British, Irish, other" ...
## $ relationship_status : chr "In a relationship/married and cohabiting" "Prefer not to say" "Prefer not to say" "In a relationship/married and cohabiting" ...
## $ age : chr NA "1 under 18" "1 under 18" "4 between 36 and 45" ...
## $ urban_rural : chr "city" "city" "city" "town" ...
## $ income : chr "3 high" NA NA "2 middle" ...
## $ education : chr "6 graduate degree or higher" "prefer not to say" "2 equivalent to high school completion" "5 undergraduate degree" ...
## $ employment : chr "3 employed" "1 high school equivalent" "1 high school equivalent" "3 employed" ...
## $ treatment : chr "no psychological disorders" "in treatment" "not in treatment" "no psychological disorders" ...
## $ health : chr "something else or not applicable" "something else or not applicable" "something else or not applicable" "two conditions" ...
## $ mhealth : chr "none or NA" "anxiety disorder" "none or NA" "none or NA" ...
## $ sleep_hours : chr "3 7-8 hours" "2 5-6 hours" "3 7-8 hours" "2 5-6 hours" ...
## $ exercise : num 0 2 3 1.5 NA 1 NA 2 2 1.7 ...
## $ pet : chr "cat" "cat" "dog" "no pets" ...
## $ covid_pos : int 0 0 0 0 0 0 0 0 0 0 ...
## $ covid_neg : int 0 0 0 0 0 0 0 0 0 0 ...
## $ big5_open : num 5.33 5.33 5 6 NA ...
## $ big5_con : num 6 3.33 5.33 5.67 NA ...
## $ big5_agr : num 4.33 4.33 6.67 4.67 NA ...
## $ big5_neu : num 6 6.67 4 4 NA ...
## $ big5_ext : num 2 1.67 6 5 NA ...
## $ pswq : num 4.94 3.36 1.86 3.94 NA ...
## $ iou : num 3.19 4 1.59 3.37 NA ...
## $ mfq_26 : num 4.2 3.35 4.65 4.65 NA 4.5 NA 4.3 5.25 4.45 ...
## $ mfq_state : num 3.62 3 5.88 4 NA ...
## $ rse : num 2.3 1.6 3.9 1.7 NA 3.9 NA 2.4 1.8 NA ...
## $ school_covid_support: num NA NA NA NA NA NA NA NA NA NA ...
## $ school_att : num NA NA NA NA NA NA NA NA NA NA ...
## $ pas_covid : num 3.22 4.56 3.33 4.22 NA ...
## $ pss : num 3.25 3.75 1 3.25 NA 2 NA 2 4 1.25 ...
## $ phq : num 1.33 3.33 1 2.33 NA ...
## $ gad : num 1.86 3.86 1.14 2 NA ...
## $ edeq12 : num 1.58 1.83 1 1.67 NA ...
## $ brs : num NA NA NA NA NA NA NA NA NA NA ...
## $ swemws : num 2.86 2.29 4.29 3.29 NA ...
## $ isolation_a : num 2.25 NA NA 2.5 NA 1.75 NA 2 1.25 NA ...
## $ isolation_c : num NA 3.5 1 NA NA NA NA NA NA 1 ...
## $ support : num 2.5 2.17 5 2.5 NA ...
## $ fake : num 1.884 7.812 0.256 4.461 NA ...
## $ pss_std : num [1:2073, 1] 0.325 0.855 -2.063 0.325 NA ...
## ..- attr(*, "scaled:center")= num 2.94
## ..- attr(*, "scaled:scale")= num 0.942
## $ phq_std : num [1:2073, 1] -0.858 1.458 -1.244 0.3 NA ...
## ..- attr(*, "scaled:center")= num 2.07
## ..- attr(*, "scaled:scale")= num 0.864
## $ support_std : num [1:2073, 1] -1.12 -1.47 1.51 -1.12 NA ...
## ..- attr(*, "scaled:center")= num 3.57
## ..- attr(*, "scaled:scale")= num 0.95
## $ swemws_std : num [1:2073, 1] -0.355 -1.033 1.34 0.153 NA ...
## ..- attr(*, "scaled:center")= num 3.16
## ..- attr(*, "scaled:scale")= num 0.843# you can use the describe() command on an entire dataframe (d) or just on a single variable
describe(d)## vars n mean sd median trimmed mad min
## X 1 2073 4982.03 2586.74 5393.00 5104.65 3106.05 1.00
## gender* 2 1885 1.44 0.87 1.00 1.27 0.00 1.00
## trans* 3 1885 1.12 0.41 1.00 1.00 0.00 1.00
## sexual_orientation* 4 1885 3.82 1.06 4.00 3.83 0.00 1.00
## ethnicity* 5 1872 7.59 2.80 9.00 8.24 0.00 1.00
## relationship_status* 6 1870 3.68 1.69 5.00 3.85 0.00 1.00
## age* 7 1410 2.10 1.64 1.00 1.88 0.00 1.00
## urban_rural* 8 1836 2.76 1.11 3.00 2.82 1.48 1.00
## income* 9 549 2.42 0.94 2.00 2.41 1.48 1.00
## education* 10 1882 3.07 2.02 2.00 2.86 1.48 1.00
## employment* 11 1885 1.81 1.24 1.00 1.59 0.00 1.00
## treatment* 12 1622 2.76 1.12 3.00 2.61 1.48 1.00
## health* 13 1716 8.60 1.37 9.00 8.91 0.00 1.00
## mhealth* 14 2073 4.63 1.39 5.00 4.88 0.00 1.00
## sleep_hours* 15 1320 2.90 0.97 3.00 2.91 1.48 1.00
## exercise 16 1153 291.55 9816.59 1.75 1.90 1.11 0.00
## pet* 17 1823 4.96 2.03 4.00 5.02 2.97 1.00
## covid_pos 18 2073 1.85 3.22 0.00 1.12 0.00 0.00
## covid_neg 19 2073 1.11 1.82 0.00 0.75 0.00 0.00
## big5_open 20 1759 5.20 1.14 5.33 5.28 0.99 1.00
## big5_con 21 1756 4.79 1.19 4.67 4.82 1.48 1.00
## big5_agr 22 1761 4.98 1.11 5.00 5.02 0.99 1.00
## big5_neu 23 1749 4.41 1.51 4.67 4.46 1.48 1.00
## big5_ext 24 1756 4.35 1.45 4.33 4.39 1.48 1.00
## pswq 25 1676 2.77 0.80 2.79 2.77 0.95 1.00
## iou 26 1510 2.59 0.92 2.44 2.53 0.99 1.00
## mfq_26 27 1442 4.28 0.70 4.35 4.31 0.67 1.00
## mfq_state 28 1291 4.12 0.98 4.25 4.17 0.93 1.00
## rse 29 1441 2.61 0.72 2.70 2.63 0.74 1.00
## school_covid_support 30 169 1.53 0.34 1.60 1.54 0.30 1.00
## school_att 31 166 3.93 1.07 4.00 3.91 1.21 1.40
## pas_covid 32 1348 3.23 0.69 3.22 3.25 0.66 1.00
## pss 33 1356 2.94 0.94 3.00 2.94 1.11 1.00
## phq 34 1312 2.07 0.86 1.89 1.99 0.99 1.00
## gad 35 1327 2.04 0.91 1.71 1.95 0.85 1.00
## edeq12 36 1251 1.89 0.73 1.75 1.82 0.74 1.00
## brs 37 436 2.68 0.88 2.67 2.68 0.99 1.00
## swemws 38 1235 3.16 0.84 3.14 3.18 0.85 1.00
## isolation_a 39 382 1.68 0.70 1.50 1.59 0.74 1.00
## isolation_c 40 925 2.34 0.82 2.25 2.35 1.11 1.00
## support 41 1298 3.57 0.95 3.67 3.62 0.99 1.00
## fake 42 1275 3.32 3.41 1.94 2.67 1.78 0.25
## pss_std 43 1356 0.00 1.00 0.06 -0.01 1.18 -2.06
## phq_std 44 1312 0.00 1.00 -0.21 -0.09 1.14 -1.24
## support_std 45 1298 0.00 1.00 0.11 0.05 1.04 -2.70
## swemws_std 46 1235 0.00 1.00 -0.02 0.03 1.01 -2.56
## max range skew kurtosis se
## X 8867.00 8866.00 -0.33 -1.15 56.81
## gender* 4.00 3.00 1.61 0.98 0.02
## trans* 3.00 2.00 3.63 12.53 0.01
## sexual_orientation* 6.00 5.00 -0.33 0.94 0.02
## ethnicity* 9.00 8.00 -1.70 1.13 0.06
## relationship_status* 5.00 4.00 -0.68 -1.32 0.04
## age* 5.00 4.00 0.96 -0.91 0.04
## urban_rural* 4.00 3.00 -0.60 -0.99 0.03
## income* 4.00 3.00 0.26 -0.84 0.04
## education* 7.00 6.00 0.62 -1.03 0.05
## employment* 6.00 5.00 1.45 1.63 0.03
## treatment* 6.00 5.00 1.26 1.99 0.03
## health* 11.00 10.00 -2.27 5.83 0.03
## mhealth* 8.00 7.00 -1.46 2.64 0.03
## sleep_hours* 5.00 4.00 0.01 -0.46 0.03
## exercise 333333.00 333333.00 33.87 1146.00 289.10
## pet* 8.00 7.00 -0.14 -1.48 0.05
## covid_pos 15.00 15.00 1.71 1.97 0.07
## covid_neg 8.00 8.00 1.42 0.79 0.04
## big5_open 7.00 6.00 -0.70 0.41 0.03
## big5_con 7.00 6.00 -0.23 -0.30 0.03
## big5_agr 7.00 6.00 -0.39 -0.02 0.03
## big5_neu 7.00 6.00 -0.31 -0.73 0.04
## big5_ext 7.00 6.00 -0.24 -0.77 0.03
## pswq 5.00 4.00 -0.01 -0.77 0.02
## iou 5.00 4.00 0.47 -0.59 0.02
## mfq_26 6.00 5.00 -0.53 0.90 0.02
## mfq_state 6.00 5.00 -0.54 0.18 0.03
## rse 4.00 3.00 -0.18 -0.71 0.02
## school_covid_support 2.00 1.00 -0.15 -1.23 0.03
## school_att 6.33 4.93 0.11 -0.64 0.08
## pas_covid 5.00 4.00 -0.18 0.00 0.02
## pss 5.00 4.00 0.07 -0.73 0.03
## phq 4.00 3.00 0.63 -0.69 0.02
## gad 4.00 3.00 0.68 -0.72 0.02
## edeq12 4.00 3.00 0.68 -0.55 0.02
## brs 5.00 4.00 0.08 -0.68 0.04
## swemws 5.00 4.00 -0.28 -0.26 0.02
## isolation_a 3.50 2.50 0.89 -0.30 0.04
## isolation_c 3.50 2.50 -0.10 -1.25 0.03
## support 5.00 4.00 -0.42 -0.56 0.03
## fake 20.00 19.75 1.87 3.64 0.10
## pss_std 2.18 4.24 0.07 -0.73 0.03
## phq_std 2.23 3.47 0.63 -0.69 0.03
## support_std 1.51 4.21 -0.42 -0.56 0.03
## swemws_std 2.19 4.75 -0.28 -0.26 0.03# also use histograms to examine your continuous variables
hist(d$pss)
hist(d$phq)
# last, use scatterplots to examine your continuous variables together
plot(d$pss, d$phq)
11 Run a Simple Regression
# to calculate standardized coefficients, we have to standardize our IV
d$pss_std <- scale(d$pss, center=T, scale=T)
hist(d$pss_std)
# use the lm() command to run the regression
# dependent/outcome variable on the left, idependent/predictor variable on the right
reg_model <- lm(phq ~ pss_std, data = d)12 Regression: Check Your Assumptions
12.1 Simple Regression Assumptions
- Should have two measurements for each participant
- Variables should be continuous and normally distributed
- Outliers should be identified and removed
- Relationship between the variables should be linear
- Residuals should be normal and have constant variance note: we will not be evaluating whether our data meets these assumptions in this lab/homework – we’ll come back to them next week when we talk about multiple linear regression
12.2 Create plots and view residuals
model.diag.metrics <- augment(reg_model)
ggplot(model.diag.metrics, aes(x = pss_std, y = phq)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
geom_segment(aes(xend = pss_std, yend = .fitted), color = "red", size = 0.3)## 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.## `geom_smooth()` using formula = 'y ~ x'
12.3 Check linearity with Residuals vs Fitted plot
This plot (below) shows the residuals for each case and the fitted line. The red line is the average residual for the specified point of the dependent variable. If the assumption of linearity is met, the red line should be horizontal. This indicates that the residuals average to around zero. You can see that for this lab, the plot shows some non-linearity because there are more datapoints below the regression line than here are above it. Thus, there are some negative residuals that don’t have positive residuals to cancel them out. However, a bit of deviation is okay – just like with skewness and kurtosis, there’s a range that we can work in before non-normality or non-linearity becomes a critical issue.
For some examples of good Residuals vs Fitted plot and ones that show serious errors, check out this page. I’ve included the images in our video and talk about them more in-depth there. But to summarize quickly, you can see the first case has a plot in which the red line sticks pretty closely to the zero line, while the other cases show some serious deviation. Ours is much closer to the ‘good’ plot than it is to the ‘serious issues’ plots. So we’ll consider our data okay and proceed with our analysis. Obviously, this is quite subjective. I’ll talk a bit about why this is in the video, but the key takeaway is that these evaluations are closely tied to the context of our sample, our data, and what we’re studying. It’s almost always a judgement call.
You’ll notice in the bottom right corner, there are some points with numbers included: these are cases or participants (indicated by row number) who have the most influence on the regression line (and so they might outliers).We’ll talk more about outliers in the next section.
To summarize: our plot suggests there is linearity, tghe plot is closer to the good plots as their are residual;s on both the negative and positive side, but they balance eachother out.
plot(reg_model, 1)
12.4 Check for outliers
The plots below both address leverage, or how much each data point is able to influence the regression line. Outliers are points that have undue influence on the regression line, the way that Bill Gates entering the room has an undue influence on the mean income.
The first plot, Cook’s distance, is a visualization of a score called (you guessed it) Cook’s distance, calculated for each case (aka row or participant) in the dataframe. Cook’s distance tells us how much the regression would change if the point was removed. Ideally, we want all points to have the same influence on the regression line, although we accept that there will be some variability. The cutoff for a high Cook’s distance score is .5 (not .05, which is our cutoff for statistical significance). For our data, some points do exert more influence than others but they’re generally equal, and none of them are close to the cutoff.
The second plot also includes the residuals in the examination of leverage. The standardized residuals are on the y-axis and leverage is on the x-axis; this shows us which points have high residuals (are far from the regression line) and high leverage. Point that have large residuals and high leverage are especially worrisome, because they are far from the regression line but are also exerting a large influence on it. The red line indicates the average residual across points with the same amount of leverage. As usual, we want this line to stay as close to the mean line (or the zero line) as possible.
Because the leverage in our plot is low, part of it actually cut off! If you check the first set of plots on this page (note that Residuals vs Leverage is the fourth in the grid) you can see there are curved red lines in the corners of the Residuals vs Leverage plots. This is the .5 cutoff for Cook’s distance, and so any points appearing past these lines is a serious outlier that needs to be removed. On this page you can also see Residuals vs Leverage plots with severe deviations from the mean line, which makes our deviations appear much less serious.
Our data doesn’t have any severe outliers. For your homework, you’ll simply need to generate these plots, assess Cook’s distance in your dataset, and then identify any potential cases that are prominent outliers. Since we have some cutoffs, that makes this process is a bit less subjective than some of the other assessments we’ve done here, which is a nice change!
# Cook's distance
plot(reg_model, 4)
# Residuals vs Leverage
plot(reg_model, 5)
12.5 Issues with My Data
Before interpreting our results, we assessed our variables to see if they met the assumptions for a simple linear regression. Analysis of a Residuals vs Fitted plot suggested that there is some minor non-linearity, but not enough to violate the assumption of linearity. We also checked Cook’s distance and a Residuals vs Leverage plot to detect outliers. there were three cases that had a large residuals and above-average leverage but was below the recommended cutoff for Cook’s distance.
13 Regression: View Test Output
summary(reg_model)##
## Call:
## lm(formula = phq ~ pss_std, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.53425 -0.39859 -0.05205 0.34817 2.56284
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.08084 0.01585 131.26 <2e-16 ***
## pss_std 0.64252 0.01577 40.74 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5729 on 1304 degrees of freedom
## (767 observations deleted due to missingness)
## Multiple R-squared: 0.56, Adjusted R-squared: 0.5597
## F-statistic: 1660 on 1 and 1304 DF, p-value: < 2.2e-16# note for section below: to type lowercase Beta below (ß) you need to hold down Alt key and type 225 on numeric keypad. If that doesn't work you should be able to copy/paste it from somewhere else14 Regression: Write Up Results
To test our hypothesis that stress (measured by the PSS) will significantly predict depression symptoms (measured by PHQ), and that the relationship will be positive, we used a simple linear regression to model the relationship between the variables. We confirmed that our data met the assumptions of a linear regression, checking the linearity of the relationship using a Residuals vs Fitted plot and checking for outliers using Cook’s distance and a Residuals vs Leverage plot. Note: we are skipping the assumptions of normality and homogeneity of variance for this assignment.
As predicted, we found that stress significantly predicted depression symptoms, Adj. R2 = .16, F(1,304) = 1660, p < .001. The relationship between stress and depression symptoms was positive, ß = .64, t(3104) = 40.74, p < .001 (refer to Figure 1). According to Cohen (1988), this constitutes a large effect size (> .50).
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.