library(psych) # for the describe() command
library(broom) # for the augment() command
library(ggplot2) # to visualize our results
##
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
##
## %+%, alpha
# import the dataset you cleaned previously
# this will be the dataset you'll use throughout the rest of the semester
# use EAMMi2 data
d <- read.csv(file="Data/mydata.csv", header=T)
We hypothesize that mental flexibility (measured by the mfq_26) will significantly predict self-esteem (measured by the rse), and that the relationship will be positive.
# 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': 1337 obs. of 6 variables:
## $ sexual_orientation : chr "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" "Heterosexual/Straight" ...
## $ relationship_status: chr "In a relationship/married and cohabiting" "Prefer not to say" "Prefer not to say" "In a relationship/married and cohabiting" ...
## $ big5_open : num 5.33 5.33 5 6 5 ...
## $ pswq : num 4.94 3.36 1.86 3.94 2.62 ...
## $ mfq_26 : num 4.2 3.35 4.65 4.65 4.5 4.3 5.25 5 4.7 4.05 ...
## $ rse : num 2.3 1.6 3.9 1.7 3.9 2.4 1.8 1.3 3.5 2.6 ...
# 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 max range
## sexual_orientation* 1 1337 3.78 1.05 4.00 3.79 0.00 1 6 5
## relationship_status* 2 1337 3.73 1.68 5.00 3.91 0.00 1 5 4
## big5_open 3 1337 5.23 1.12 5.33 5.31 0.99 1 7 6
## pswq 4 1337 2.75 0.79 2.79 2.75 0.95 1 5 4
## mfq_26 5 1337 4.28 0.70 4.35 4.31 0.67 1 6 5
## rse 6 1337 2.61 0.72 2.70 2.62 0.74 1 4 3
## skew kurtosis se
## sexual_orientation* -0.44 0.99 0.03
## relationship_status* -0.74 -1.24 0.05
## big5_open -0.69 0.33 0.03
## pswq 0.01 -0.77 0.02
## mfq_26 -0.55 0.96 0.02
## rse -0.17 -0.72 0.02
# also use histograms to examine your continuous variables
hist(d$mfq_26)
hist(d$rse)
# last, use scatterplots to examine your continuous variables together
plot(d$mfq_26, d$rse)
# to calculate standardized coefficients, we have to standardize our IV
d$mfq_26 <- scale(d$mfq_26, center=T, scale=T)
d$rse <- scale(d$rse, center=T, scale=T)
# use the lm() command to run the regression
# dependent/outcome variable on the left, idependent/predictor variable on the right
reg_model <- lm(rse ~ mfq_26, data = d)
model.diag.metrics <- augment(reg_model)
ggplot(model.diag.metrics, aes(x = mfq_26, y = rse)) +
geom_point() +
stat_smooth(method = lm, se = FALSE) +
geom_segment(aes(xend = mfq_26, 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'
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 homework, the plot shows minor 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 top left 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 a little minor non-linearity. However, it is still closer to the ‘good’ plots.
plot(reg_model, 1)
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 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 high, 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.
Even though our data may have some sever outliers that pull on the red line which makes it high, none of the points are close to the cutoff for Cook’s distance.
# Cook's distance
plot(reg_model, 4)
# Residuals vs Leverage
plot(reg_model, 5)
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 little 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. Two cases had a large residuals and above-average leverage but were below the recommended cutoff for Cook’s distance.
Trivial: Less than 0.10 Small: 0.10–0.29 Medium: 0.30–0.49 Large: 0.50 or greater
summary(reg_model)
##
## Call:
## lm(formula = rse ~ mfq_26, data = d)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5508 -0.5732 0.0485 0.5678 3.3716
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.404e-16 2.216e-02 0.00 1
## mfq_26 5.864e-01 2.217e-02 26.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8103 on 1335 degrees of freedom
## Multiple R-squared: 0.3439, Adjusted R-squared: 0.3434
## F-statistic: 699.7 on 1 and 1335 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 else
To test our hypothesis that mental flexibility (measured by the mfq_26) will significantly predict self-esteem (measured by the rse), 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 mental flexibility significantly predicted self-esteem, Adj. R2 = .34, F(1,1335) = 699.7, p < .001. The relationship between mental flexibility and self-esteem was positive, ß = .59, t(1335) = 26.45, 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.