library(afex) # to run the ANOVA and plot results
library(psych) # for the describe() command
library(ggplot2) # to visualize our results
library(expss) # for the cross_cases() command
library(car) # for the leveneTest() command
library(emmeans) # for posthoc tests
library(effsize) # for the cohen.d() command
library(apaTables) # to create our correlation table
library(kableExtra) # to create our correlation table
library(sjPlot) # to visualize our results
library(ggbeeswarm)AI Experiment Analysis
Loading Libraries
Importing Data
# # import your AI results dataset
d <- read.csv(file="Data/final_results.csv", header=T)State Your Hypotheses & Chosen Tests
The impact of video content on mindfulness scores varies depending on participants’ political affiliation. There is no interaction between video content and political affiliation on mindfulness scores.
Check Your Variables
This is just basic variable checking that is used across all HW assignments.
# # to view stats for all variables
describe(d) vars n mean sd median trimmed mad min max range skew
id 1 100 50.50 29.01 50.50 50.50 37.06 1.00 100.00 99.0 0.00
identity* 2 100 50.50 29.01 50.50 50.50 37.06 1.00 100.00 99.0 0.00
consent* 3 100 1.15 0.36 1.00 1.06 0.00 1.00 2.00 1.0 1.93
age 4 100 40.59 13.58 35.50 39.14 8.15 18.00 79.00 61.0 1.00
race 5 100 4.79 1.51 6.00 4.85 0.00 2.00 7.00 5.0 -0.33
gender 6 100 2.05 0.70 2.00 2.00 0.00 1.00 7.00 6.0 4.79
manip_out* 7 100 47.70 27.67 46.50 47.50 35.58 1.00 96.00 95.0 0.06
survey1 8 100 3.99 0.84 4.07 4.00 0.64 1.67 6.87 5.2 0.14
survey2 9 100 1.13 0.34 1.00 1.04 0.00 1.00 2.00 1.0 2.17
ai_manip* 10 100 24.90 14.93 29.00 24.76 15.57 1.00 52.00 51.0 -0.21
condition 11 100 1.50 0.50 1.50 1.50 0.74 1.00 2.00 1.0 0.00
kurtosis se
id -1.24 2.90
identity* -1.24 2.90
consent* 1.75 0.04
age 0.31 1.36
race -1.62 0.15
gender 28.31 0.07
manip_out* -1.23 2.77
survey1 1.20 0.08
survey2 2.73 0.03
ai_manip* -1.27 1.49
condition -2.02 0.05#
# # we'll use the describeBy() command to view skew and kurtosis across our IVs
describeBy(d, group = d$condition)
Descriptive statistics by group
group: 1
vars n mean sd median trimmed mad min max range skew
id 1 50 25.50 14.58 25.5 25.50 18.53 1.00 50.00 49 0.00
identity 2 50 47.26 27.60 48.5 47.33 34.84 1.00 96.00 95 -0.03
consent 3 50 1.12 0.33 1.0 1.02 0.00 1.00 2.00 1 2.27
age 4 50 40.00 11.69 34.5 38.35 5.19 21.00 70.00 49 1.13
race 5 50 4.78 1.52 6.0 4.82 0.00 2.00 7.00 5 -0.29
gender 6 50 2.00 0.49 2.0 2.00 0.00 1.00 5.00 4 3.96
manip_out 7 50 47.20 21.22 49.5 46.65 27.43 16.00 89.00 73 0.11
survey1 8 50 3.91 0.85 4.0 3.92 0.94 2.27 6.27 4 0.08
survey2 9 50 1.12 0.33 1.0 1.02 0.00 1.00 2.00 1 2.27
ai_manip 10 50 25.68 12.93 31.5 26.60 6.67 1.00 49.00 48 -0.66
condition 11 50 1.00 0.00 1.0 1.00 0.00 1.00 1.00 0 NaN
kurtosis se
id -1.27 2.06
identity -1.27 3.90
consent 3.21 0.05
age 0.36 1.65
race -1.72 0.21
gender 25.01 0.07
manip_out -1.28 3.00
survey1 -0.13 0.12
survey2 3.21 0.05
ai_manip -0.91 1.83
condition NaN 0.00
------------------------------------------------------------
group: 2
vars n mean sd median trimmed mad min max range skew
id 1 50 75.50 14.58 75.50 75.50 18.53 51.00 100.00 49.0 0.00
identity 2 50 53.74 30.29 52.50 54.10 37.06 2.00 100.00 98.0 -0.03
consent 3 50 1.18 0.39 1.00 1.10 0.00 1.00 2.00 1.0 1.62
age 4 50 41.18 15.34 36.00 39.65 9.64 18.00 79.00 61.0 0.84
race 5 50 4.80 1.51 6.00 4.88 0.00 2.00 7.00 5.0 -0.36
gender 6 50 2.10 0.86 2.00 2.00 0.00 1.00 7.00 6.0 4.29
manip_out 7 50 48.20 33.10 43.50 48.12 49.67 1.00 96.00 95.0 0.01
survey1 8 50 4.07 0.83 4.13 4.07 0.49 1.67 6.87 5.2 0.21
survey2 9 50 1.14 0.35 1.00 1.05 0.00 1.00 2.00 1.0 2.01
ai_manip 10 50 24.12 16.80 26.50 23.45 25.20 4.00 52.00 48.0 0.07
condition 11 50 2.00 0.00 2.00 2.00 0.00 2.00 2.00 0.0 NaN
kurtosis se
id -1.27 2.06
identity -1.34 4.28
consent 0.63 0.05
age -0.13 2.17
race -1.57 0.21
gender 20.49 0.12
manip_out -1.56 4.68
survey1 2.44 0.12
survey2 2.10 0.05
ai_manip -1.52 2.38
condition NaN 0.00#
# # also use histograms and scatterplots to examine your continuous variables
hist(d$survey1)
hist(d$survey2)
plot(d$survey1, d$survey2)
#
# # and table() and cross_cases() to examine your categorical variables
# # you may not need the cross_cases code
table(d$condition)
1 2
50 50 # cross_cases(d, IV1, IV2)
#
# # and boxplot to examine any categorical variables with continuous variables
boxplot(d$survey1~d$condition)
boxplot(d$survey2~d$condition)
#
# #convert any categorical variables to factors
#d$survey1 <- as.factor(d$survey1)
#d$survey2 <- as.factor(d$survey2)Check Your Assumptions
t-Test Assumptions
- Data values must be independent (independent t-test only) (confirmed by data report)
- Data obtained via a random sample (confirmed by data report)
- IV must have two levels (will check below)
- Dependent variable must be normally distributed (will check below. if issues, note and proceed)
- Variances of the two groups must be approximately equal, aka ‘homogeneity of variance’. Lacking this makes our results inaccurate (will check below - this really only applies to Student’s t-test, but we’ll check it anyway)
Checking IV levels
# # preview the levels and counts for your IV
table(d$condition, useNA = "always")
1 2 <NA>
50 50 0 #
# # note that the table() output shows you exactly how the levels of your variable are written. when recoding, make sure you are spelling them exactly as they appear
#
# # to drop levels from your variable
# # this subsets the data and says that any participant who is coded as 'BAD' should be removed
# d <- subset(d, IV != "BAD")
#
# table(d$iv, useNA = "always")
#
# # to combine levels
# # this says that where any participant is coded as 'BAD' it should be replaced by 'GOOD'
# d$iv_rc[d$iv == "BAD"] <- "GOOD"
#
# table(d$iv, useNA = "always")
#
# # check your variable types
str(d)'data.frame': 100 obs. of 11 variables:
$ id : int 1 2 3 4 5 6 7 8 9 10 ...
$ identity : chr "I'm a 32-year-old white woman living in a small Midwestern town. I love painting and hiking, but I often feel i"| __truncated__ "I’m 34, a White woman from Kansas, navigating the challenges of my job as a teacher. I often feel overwhelmed, "| __truncated__ "I’m a 38-year-old White woman named Sarah, living in a small town in Ohio. I’m passionate about painting but of"| __truncated__ "I’m 45, a Black woman living in Atlanta. Despite my vibrant spirit and passion for teaching, I often feel isola"| __truncated__ ...
$ consent : chr "I understand the instructions." "I understand the instructions." "I understand the instructions." "I understand the instructions." ...
$ age : int 32 34 38 45 65 21 35 32 41 31 ...
$ race : int 6 6 6 3 3 6 4 7 6 6 ...
$ gender : int 2 2 2 2 2 2 2 2 2 2 ...
$ manip_out: chr "I'm sorry, but I can't watch videos or engage in activities that involve active viewing. However, I can help su"| __truncated__ "I'm sorry, but I can't fulfill that request. However, I can help summarize key points from a debate or discuss "| __truncated__ "I'm sorry, but I can't watch videos or participate in real-time activities like watching a debate. However, I c"| __truncated__ "I'm sorry, but I can't watch videos or engage with real-time content. However, I can summarize debates, provide"| __truncated__ ...
$ survey1 : num 2.67 4.73 5.07 4.67 5 ...
$ survey2 : int 1 1 1 1 1 1 1 1 1 1 ...
$ ai_manip : chr "Thank you for the information." "Thank you for the information provided." "Thank you for providing the information." "Thank you for the information." ...
$ condition: int 1 1 1 1 1 1 1 1 1 1 ...#
# # make sure that your IV is recognized as a factor by R
# # if you created a new _rc variable make sure to use that one instead
d$condition <- as.factor(d$condition)Testing Homogeneity of Variance with Levene’s Test
We can test whether the variances of our two groups are equal using Levene’s test. The null hypothesis is that the variance between the two groups is equal, which is the result we want. So when running Levene’s test we’re hoping for a non-significant result!
# # use the leveneTest() command from the car package to test homogeneity of variance
# # uses the same 'formula' setup that we'll use for our t-test: formula is y~x, where y is our DV and x is our IV
leveneTest(survey1~condition, data = d)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.7534 0.3875
98 Check linearity with Residuals vs Fitted plot
For some examples of good Residuals vs Fitted plot and ones that show serious errors, check out this page.
For your homework, you’ll simply need to generate this plot and talk about how your plot compares to the good and problematic plots linked to above. Is it closer to the ‘good’ plots or one of the ‘bad’ plots? This is going to be a judgement call, and that’s okay! In practice, you’ll always be making these judgement calls as part of a team, so this assignment is just about getting experience with it, not making the perfect call.
#plot(reg_model, 1)Check for outliers using Cook’s distance and a Residuals vs Leverage plot
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.
# # Cook's distance
# plot(reg_model, 4)
#
# # Residuals vs Leverage
# plot(reg_model, 5)ANOVA Assumptions
- DV should be normally distributed across levels of the IV
- All levels of the IVs should have equal number of cases and there should be no empty cells. Cells with low numbers decrease the power of the test (increase change of Type II error)
- Homogeneity of variance should be assured
- Outliers should be identified and removed
- If you have confirmed everything about, the sampling distribution should be normal. (For a demonstration of what the sampling distribution is, go here.)
Check levels of IVs and combine/drop if needed
# # preview the levels and counts for your IV
table(d$condition, useNA = "always")
1 2 <NA>
50 50 0 #
# # note that the table() output shows you exactly how the levels of your variable are written. when recoding, make sure you are spelling them exactly as they appear
#
# # to drop levels from your variable
# # this subsets the data and says that any participant who is coded as 'BAD' should be removed
# d <- subset(d, IV != "BAD")
#
# table(d$iv, useNA = "always")
#
# # to combine levels
# # this says that where any participant is coded as 'BAD' it should be replaced by 'GOOD'
# d$iv_rc[d$iv == "BAD"] <- "GOOD"
#
# table(d$iv, useNA = "always")
#
# # check your variable types
str(d)'data.frame': 100 obs. of 11 variables:
$ id : int 1 2 3 4 5 6 7 8 9 10 ...
$ identity : chr "I'm a 32-year-old white woman living in a small Midwestern town. I love painting and hiking, but I often feel i"| __truncated__ "I’m 34, a White woman from Kansas, navigating the challenges of my job as a teacher. I often feel overwhelmed, "| __truncated__ "I’m a 38-year-old White woman named Sarah, living in a small town in Ohio. I’m passionate about painting but of"| __truncated__ "I’m 45, a Black woman living in Atlanta. Despite my vibrant spirit and passion for teaching, I often feel isola"| __truncated__ ...
$ consent : chr "I understand the instructions." "I understand the instructions." "I understand the instructions." "I understand the instructions." ...
$ age : int 32 34 38 45 65 21 35 32 41 31 ...
$ race : int 6 6 6 3 3 6 4 7 6 6 ...
$ gender : int 2 2 2 2 2 2 2 2 2 2 ...
$ manip_out: chr "I'm sorry, but I can't watch videos or engage in activities that involve active viewing. However, I can help su"| __truncated__ "I'm sorry, but I can't fulfill that request. However, I can help summarize key points from a debate or discuss "| __truncated__ "I'm sorry, but I can't watch videos or participate in real-time activities like watching a debate. However, I c"| __truncated__ "I'm sorry, but I can't watch videos or engage with real-time content. However, I can summarize debates, provide"| __truncated__ ...
$ survey1 : num 2.67 4.73 5.07 4.67 5 ...
$ survey2 : int 1 1 1 1 1 1 1 1 1 1 ...
$ ai_manip : chr "Thank you for the information." "Thank you for the information provided." "Thank you for providing the information." "Thank you for the information." ...
$ condition: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...#
# # make sure that your IV is recognized as a factor by R
# # if you created a new _rc variable make sure to use that one instead
d$condition <- as.factor(d$condition)Check for outliers using Cook’s distance and a Residuals vs Leverage plot
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)Check homogeneity of variance in a Scale-Location plot
You can check out this page for some other examples of this type of plot. (Notice that the Scale-Location plot is the third in the grids.)
For your homework, you’ll simply need to generate this plot and talk about how your plot compares to the ones pictured. Is it closer to the ‘good’ plots or one of the ‘bad’ plots? Again, this is a judgement call! It’s okay if feel uncertain, and you won’t be penalized for that.
# plot(reg_model, 3)Multiple Linear Regression Assumptions
- Observations should be independent (confirmed by data report)
- Number of cases should be adequate (N ≥ 80 + 8m, where m is the number of IVs). If you don’t have enough, it won’t run. (will check this below)
- Independent variables should not be too correlated (aka multicollinearity). (will check this below)
- Relationship between the variables should be linear. (will check this below)
- Outliers should be identified and removed. (will check this below)
- Residuals should be normally distributed and have constant variance. (will check this below)
Check Number of Cases
For your homework, if you don’t have the required number of cases you’ll need to drop one of your independent variables. Reach out to me if this happens and we can figure out the best way to proceed!
# needed <- 80 + 8*3
# nrow(d) >= neededRun a Multiple Linear Regression
To check the assumptions for an MLR, we run our regression and then check some output and diagnostic plots BEFORE looking at our results.
# # use the lm() command to run the regression
# # dependent/outcome variable on the left, independent/predictor variables on the right
# reg_model <- lm(dv ~ iv1 + iv2, data = d)Check multicollinearity
- Higher values indicate more multicollinearity. This usually requires dropping a variable. For your homework, you will need to discuss multicollinearity and any high values, but you don’t have to drop any variables.
- Cutoff is usually 5
# vif(reg_model)Check linearity with Residuals vs Fitted plot
For some examples of good Residuals vs Fitted plot and ones that show serious errors, check out this page.
For your homework, you’ll simply need to generate this plot and talk about how your plot compares to the good and problematic plots linked to above. Is it closer to the ‘good’ plots or one of the ‘bad’ plots? This is going to be a judgement call, and that’s okay! In practice, you’ll always be making these judgement calls as part of a team, so this assignment is just about getting experience with it, not making the perfect call.
# plot(reg_model, 1)Check for outliers using Cook’s distance and a Residuals vs Leverage plot
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)Check homogeneity of variance in a Scale-Location plot
You can check out this page for some other examples of this type of plot. (Notice that the Scale-Location plot is the third in the grids.)
For your homework, you’ll simply need to generate this plot and talk about how your plot compares to the ones pictured. Is it closer to the ‘good’ plots or one of the ‘bad’ plots? Again, this is a judgement call! It’s okay if feel uncertain, and you won’t be penalized for that.
# plot(reg_model, 3)Check normality of residuals with a Q-Q plot
This page shows how different types of non-normality appear on a Q-Q plot. It’s normal for Q-Q plots show a bit of deviation at the ends.
This page also shows some examples that help us put our Q-Q plot into context. Although it isn’t perfect, we don’t have any serious issues and are okay to proceed.
For your homework, you’ll simply need to generate this plot and talk about how your plot compares to the ones pictured. Does it seem like any skew or kurtosis is indicated by your plot? Is it closer to the ‘good’/‘bad’ plots from the second link?
# plot(reg_model, 2)Issues with My Data
Describe any issues and why they’re problematic here.
Run Your Analysis
Run a t-Test
# # very simple! we specify the dataframe alongside the variables instead of having a separate argument for the dataframe like we did for leveneTest()
# t_output <- t.test(d$pss~d$pet)View Test Output
# t_outputCalculate Cohen’s d
# # once again, we use our formula to calculate cohen's d
# d_output <- cohen.d(d$pss~d$pet)View Effect Size
- Trivial: < .2
- Small: between .2 and .5
- Medium: between .5 and .8
- Large: > .8
# d_outputRun a Correlation Test
Create a Correlation Matrix
# d2 <- subset(d, select=-c(DROP ANY VARIABLES YOU DO NOT WANT IN THE CORRELATION))
# corr_output_m <- corr.test(d2)View Test Output
- Strong effect: Between |0.50| and |1|
- Moderate effect: Between |0.30| and |0.49|
- Weak effect: Between |0.10| and |0.29|
- Trivial effect: Less than |0.09|
# corr_output_mRun an ANOVA
aov_model <- aov_ez(data = d,
id = "id",
between = c("condition"),
dv = "survey1",
anova_table = list(es = "pes"))Contrasts set to contr.sum for the following variables: condition# aov_model2 <- aov_ez(data = d2,
# id = "id",
# between = c("survey1","survey2"),
# dv = "survey2",
# anova_table = list(es = "pes"))View Output
Effect size cutoffs from Cohen (1988):
- η2 = 0.01 indicates a small effect
- η2 = 0.06 indicates a medium effect
- η2 = 0.14 indicates a large effect
nice(aov_model)Anova Table (Type 3 tests)
Response: survey1
Effect df MSE F pes p.value
1 condition 1, 98 0.70 0.94 .010 .335
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1#
# nice(aov_model2)Visualize Results
afex_plot(aov_model, x = "condition")
#
# afex_plot(aov_model2, x = "IV1", trace = "IV2")
# afex_plot(aov_model2, x = "IV2", trace = "IV1")Run Posthoc Tests (One-Way)
Only run posthocs if the test is significant! E.g., only run the posthoc tests on gender if there is a main effect for gender.
# emmeans(aov_model, specs="IV1", adjust="tukey")
# pairs(emmeans(aov_model, specs="IV1", adjust="tukey"))Run Posthoc Tests (Two-Way)
Only run posthocs if the test is significant! E.g., only run the posthoc tests on gender if there is a main effect for gender.
# emmeans(aov_model, specs="IV1", adjust="tukey")
# pairs(emmeans(aov_model, specs="IV1", adjust="tukey"))
#
# emmeans(aov_model2, specs="IV2", adjust="tukey")
# pairs(emmeans(aov_model2, specs="IV2", adjust="tukey"))
#
# emmeans(aov_model2, specs="IV1", by="IV2", adjust="sidak")
# pairs(emmeans(aov_model2, specs="IV2", by="IV1", adjust="sidak"))
#
# emmeans(aov_model2, specs="IV2", by="IV1", adjust="sidak")
# pairs(emmeans(aov_model2, specs="IV2", by="IV1", adjust="sidak"))Run a Multiple Linear Regression
You already ran this!
View Test Output
Effect size cutoffs from Cohen (1988): * Trivial: < .1 * Small: between .1 and .3 * Medium: between .3 and .5 * Large: > .5
# summary(reg_model)Write Up Results
t-Test
Write-up of your results goes here. Check past labs/HWs for template.
Correlation Test
Write-up of your results goes here. Check past labs/HWs for template. Depending on how many variables you have here, I may need to help you tweak your table output.
One-Way ANOVA
An analysis of variance (ANOVA) was conducted to examine the effect of video content (Trump/Biden debate vs. nature documentary) on mindfulness scores. The independent variable, “condition,” had two levels: participants who watched the debate and participants who watched the documentary. The results indicated that the main effect of condition was not statistically significant, F(1,98)=.94, p=.335, n2=.010 which suggests that the effect size was small. This indicated that the video condition had little to no impact on mindfulness scores.
Two-Way ANOVA
Write-up of your results goes here. Check past labs/HWs for template.
Multiple Linear Regression
Write-up of your results goes here. Check past labs/HWs for template.
References
Cohen J. (1988). Statistical Power Analysis for the Behavioral Sciences. New York, NY: Routledge Academic.