PASA Analysis for Equanimity Study
Megan Lipsett
Fall 2020
Analysis of PASA Scale for Equanimity participants
#Analysis of the Primary and Secondary Appraisal Scale (continuous) for the Equanimity Study (categorical: condition 1 and condition 2)
df_eq_pasa <- read.csv(file = "/Users/mlipsett/Desktop/Equanimity/Database/Paper Measures/Eq 2.7.20 JS.csv", stringsAsFactors = FALSE, header = TRUE) #Reading in the csv file with PASA (and other) data for the equanimity study
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondf_eq_pasa <- dplyr::select(df_eq_pasa, PID, Condition, PASA_no_threat, PASA_task_import, PASA_I_can, PASA_depends_experts, PASA_unpleasant, PASA_do_not_care, PASA_no_idea, PASA_protect, PASA_no_worry, PASA_not_challenge, PASA_solutions, PASA_determine, PASA_scared, PASA_can_solve, PASA_my_effort) #Subsetting the data to include only PASA data #Imputing the mean for PID 355's missing item "PASA_my_effort"
my_effort_mean <- mean(df_eq_pasa$PASA_my_effort, na.rm = TRUE) #Finding the mean across groups for the missing item.
df_eq_pasa$PASA_my_effort[30] = my_effort_mean #Replacing missing item with the mean for PID 355df_eq_pasa <- df_eq_pasa[-c(12, 14, 17, 56, 96),] #Removing participants for whom we are missing the entire PASA measure
#avg_age <- mean(df_eq_pasa$DEMO_age, na.rm = TRUE) #Find mean age of analyzed participants
#avg_age#Reverse Coding: PASA_no_threat, PASA_do_not_care, PASA_no_idea, PASA_no_worry, PASA_not_challenge (1 becomes 6, 2 becomes 5, 3 becomes 4 etc.) the following items:1, 6, 7, 9, 10
df_eq_pasa$PASA_no_threat_r <- 7-df_eq_pasa$PASA_no_threat
df_eq_pasa$PASA_do_not_care_r <- 7-df_eq_pasa$PASA_do_not_care
df_eq_pasa$PASA_no_idea_r <- 7-df_eq_pasa$PASA_no_idea
df_eq_pasa$PASA_no_worry_r <- 7-df_eq_pasa$PASA_no_worry
df_eq_pasa$PASA_not_challenge_r <- 7-df_eq_pasa$PASA_not_challenge
df_eq_pasa <- dplyr::select(df_eq_pasa, PID, Condition, PASA_no_threat_r, PASA_task_import, PASA_I_can, PASA_depends_experts, PASA_unpleasant, PASA_do_not_care_r, PASA_no_idea_r, PASA_protect, PASA_no_worry_r, PASA_not_challenge_r, PASA_solutions, PASA_determine, PASA_scared, PASA_can_solve, PASA_my_effort) #Subsetting the data to remove original items that were reverse codedCorrelation Matrix - All items
# All Items
sapply(df_eq_pasa, is.factor)## PID Condition PASA_no_threat_r
## FALSE FALSE FALSE
## PASA_task_import PASA_I_can PASA_depends_experts
## FALSE FALSE FALSE
## PASA_unpleasant PASA_do_not_care_r PASA_no_idea_r
## FALSE FALSE FALSE
## PASA_protect PASA_no_worry_r PASA_not_challenge_r
## FALSE FALSE FALSE
## PASA_solutions PASA_determine PASA_scared
## FALSE FALSE FALSE
## PASA_can_solve PASA_my_effort
## FALSE FALSEcor <- cor(df_eq_pasa[sapply(df_eq_pasa, function(x) !is.factor(x))], use = "complete.obs") #Correlation matrix of indices
cor## PID Condition PASA_no_threat_r PASA_task_import
## PID 1.000000000 -0.03244415 0.006026812 -0.04619844
## Condition -0.032444151 1.00000000 -0.064448450 -0.02218839
## PASA_no_threat_r 0.006026812 -0.06444845 1.000000000 0.13397340
## PASA_task_import -0.046198438 -0.02218839 0.133973403 1.00000000
## PASA_I_can -0.120599659 -0.24000561 -0.373161964 -0.04787633
## PASA_depends_experts -0.033310761 -0.01312114 0.151269957 0.09433738
## PASA_unpleasant 0.026008553 0.13408774 0.460912052 0.15800210
## PASA_do_not_care_r -0.126475496 -0.10657155 0.153278435 0.55290912
## PASA_no_idea_r -0.063230612 -0.27552453 -0.295166142 0.03696563
## PASA_protect -0.034519304 0.02217707 -0.089099327 -0.17675303
## PASA_no_worry_r 0.006541330 -0.01528477 0.636650562 0.29475544
## PASA_not_challenge_r 0.034410071 -0.03145482 0.525610466 0.15559796
## PASA_solutions 0.067195184 -0.11339297 -0.209538914 0.04463128
## PASA_determine 0.116501174 -0.16086444 -0.065460841 0.01991743
## PASA_scared 0.006276489 -0.00912133 0.565005103 0.23679709
## PASA_can_solve -0.007790543 -0.08891637 -0.225227677 -0.07546448
## PASA_my_effort -0.004289245 -0.10359774 0.113606854 0.21022604
## PASA_I_can PASA_depends_experts PASA_unpleasant
## PID -0.12059966 -0.03331076 0.02600855
## Condition -0.24000561 -0.01312114 0.13408774
## PASA_no_threat_r -0.37316196 0.15126996 0.46091205
## PASA_task_import -0.04787633 0.09433738 0.15800210
## PASA_I_can 1.00000000 -0.09173506 -0.41094829
## PASA_depends_experts -0.09173506 1.00000000 0.27115518
## PASA_unpleasant -0.41094829 0.27115518 1.00000000
## PASA_do_not_care_r 0.03625926 0.01297199 -0.04878697
## PASA_no_idea_r 0.49336921 -0.18974362 -0.27764857
## PASA_protect -0.01220816 0.02913941 -0.16192185
## PASA_no_worry_r -0.35967026 0.08655227 0.51299073
## PASA_not_challenge_r -0.40669496 0.22316793 0.54626874
## PASA_solutions 0.46206318 -0.36614810 -0.38981525
## PASA_determine -0.02749809 -0.07092212 -0.07326456
## PASA_scared -0.37860101 0.27846598 0.73384689
## PASA_can_solve 0.55957293 -0.13523626 -0.38196845
## PASA_my_effort 0.06132531 0.20578019 0.06476351
## PASA_do_not_care_r PASA_no_idea_r PASA_protect
## PID -0.12647550 -0.06323061 -0.03451930
## Condition -0.10657155 -0.27552453 0.02217707
## PASA_no_threat_r 0.15327843 -0.29516614 -0.08909933
## PASA_task_import 0.55290912 0.03696563 -0.17675303
## PASA_I_can 0.03625926 0.49336921 -0.01220816
## PASA_depends_experts 0.01297199 -0.18974362 0.02913941
## PASA_unpleasant -0.04878697 -0.27764857 -0.16192185
## PASA_do_not_care_r 1.00000000 0.04028965 0.14501118
## PASA_no_idea_r 0.04028965 1.00000000 -0.03859003
## PASA_protect 0.14501118 -0.03859003 1.00000000
## PASA_no_worry_r 0.15540535 -0.23143302 -0.28012685
## PASA_not_challenge_r 0.17233069 -0.31923177 -0.08379568
## PASA_solutions 0.08798373 0.35606826 0.13100219
## PASA_determine 0.06561301 0.14497227 0.22397681
## PASA_scared 0.12153035 -0.34350795 -0.05634901
## PASA_can_solve 0.02816083 0.39175166 0.03502423
## PASA_my_effort 0.10742178 -0.06590094 0.03015462
## PASA_no_worry_r PASA_not_challenge_r PASA_solutions
## PID 0.00654133 0.03441007 0.06719518
## Condition -0.01528477 -0.03145482 -0.11339297
## PASA_no_threat_r 0.63665056 0.52561047 -0.20953891
## PASA_task_import 0.29475544 0.15559796 0.04463128
## PASA_I_can -0.35967026 -0.40669496 0.46206318
## PASA_depends_experts 0.08655227 0.22316793 -0.36614810
## PASA_unpleasant 0.51299073 0.54626874 -0.38981525
## PASA_do_not_care_r 0.15540535 0.17233069 0.08798373
## PASA_no_idea_r -0.23143302 -0.31923177 0.35606826
## PASA_protect -0.28012685 -0.08379568 0.13100219
## PASA_no_worry_r 1.00000000 0.64466615 -0.29096932
## PASA_not_challenge_r 0.64466615 1.00000000 -0.32193557
## PASA_solutions -0.29096932 -0.32193557 1.00000000
## PASA_determine -0.04922000 -0.11547278 0.12200164
## PASA_scared 0.56266621 0.61452758 -0.31906180
## PASA_can_solve -0.34515638 -0.41359191 0.40890577
## PASA_my_effort 0.10215267 0.01436219 0.18212188
## PASA_determine PASA_scared PASA_can_solve PASA_my_effort
## PID 0.116501174 0.006276489 -0.007790543 -0.004289245
## Condition -0.160864439 -0.009121330 -0.088916367 -0.103597742
## PASA_no_threat_r -0.065460841 0.565005103 -0.225227677 0.113606854
## PASA_task_import 0.019917430 0.236797092 -0.075464480 0.210226037
## PASA_I_can -0.027498094 -0.378601006 0.559572931 0.061325313
## PASA_depends_experts -0.070922120 0.278465975 -0.135236264 0.205780187
## PASA_unpleasant -0.073264558 0.733846890 -0.381968452 0.064763507
## PASA_do_not_care_r 0.065613009 0.121530347 0.028160831 0.107421780
## PASA_no_idea_r 0.144972266 -0.343507949 0.391751659 -0.065900939
## PASA_protect 0.223976807 -0.056349015 0.035024230 0.030154615
## PASA_no_worry_r -0.049219999 0.562666215 -0.345156382 0.102152668
## PASA_not_challenge_r -0.115472777 0.614527577 -0.413591905 0.014362190
## PASA_solutions 0.122001639 -0.319061803 0.408905767 0.182121880
## PASA_determine 1.000000000 -0.003613629 -0.105924948 0.128253763
## PASA_scared -0.003613629 1.000000000 -0.428473309 0.130147596
## PASA_can_solve -0.105924948 -0.428473309 1.000000000 -0.010316919
## PASA_my_effort 0.128253763 0.130147596 -0.010316919 1.000000000Correlation Matrix Heat Maps
Correlation Matrix Heat Maps
#All items
df_eq_pasa_cor <- dplyr::select(df_eq_pasa, PASA_no_threat_r:PASA_my_effort)
sapply(df_eq_pasa_cor, is.factor)## PASA_no_threat_r PASA_task_import PASA_I_can
## FALSE FALSE FALSE
## PASA_depends_experts PASA_unpleasant PASA_do_not_care_r
## FALSE FALSE FALSE
## PASA_no_idea_r PASA_protect PASA_no_worry_r
## FALSE FALSE FALSE
## PASA_not_challenge_r PASA_solutions PASA_determine
## FALSE FALSE FALSE
## PASA_scared PASA_can_solve PASA_my_effort
## FALSE FALSE FALSEcor <- round(cor(df_eq_pasa_cor [sapply(df_eq_pasa_cor, function(x) !is.factor(x))], use = "complete.obs"),2) #Correlation matrix of indices
#cor
library(reshape2)
melted_cormat <- melt(cor)
#head(melted_cormat)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
cormat[upper.tri(cormat)] <- NA
return(cormat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
cormat[lower.tri(cormat)]<- NA
return(cormat)
}
upper_tri <- get_upper_tri(cor)
#upper_tri
# Melt the correlation matrix
library(reshape2)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cor){
# Use correlation between variables as distance
dd <- as.dist((1-cor)/2)
hc <- hclust(dd)
cor <-cor[hc$order, hc$order]
}
# Reorder the correlation matrix
cor <- reorder_cormat(cor)
upper_tri <- get_upper_tri(cor)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
library(ggplot2)
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
cor_heatmap <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(1, 0),
legend.position = c(0.6, 0.7),
legend.direction = "horizontal")+
guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
title.position = "top", title.hjust = 0.5))
cor_heatmap
# Threat items
df_eq_pasa_threat_cor <- dplyr::select(df_eq_pasa, PASA_no_threat_r, PASA_unpleasant, PASA_no_worry_r, PASA_scared, PASA_not_challenge_r)
sapply(df_eq_pasa_threat_cor, is.factor)## PASA_no_threat_r PASA_unpleasant PASA_no_worry_r
## FALSE FALSE FALSE
## PASA_scared PASA_not_challenge_r
## FALSE FALSEcor <- round(cor(df_eq_pasa_threat_cor [sapply(df_eq_pasa_threat_cor, function(x) !is.factor(x))], use = "complete.obs"),2) #Correlation matrix of indices
#cor
library(reshape2)
melted_cormat <- melt(cor)
#head(melted_cormat)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
cormat[upper.tri(cormat)] <- NA
return(cormat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
cormat[lower.tri(cormat)]<- NA
return(cormat)
}
upper_tri <- get_upper_tri(cor)
#upper_tri
# Melt the correlation matrix
library(reshape2)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cor){
# Use correlation between variables as distance
dd <- as.dist((1-cor)/2)
hc <- hclust(dd)
cor <-cor[hc$order, hc$order]
}
# Reorder the correlation matrix
cor <- reorder_cormat(cor)
upper_tri <- get_upper_tri(cor)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
library(ggplot2)
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
cor_heatmap <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(1, 0),
legend.position = c(0.6, 0.7),
legend.direction = "horizontal")+
guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
title.position = "top", title.hjust = 0.5))
cor_heatmap
# Challenge items
df_eq_pasa_challenge_cor <- dplyr::select(df_eq_pasa, PASA_task_import, PASA_do_not_care_r, PASA_not_challenge_r)
# All Items
sapply(df_eq_pasa_challenge_cor, is.factor)## PASA_task_import PASA_do_not_care_r PASA_not_challenge_r
## FALSE FALSE FALSEcor <- round(cor(df_eq_pasa_challenge_cor[sapply(df_eq_pasa_challenge_cor, function(x) !is.factor(x))], use = "complete.obs"),2) #Correlation matrix of indices
#cor
library(reshape2)
melted_cormat <- melt(cor)
#head(melted_cormat)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
cormat[upper.tri(cormat)] <- NA
return(cormat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
cormat[lower.tri(cormat)]<- NA
return(cormat)
}
upper_tri <- get_upper_tri(cor)
#upper_tri
# Melt the correlation matrix
library(reshape2)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cor){
# Use correlation between variables as distance
dd <- as.dist((1-cor)/2)
hc <- hclust(dd)
cor <-cor[hc$order, hc$order]
}
# Reorder the correlation matrix
cor <- reorder_cormat(cor)
upper_tri <- get_upper_tri(cor)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
library(ggplot2)
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
cor_heatmap <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(1, 0),
legend.position = c(0.6, 0.7),
legend.direction = "horizontal")+
guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
title.position = "top", title.hjust = 0.5))
cor_heatmap
# Controllability items
df_eq_pasa_control_cor <- dplyr::select(df_eq_pasa, PASA_depends_experts, PASA_protect, PASA_determine, PASA_my_effort)
# All Items
sapply(df_eq_pasa_control_cor , is.factor)## PASA_depends_experts PASA_protect PASA_determine
## FALSE FALSE FALSE
## PASA_my_effort
## FALSEcor <- round(cor(df_eq_pasa_control_cor [sapply(df_eq_pasa_control_cor , function(x) !is.factor(x))], use = "complete.obs"),2) #Correlation matrix of indices
#cor
library(reshape2)
melted_cormat <- melt(cor)
#head(melted_cormat)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
cormat[upper.tri(cormat)] <- NA
return(cormat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
cormat[lower.tri(cormat)]<- NA
return(cormat)
}
upper_tri <- get_upper_tri(cor)
#upper_tri
# Melt the correlation matrix
library(reshape2)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cor){
# Use correlation between variables as distance
dd <- as.dist((1-cor)/2)
hc <- hclust(dd)
cor <-cor[hc$order, hc$order]
}
# Reorder the correlation matrix
cor <- reorder_cormat(cor)
upper_tri <- get_upper_tri(cor)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
library(ggplot2)
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
cor_heatmap <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(1, 0),
legend.position = c(0.6, 0.7),
legend.direction = "horizontal")+
guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
title.position = "top", title.hjust = 0.5))
cor_heatmap
# Ability items
df_eq_pasa_ability_cor <- dplyr::select(df_eq_pasa, PASA_I_can, PASA_no_idea_r, PASA_solutions, PASA_can_solve)
# All Items
sapply(df_eq_pasa_ability_cor, is.factor)## PASA_I_can PASA_no_idea_r PASA_solutions PASA_can_solve
## FALSE FALSE FALSE FALSEcor <- round(cor(df_eq_pasa_ability_cor [sapply(df_eq_pasa_ability_cor, function(x) !is.factor(x))], use = "complete.obs"),2) #Correlation matrix of indices
#cor
library(reshape2)
melted_cormat <- melt(cor)
#head(melted_cormat)
# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
cormat[upper.tri(cormat)] <- NA
return(cormat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
cormat[lower.tri(cormat)]<- NA
return(cormat)
}
upper_tri <- get_upper_tri(cor)
#upper_tri
# Melt the correlation matrix
library(reshape2)
melted_cormat <- melt(upper_tri, na.rm = TRUE)
reorder_cormat <- function(cor){
# Use correlation between variables as distance
dd <- as.dist((1-cor)/2)
hc <- hclust(dd)
cor <-cor[hc$order, hc$order]
}
# Reorder the correlation matrix
cor <- reorder_cormat(cor)
upper_tri <- get_upper_tri(cor)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)
# Create a ggheatmap
library(ggplot2)
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
geom_tile(color = "white")+
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limit = c(-1,1), space = "Lab",
name="Pearson\nCorrelation") +
theme_minimal()+ # minimal theme
theme(axis.text.x = element_text(angle = 45, vjust = 1,
size = 12, hjust = 1))+
coord_fixed()
cor_heatmap <- ggheatmap +
geom_text(aes(Var2, Var1, label = value), color = "black", size = 4) +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
panel.background = element_blank(),
axis.ticks = element_blank(),
legend.justification = c(1, 0),
legend.position = c(0.6, 0.7),
legend.direction = "horizontal")+
guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
title.position = "top", title.hjust = 0.5))
cor_heatmap
#Creating composite scores for PASA measure items - Four anticipatory cognitive appraisal processes: Threat Items: PASA_no_threat_r + PASA_unpleasant + PASA_no_worry_r + PASA_scared - Higher scores reflect higher perceptions of task-as-threat
Challenge Items PASA_task_import + PASA_do_not_care_r + PASA_not_challenge_r + PASA_not_challenge_r #Second is repeated to reflect missing item - Higher scores reflect higher perceptions of task-as-challenge
Self-perception of Ability Items PASA_I_can + PASA_no_idea_r + PASA_solutions + PASA_can_solve - Higher scores reflect higher self-perceptions of ability.
Control Belief Items: PASA_depends_experts + PASA_protect + PASA_determine + PASA_my_effort - Higher scores reflect higher self-perceptions of control.
Composite for Threat Items
#Checking for internal validity
# Note: adding PASA_not_challenge_r because there is internal validity (though previous recommendations suggest to put this in the "challenge" item)
df_eq_pasa_threat <- dplyr::select(df_eq_pasa, PASA_no_threat_r, PASA_unpleasant, PASA_no_worry_r, PASA_scared, PASA_not_challenge_r) #should I be including "PID, Condition" for this function?
psych::alpha(df_eq_pasa_threat) #calculating chronbach's alpha for threat items##
## Reliability analysis
## Call: psych::alpha(x = df_eq_pasa_threat)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.87 0.87 0.86 0.58 6.9 0.019 3.1 1.1 0.56
##
## lower alpha upper 95% confidence boundaries
## 0.83 0.87 0.91
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## PASA_no_threat_r 0.86 0.86 0.84 0.60 6.1 0.023 0.0064
## PASA_unpleasant 0.85 0.85 0.82 0.59 5.8 0.023 0.0023
## PASA_no_worry_r 0.84 0.84 0.82 0.57 5.4 0.025 0.0086
## PASA_scared 0.83 0.83 0.80 0.55 5.0 0.026 0.0053
## PASA_not_challenge_r 0.85 0.85 0.83 0.58 5.5 0.024 0.0092
## med.r
## PASA_no_threat_r 0.59
## PASA_unpleasant 0.59
## PASA_no_worry_r 0.56
## PASA_scared 0.54
## PASA_not_challenge_r 0.56
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PASA_no_threat_r 109 0.79 0.78 0.70 0.66 3.5 1.4
## PASA_unpleasant 109 0.80 0.80 0.74 0.68 2.7 1.3
## PASA_no_worry_r 109 0.83 0.82 0.77 0.71 3.6 1.4
## PASA_scared 109 0.86 0.85 0.82 0.76 2.9 1.4
## PASA_not_challenge_r 109 0.81 0.82 0.75 0.71 2.7 1.2
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## PASA_no_threat_r 0.05 0.24 0.26 0.17 0.23 0.06 0
## PASA_unpleasant 0.20 0.28 0.24 0.15 0.12 0.01 0
## PASA_no_worry_r 0.05 0.21 0.24 0.25 0.15 0.11 0
## PASA_scared 0.17 0.28 0.27 0.14 0.12 0.04 0
## PASA_not_challenge_r 0.14 0.37 0.26 0.17 0.05 0.02 0alpha_threat <- psych::alpha(df_eq_pasa_threat)$total$std.alpha #extracting just the alpha value for threat items
alpha_threat <- round(alpha_threat, 2)
View(alpha_threat) #Good internal consistency
#Creating composite for Threat Items
df_eq_pasa$pasa_threat <- (df_eq_pasa$PASA_no_threat_r + df_eq_pasa$PASA_unpleasant + df_eq_pasa$PASA_no_worry_r + df_eq_pasa$PASA_scared + df_eq_pasa$PASA_not_challenge_r)/5 #Composite Score - Threat Items
#FUTURE TO DO: Standardizing the composite scores.Composite for Challenge Items
#Checking internal consistency
df_eq_pasa_challenge <- dplyr::select(df_eq_pasa, PASA_task_import, PASA_do_not_care_r, PASA_not_challenge_r) #should I be including "PID, Condition" for this function?
psych::alpha(df_eq_pasa_challenge) #calculating chronbach's alpha ##
## Reliability analysis
## Call: psych::alpha(x = df_eq_pasa_challenge)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.56 0.55 0.51 0.29 1.2 0.072 3.2 0.9 0.17
##
## lower alpha upper 95% confidence boundaries
## 0.42 0.56 0.7
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## PASA_task_import 0.29 0.29 0.17 0.17 0.42 0.135 NA
## PASA_do_not_care_r 0.27 0.27 0.16 0.16 0.37 0.139 NA
## PASA_not_challenge_r 0.71 0.71 0.55 0.55 2.47 0.055 NA
## med.r
## PASA_task_import 0.17
## PASA_do_not_care_r 0.16
## PASA_not_challenge_r 0.55
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PASA_task_import 109 0.80 0.78 0.65 0.47 3.6 1.3
## PASA_do_not_care_r 109 0.80 0.79 0.67 0.49 3.2 1.2
## PASA_not_challenge_r 109 0.58 0.61 0.23 0.19 2.7 1.2
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## PASA_task_import 0.02 0.19 0.29 0.22 0.18 0.09 0
## PASA_do_not_care_r 0.06 0.21 0.39 0.19 0.09 0.06 0
## PASA_not_challenge_r 0.14 0.37 0.26 0.17 0.05 0.02 0alpha_challenge <- psych::alpha(df_eq_pasa_challenge)$total$std.alpha #extracting just the alpha value
alpha_challenge <- round(alpha_challenge, 2)
View(alpha_challenge) #Poor internal consistency
#Creating composite measure
df_eq_pasa$pasa_challenge <- (df_eq_pasa$PASA_task_import + df_eq_pasa$PASA_do_not_care_r + df_eq_pasa$PASA_not_challenge_r)/3 #Composite Score - Challenge Items. NOTE: Item #14 was missing from the measure provided to the participants. This composite would typically be the average of the four items, but is now composed of only 3 items.
#NEED TO CONSIDER HOW TO HANDLE THIS POOR ALPHAComposite for Ability Items
#Checking internal consistency
df_eq_pasa_ability <- dplyr::select(df_eq_pasa, PASA_I_can, PASA_no_idea_r, PASA_solutions, PASA_can_solve) #should I be including "PID, Condition" for this function?
psych::alpha(df_eq_pasa_ability) #calculating chronbach's alpha ##
## Reliability analysis
## Call: psych::alpha(x = df_eq_pasa_ability)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.74 0.76 0.71 0.45 3.2 0.039 3 0.81 0.44
##
## lower alpha upper 95% confidence boundaries
## 0.67 0.74 0.82
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
## PASA_I_can 0.63 0.65 0.56 0.39 1.9 0.059 0.00073 0.39
## PASA_no_idea_r 0.73 0.73 0.65 0.48 2.7 0.044 0.00584 0.46
## PASA_solutions 0.71 0.74 0.66 0.48 2.8 0.047 0.00715 0.49
## PASA_can_solve 0.68 0.70 0.62 0.44 2.3 0.051 0.00518 0.46
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PASA_I_can 109 0.81 0.82 0.76 0.65 2.9 0.98
## PASA_no_idea_r 109 0.80 0.73 0.59 0.51 3.1 1.40
## PASA_solutions 109 0.69 0.73 0.58 0.50 2.7 0.89
## PASA_can_solve 109 0.75 0.77 0.66 0.56 3.3 0.95
##
## Non missing response frequency for each item
## 1 2 3 4 5 6 miss
## PASA_I_can 0.06 0.29 0.46 0.14 0.04 0.02 0
## PASA_no_idea_r 0.14 0.25 0.26 0.20 0.09 0.06 0
## PASA_solutions 0.06 0.37 0.47 0.06 0.04 0.01 0
## PASA_can_solve 0.02 0.18 0.40 0.30 0.08 0.01 0alpha_ability <- psych::alpha(df_eq_pasa_ability)$total$std.alpha #extracting just the alpha value
alpha_ability <- round(alpha_ability, 2)
View(alpha_ability) #Acceptable internal consistency
#Creating composite measure
df_eq_pasa$pasa_ability <- (df_eq_pasa$PASA_I_can + df_eq_pasa$PASA_no_idea_r + df_eq_pasa$PASA_solutions + df_eq_pasa$PASA_can_solve)/4 #Composite Score - Self-perception of Ability ItemsComposite for Control Belief Items
#Checking internal consistency
df_eq_pasa_control <- dplyr::select(df_eq_pasa, PASA_depends_experts, PASA_protect, PASA_determine, PASA_my_effort) #should I be including "PID, Condition" for this function?
df_eq_pasa_control <- dplyr::select(df_eq_pasa, PASA_depends_experts, PASA_protect, PASA_determine, PASA_my_effort) #should I be including "PID, Condition" for this function?
df_eq_pasa_control <- dplyr::select(df_eq_pasa, PASA_depends_experts, PASA_protect, PASA_determine, PASA_my_effort) #should I be including "PID, Condition" for this function?
df_eq_pasa_control <- dplyr::select(df_eq_pasa, PASA_depends_experts, PASA_protect, PASA_determine, PASA_my_effort) #should I be including "PID, Condition" for this function?
df_eq_pasa_control <- dplyr::select(df_eq_pasa, PASA_depends_experts, PASA_protect, PASA_determine, PASA_my_effort) #should I be including "PID, Condition" for this function?
psych::alpha(df_eq_pasa_control) #calculating chronbach's alpha ##
## Reliability analysis
## Call: psych::alpha(x = df_eq_pasa_control)
##
## raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
## 0.28 0.29 0.26 0.091 0.4 0.11 2.7 0.6 0.079
##
## lower alpha upper 95% confidence boundaries
## 0.06 0.28 0.5
##
## Reliability if an item is dropped:
## raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r
## PASA_depends_experts 0.31 0.30 0.24 0.127 0.44 0.11 0.0094
## PASA_protect 0.22 0.22 0.19 0.088 0.29 0.13 0.0204
## PASA_determine 0.23 0.23 0.17 0.088 0.29 0.13 0.0103
## PASA_my_effort 0.14 0.16 0.14 0.061 0.19 0.14 0.0225
## med.r
## PASA_depends_experts 0.128
## PASA_protect 0.128
## PASA_determine 0.030
## PASA_my_effort 0.029
##
## Item statistics
## n raw.r std.r r.cor r.drop mean sd
## PASA_depends_experts 109 0.56 0.52 0.19 0.083 3.1 1.17
## PASA_protect 109 0.53 0.57 0.29 0.149 2.4 0.98
## PASA_determine 109 0.56 0.57 0.31 0.134 2.7 1.08
## PASA_my_effort 109 0.60 0.60 0.37 0.208 2.6 1.05
##
## Non missing response frequency for each item
## 1 2 2.55555555555556 3 4 5 6 miss
## PASA_depends_experts 0.07 0.22 0.00 0.39 0.20 0.08 0.04 0
## PASA_protect 0.15 0.41 0.00 0.33 0.06 0.05 0.00 0
## PASA_determine 0.13 0.29 0.00 0.37 0.16 0.05 0.01 0
## PASA_my_effort 0.17 0.33 0.01 0.31 0.15 0.04 0.00 0alpha_control <- psych::alpha(df_eq_pasa_control)$total$std.alpha #extracting just the alpha value
alpha_control <- round(alpha_control, 2)
View(alpha_control) #Unacceptable internal consistency
#Creating composite measure
df_eq_pasa$pasa_control <- (df_eq_pasa$PASA_depends_experts + df_eq_pasa$PASA_protect + df_eq_pasa$PASA_determine + df_eq_pasa$PASA_my_effort)/4 #Composite Score - Control Belief Items
#NEED TO CONSIDER HOW TO HANDLE THIS UNACCEPTABLE ALPHAExploratory Factor Analysis
The following is the suggested method for creating composite scores for PASA measure items - Primary Appraisal, Secondary Appraisal, and Tertiary Appraisal. However, we did not find reliability of scale items for this method.
- Primary Appraisal (threat + challenge / 2): a person’s judgment about the significance of an event as stressful, positive, controllable, challenging, or irrelevant.
- Higher scores reflect higher appraisal of significance
2.Secondary Appraisal (self-concept + control belief / 2): a person’s assessment of available coping resources and options when faced with a stressor - Higher scores reflect higher perception of coping resources
3.Tertiary Appraisal, aka “global PASA scale” / “stress index”. (stress index = Primary Appraisal - Secondary Appraisal): combines “primary appraisal” and “secondary appraisal” scales to provide an integrated measure of transactional stress perception. - Higher score represents a higher general stress perception.
We will analyze the composite measure as proposed by the measure creators, but because of the low psychometrics of the Challenge Items (Chronbach’s alpha = 0.55) and Control Belief Item (Chronbach’s alpha = 0.29), we will look at the the following subscales based on high inter-item reliability: 1. Threat Item (adding the PASA_not_challenge_r item)(Chronbach’s alpha = 0.87) 2. Ability Item (Chronbach’s alpha = 0.76) 3. Resilience Appraisal Item: Ability item - Threat Item
#df_eq_pasa$pasa_primary <- (df_eq_pasa$pasa_threat + df_eq_pasa$pasa_challenge)/2 # Creating composite for Primary Appraisal (threat + challenge / 2)
#df_eq_pasa$pasa_primary <- round(df_eq_pasa$pasa_primary, 2) #rounding to hundredths
#df_eq_pasa$pasa_secondary <- ( + df_eq_pasa$pasa_control) /2 # Creating composite for Secondary Appraisal (self-concept + control belief / 2)
#df_eq_pasa$pasa_secondary <- round(df_eq_pasa$pasa_secondary, 2) #rounding to hundredths
#df_eq_pasa$pasa_stress <- (df_eq_pasa$pasa_primary-df_eq_pasa$pasa_secondary) # Creating composite for Tertiary Appraisal (stress index = Primary Appraisal - Secondary Appraisal)
#df_eq_pasa$pasa_stress <- round(df_eq_pasa$pasa_stress, 2) #rounding to hundredths
#...Instead, we will use only
# Construing Item as Stress Appraisal
df_eq_pasa$pasa_stress <- (df_eq_pasa$pasa_threat-df_eq_pasa$pasa_ability) # Creating composite for Tertiary Appraisal (stress index = Threat Appraisal - Ability Appraisal)
df_eq_pasa$pasa_stress <- round(df_eq_pasa$pasa_stress, 2) #rounding to hundredths
# Construing Item as Resilience Appraisal
df_eq_pasa$pasa_resilience <- (df_eq_pasa$pasa_ability-df_eq_pasa$pasa_threat) # Creating composite for Resilience Appraisal Item (stress index = Threat Appraisal - Ability Appraisal)
df_eq_pasa$pasa_resilience <- round(df_eq_pasa$pasa_resilience, 2) #rounding to hundredthsdf_eq_pasa <- dplyr::select(df_eq_pasa, PID, Condition, pasa_threat, pasa_ability,
#pasa_challenge, pasa_control, pasa_primary, pasa_secondary, #not including bc of internal validity issues
pasa_stress, pasa_resilience) #Subsetting Data to include only composite items
df_eq_pasa <- df_eq_pasa %>%
mutate(group_fct = factor(Condition,
labels = c("Group 1", "Group 2"))) #Recoding Condition as a factor
data(df_eq_pasa)## Warning in data(df_eq_pasa): data set 'df_eq_pasa' not foundhead(df_eq_pasa)## PID Condition pasa_threat pasa_ability pasa_stress pasa_resilience group_fct
## 1 103 2 3.0 3.50 -0.50 0.50 Group 2
## 2 135 2 4.4 2.50 1.90 -1.90 Group 2
## 3 136 2 1.8 4.50 -2.70 2.70 Group 2
## 4 165 1 1.0 5.50 -4.50 4.50 Group 1
## 5 166 1 1.6 4.25 -2.65 2.65 Group 1
## 6 167 1 3.0 2.75 0.25 -0.25 Group 1Describe the data
psych::describeBy(df_eq_pasa$pasa_threat, df_eq_pasa$Condition, mat = TRUE) #- Higher score represents a higher general threat perception.## item group1 vars n mean sd median trimmed mad min max range
## X11 1 1 1 56 3.064286 1.123052 3.0 3.043478 1.48260 1.0 5.4 4.4
## X12 2 2 1 53 3.071698 1.035593 2.8 3.009302 0.88956 1.6 5.8 4.2
## skew kurtosis se
## X11 0.1682140 -1.066960 0.1500742
## X12 0.5853783 -0.392869 0.1422496psych::describeBy(df_eq_pasa$pasa_ability, df_eq_pasa$Condition, mat = TRUE) #- Higher score represents a higher general ability perception.## item group1 vars n mean sd median trimmed mad min max
## X11 1 1 1 56 3.165179 0.8558313 3.125 3.146739 0.926625 1.25 5.5
## X12 2 2 1 53 2.764151 0.7026993 2.750 2.750000 0.741300 1.25 4.5
## range skew kurtosis se
## X11 4.25 0.3032316 0.41564881 0.11436527
## X12 3.25 0.2695290 0.08905955 0.09652317psych::describeBy(df_eq_pasa$pasa_stress, df_eq_pasa$Condition, mat = TRUE) #- Higher score represents a higher general stress perception.## item group1 vars n mean sd median trimmed mad min
## X11 1 1 1 56 -0.1008929 1.756766 0.175 -0.07717391 1.92738 -4.5
## X12 2 2 1 53 0.3075472 1.523861 0.250 0.26162791 1.55673 -2.7
## max range skew kurtosis se
## X11 3.65 8.15 -0.1034715 -0.56741931 0.2347577
## X12 4.55 7.25 0.4527281 0.03518629 0.2093185psych::describeBy(df_eq_pasa$pasa_resilience, df_eq_pasa$Condition, mat = TRUE) #- Higher score represents a higher general resilience perception.## item group1 vars n mean sd median trimmed mad min
## X11 1 1 1 56 0.1008929 1.756766 -0.175 0.07717391 1.92738 -3.65
## X12 2 2 1 53 -0.3075472 1.523861 -0.250 -0.26162791 1.55673 -4.55
## max range skew kurtosis se
## X11 4.5 8.15 0.1034715 -0.56741931 0.2347577
## X12 2.7 7.25 -0.4527281 0.03518629 0.2093185Visualizing Overall Threat Scores
df_eq_pasa_1 <- subset(df_eq_pasa, Condition > 1) #Subset of Condition 1 only
df_eq_pasa_2 <- subset(df_eq_pasa, Condition < 2) #Subset of Condition 2 only
library(ggplot2)
hist_1 <- df_eq_pasa_1 %>%
ggplot(aes(x = pasa_threat)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_1 #Histogram of Threat Index Scores for Condition 1
hist_2 <- df_eq_pasa_2 %>%
ggplot(aes(x = pasa_threat)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_2 #Histogram of Threat Index Scores for Condition 2
hist <- df_eq_pasa %>%
ggplot(aes(x = pasa_threat, fill = group_fct)) +
geom_density(lwd = 1, alpha = .3)
hist #Histogram of Threat Index Scores for both conditions
Visualizing Overall Ability Scores
df_eq_pasa_1 <- subset(df_eq_pasa, Condition > 1) #Subset of Condition 1 only
df_eq_pasa_2 <- subset(df_eq_pasa, Condition < 2) #Subset of Condition 2 only
library(ggplot2)
hist_1 <- df_eq_pasa_1 %>%
ggplot(aes(x = pasa_ability)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_1 #Histogram of Ability Index Scores for Condition 1
hist_2 <- df_eq_pasa_2 %>%
ggplot(aes(x = pasa_ability)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_2 #Histogram of Ability Index Scores for Condition 2
hist <- df_eq_pasa %>%
ggplot(aes(x = pasa_ability, fill = group_fct)) +
geom_density(lwd = 1, alpha = .3)
hist #Histogram of Ability Index Scores for both conditions
Visualizing Overall Stress Scores
df_eq_pasa_1 <- subset(df_eq_pasa, Condition > 1) #Subset of Condition 1 only
df_eq_pasa_2 <- subset(df_eq_pasa, Condition < 2) #Subset of Condition 2 only
library(ggplot2)
hist_1 <- df_eq_pasa_1 %>%
ggplot(aes(x = pasa_stress)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_1 #Histogram of Stress Index Scores for Condition 1
hist_2 <- df_eq_pasa_2 %>%
ggplot(aes(x = pasa_stress)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_2 #Histogram of Stress Index Scores for Condition 2
hist <- df_eq_pasa %>%
ggplot(aes(x = pasa_stress, fill = group_fct)) +
geom_density(lwd = 1, alpha = .3)
hist #Histogram of Stress Index Scores for both conditions
Visualizing Overall Resilience Scores
df_eq_pasa_1 <- subset(df_eq_pasa, Condition > 1) #Subset of Condition 1 only
df_eq_pasa_2 <- subset(df_eq_pasa, Condition < 2) #Subset of Condition 2 only
library(ggplot2)
hist_1 <- df_eq_pasa_1 %>%
ggplot(aes(x = pasa_resilience)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_1 #Histogram of resilience Index Scores for Condition 1
hist_2 <- df_eq_pasa_2 %>%
ggplot(aes(x = pasa_resilience)) +
geom_histogram(fill = "turquoise",
colour = "black",
alpha = .6,
bins = 30) +
geom_density(colour = "darkorchid",
lwd = 1)
hist_2 #Histogram of resilience Index Scores for Condition 2
hist <- df_eq_pasa %>%
ggplot(aes(x = pasa_resilience, fill = group_fct)) +
geom_density(lwd = 1, alpha = .3)
hist #Histogram of resilience Index Scores for both conditions
#Checking the Assumptions
#(grab code from class)pasa_box <- boxplot(pasa_stress~Condition,data=df_eq_pasa, main="Stress Index Scores by Condition",
xlab="Condition", ylab="Stress index") #Boxplot for Stress Index by condition
pasa_box ## $stats
## [,1] [,2]
## [1,] -4.500 -2.70
## [2,] -1.400 -0.80
## [3,] 0.175 0.25
## [4,] 1.200 1.30
## [5,] 3.650 4.15
##
## $n
## [1] 56 53
##
## $conf
## [,1] [,2]
## [1,] -0.3739546 -0.2057624
## [2,] 0.7239546 0.7057624
##
## $out
## [1] 4.55
##
## $group
## [1] 2
##
## $names
## [1] "1" "2"#Detecting Outliers
#IF outliers had been detected: A common method for removing outliers that are Q1 (first quartile of the variable distribution) − 1.5 × IQR (IQR=Q3−Q1; the interquartile range). Remove outliers that are above Q3 (third quartile of the variable distribution) + 1.5 × IQR (IQR=Q3−Q1; the interquartile range). If you remove the outliers, report with and without. Option 2 = windsorize - aka collapse them down.
#quantile(df_eq_pasa$pasa_stress, na.rm=TRUE)
#IQR <- .48*2Threat Index
pasa_box_threat <- boxplot(pasa_threat~Condition,data=df_eq_pasa, main="Threat Index Scores by Condition",
xlab="Condition", ylab="Threat index") #Boxplot for Threat Index by condition
pasa_box_threat## $stats
## [,1] [,2]
## [1,] 1.0 1.6
## [2,] 2.2 2.2
## [3,] 3.0 2.8
## [4,] 4.0 3.8
## [5,] 5.4 5.8
##
## $n
## [1] 56 53
##
## $conf
## [,1] [,2]
## [1,] 2.619955 2.452752
## [2,] 3.380045 3.147248
##
## $out
## numeric(0)
##
## $group
## numeric(0)
##
## $names
## [1] "1" "2"Running the regression of Threat Index
# load required packages
library(lme4) # fitting models## Loading required package: Matrixlibrary(lmerTest)##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## steplibrary(sjPlot) # table functions
library(sjmisc) # sample data## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!model_ThreatIndex <- lm(pasa_threat ~ Condition, # outcome (Threat Index) ~ predictor (Condition)
data = df_eq_pasa)
summary(model_ThreatIndex)##
## Call:
## lm(formula = pasa_threat ~ Condition, data = df_eq_pasa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0643 -0.8643 -0.0717 0.9283 2.7283
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.056873 0.324964 9.407 1.12e-15 ***
## Condition 0.007412 0.207243 0.036 0.972
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.081 on 107 degrees of freedom
## Multiple R-squared: 1.196e-05, Adjusted R-squared: -0.009334
## F-statistic: 0.001279 on 1 and 107 DF, p-value: 0.9715str(model_ThreatIndex)## List of 12
## $ coefficients : Named num [1:2] 3.05687 0.00741
## ..- attr(*, "names")= chr [1:2] "(Intercept)" "Condition"
## $ residuals : Named num [1:109] -0.0717 1.3283 -1.2717 -2.0643 -1.4643 ...
## ..- attr(*, "names")= chr [1:109] "1" "2" "3" "4" ...
## $ effects : Named num [1:109] -32.0297 -0.0387 -1.3755 -1.9334 -1.3334 ...
## ..- attr(*, "names")= chr [1:109] "(Intercept)" "Condition" "" "" ...
## $ rank : int 2
## $ fitted.values: Named num [1:109] 3.07 3.07 3.07 3.06 3.06 ...
## ..- attr(*, "names")= chr [1:109] "1" "2" "3" "4" ...
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:109, 1:2] -10.4403 0.0958 0.0958 0.0958 0.0958 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:109] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "(Intercept)" "Condition"
## .. ..- attr(*, "assign")= int [1:2] 0 1
## ..$ qraux: num [1:2] 1.1 1.09
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 107
## $ xlevels : Named list()
## $ call : language lm(formula = pasa_threat ~ Condition, data = df_eq_pasa)
## $ terms :Classes 'terms', 'formula' language pasa_threat ~ Condition
## .. ..- attr(*, "variables")= language list(pasa_threat, Condition)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "pasa_threat" "Condition"
## .. .. .. ..$ : chr "Condition"
## .. ..- attr(*, "term.labels")= chr "Condition"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(pasa_threat, Condition)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
## .. .. ..- attr(*, "names")= chr [1:2] "pasa_threat" "Condition"
## $ model :'data.frame': 109 obs. of 2 variables:
## ..$ pasa_threat: num [1:109] 3 4.4 1.8 1 1.6 3 1.8 3.4 2.6 1.6 ...
## ..$ Condition : int [1:109] 2 2 2 1 1 1 2 1 2 1 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language pasa_threat ~ Condition
## .. .. ..- attr(*, "variables")= language list(pasa_threat, Condition)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "pasa_threat" "Condition"
## .. .. .. .. ..$ : chr "Condition"
## .. .. ..- attr(*, "term.labels")= chr "Condition"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(pasa_threat, Condition)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
## .. .. .. ..- attr(*, "names")= chr [1:2] "pasa_threat" "Condition"
## - attr(*, "class")= chr "lm"broom::tidy(model_ThreatIndex)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.06 0.325 9.41 1.12e-15
## 2 Condition 0.00741 0.207 0.0358 9.72e- 1tab_model(model_ThreatIndex)| pasa threat | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 3.06 | 2.41 – 3.70 | <0.001 |
| Condition | 0.01 | -0.40 – 0.42 | 0.972 |
| Observations | 109 | ||
| R2 / R2 adjusted | 0.000 / -0.009 | ||
# Threat Index model results
#model_ThreatIndex_results <- apa_print(model_ThreatIndex)
#model_ThreatIndex_results$table %>%
#apa_table(caption = "Model Regressing Threat Index Score on Condition",
# note = "* p < 0.05; ** p < 0.01; *** p < 0.001",
# escape = TRUE)
# Get r squared
library(broom)
model_info = augment(model_ThreatIndex)
head(model_info)## # A tibble: 6 x 9
## pasa_threat Condition .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
## <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3 2 3.07 0.149 -0.0717 0.0189 1.09 4.31e-5 -0.0669
## 2 4.4 2 3.07 0.149 1.33 0.0189 1.08 1.48e-2 1.24
## 3 1.8 2 3.07 0.149 -1.27 0.0189 1.08 1.36e-2 -1.19
## 4 1 1 3.06 0.145 -2.06 0.0179 1.07 3.37e-2 -1.93
## 5 1.6 1 3.06 0.145 -1.46 0.0179 1.08 1.70e-2 -1.37
## 6 3 1 3.06 0.145 -0.0643 0.0179 1.09 3.27e-5 -0.0600summary(model_ThreatIndex)$r.squared## [1] 1.195548e-05Write-Up of Threat Index
A simple linear regression was performed to predict threat appraisal based on condition. The mean threat appraisal score for Condition 1 (M = 3.064, SE = 0.150) was not significantly different than Condition 2 (M = 3.072, SE = 0.142) (F(107) = 0.001, p = 0.972), with an R2 of 0.000.
Ability Index
pasa_box_ability <- boxplot(pasa_ability~Condition,data=df_eq_pasa, main="Ability Index Scores by Condition",
xlab="Condition", ylab="Ability index") #Boxplot for Ability Index by condition
pasa_box_ability## $stats
## [,1] [,2]
## [1,] 1.250 1.25
## [2,] 2.625 2.25
## [3,] 3.125 2.75
## [4,] 3.750 3.25
## [5,] 4.500 4.50
##
## $n
## [1] 56 53
##
## $conf
## [,1] [,2]
## [1,] 2.887472 2.53297
## [2,] 3.362528 2.96703
##
## $out
## [1] 5.5 5.5
##
## $group
## [1] 1 1
##
## $names
## [1] "1" "2"Running the regression of Self-perceived Ability Index
model_AbilityIndex <- lm(pasa_ability ~ group_fct, # outcome (Ability Index) ~ predictor (Condition)
data = df_eq_pasa)
summary(model_AbilityIndex)##
## Call:
## lm(formula = pasa_ability ~ group_fct, data = df_eq_pasa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.91518 -0.51415 -0.01415 0.48585 2.33482
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1652 0.1049 30.167 < 2e-16 ***
## group_fctGroup 2 -0.4010 0.1505 -2.665 0.00888 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7852 on 107 degrees of freedom
## Multiple R-squared: 0.06226, Adjusted R-squared: 0.05349
## F-statistic: 7.104 on 1 and 107 DF, p-value: 0.008885str(model_AbilityIndex)## List of 13
## $ coefficients : Named num [1:2] 3.165 -0.401
## ..- attr(*, "names")= chr [1:2] "(Intercept)" "group_fctGroup 2"
## $ residuals : Named num [1:109] 0.736 -0.264 1.736 2.335 1.085 ...
## ..- attr(*, "names")= chr [1:109] "1" "2" "3" "4" ...
## $ effects : Named num [1:109] -31.01 2.09 1.7 2.24 0.99 ...
## ..- attr(*, "names")= chr [1:109] "(Intercept)" "group_fctGroup 2" "" "" ...
## $ rank : int 2
## $ fitted.values: Named num [1:109] 2.76 2.76 2.76 3.17 3.17 ...
## ..- attr(*, "names")= chr [1:109] "1" "2" "3" "4" ...
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:109, 1:2] -10.4403 0.0958 0.0958 0.0958 0.0958 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:109] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "(Intercept)" "group_fctGroup 2"
## .. ..- attr(*, "assign")= int [1:2] 0 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ group_fct: chr "contr.treatment"
## ..$ qraux: num [1:2] 1.1 1.09
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 107
## $ contrasts :List of 1
## ..$ group_fct: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ group_fct: chr [1:2] "Group 1" "Group 2"
## $ call : language lm(formula = pasa_ability ~ group_fct, data = df_eq_pasa)
## $ terms :Classes 'terms', 'formula' language pasa_ability ~ group_fct
## .. ..- attr(*, "variables")= language list(pasa_ability, group_fct)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "pasa_ability" "group_fct"
## .. .. .. ..$ : chr "group_fct"
## .. ..- attr(*, "term.labels")= chr "group_fct"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(pasa_ability, group_fct)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## .. .. ..- attr(*, "names")= chr [1:2] "pasa_ability" "group_fct"
## $ model :'data.frame': 109 obs. of 2 variables:
## ..$ pasa_ability: num [1:109] 3.5 2.5 4.5 5.5 4.25 2.75 3 3.5 2.75 3 ...
## ..$ group_fct : Factor w/ 2 levels "Group 1","Group 2": 2 2 2 1 1 1 2 1 2 1 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language pasa_ability ~ group_fct
## .. .. ..- attr(*, "variables")= language list(pasa_ability, group_fct)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "pasa_ability" "group_fct"
## .. .. .. .. ..$ : chr "group_fct"
## .. .. ..- attr(*, "term.labels")= chr "group_fct"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(pasa_ability, group_fct)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## .. .. .. ..- attr(*, "names")= chr [1:2] "pasa_ability" "group_fct"
## - attr(*, "class")= chr "lm"broom::tidy(model_AbilityIndex)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.17 0.105 30.2 3.88e-54
## 2 group_fctGroup 2 -0.401 0.150 -2.67 8.88e- 3tab_model(model_AbilityIndex)| pasa ability | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 3.17 | 2.96 – 3.37 | <0.001 |
| group_fct [Group 2] | -0.40 | -0.70 – -0.10 | 0.009 |
| Observations | 109 | ||
| R2 / R2 adjusted | 0.062 / 0.053 | ||
# Ability Index model results
#model_AbilityIndex_results <- apa_print(model_AbilityIndex)
#Calculate a cohen's D effect size (to compare to multiple R squared below).
library(lsr)
cohensD(pasa_ability ~ group_fct, data = df_eq_pasa)## [1] 0.5107648#This (.51) is a medium effect size. The multiple R squared is .06. Corre of about .25...medium Write-Up of Ability Index
A simple linear regression was performed to predict self-perceived ability based on condition. The mean self-perceived ability score for Condition 1 (M = 3.165, SE = 0.114) was significantly different than Condition 2 (M = 2.764 , SE = 0.097) (F(107) = 7.104, p = 0.009), with an R2 of 0.062. Removing outliers did not impact this effect (F(105) = 5.084, p = 0.026), with an R2 of 0.046.
Detecting Outliers
#Obtaining values of outliers
library(dplyr)
library(ggplot2)
is_outlier <- function(x) {
return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}
df_eq_pasa %>%
group_by(Condition) %>%
mutate(outlier = ifelse(is_outlier(pasa_ability), pasa_ability, as.numeric(NA))) %>%
ggplot(., aes(x = factor(Condition), y = pasa_ability)) +
geom_boxplot() +
geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)
# Removing Outliers
df_eq_pasa_out1 <- df_eq_pasa %>%
filter(Condition == 1, pasa_ability < 5.4)
df_eq_pasa_out2 <- df_eq_pasa %>%
filter(Condition == 2, pasa_ability > 0)
df_eq_pasa_out <- rbind(df_eq_pasa_out1, df_eq_pasa_out2)Re-running the regression of Self-perceived Ability Index (outliers removed)
model_AbilityIndex_out <- lm(pasa_ability ~ group_fct, # outcome (Ability Index) ~ predictor (Condition)
data = df_eq_pasa_out)
summary(model_AbilityIndex_out)##
## Call:
## lm(formula = pasa_ability ~ group_fct, data = df_eq_pasa_out)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.82870 -0.51415 -0.01415 0.48585 1.73585
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.07870 0.09818 31.358 <2e-16 ***
## group_fctGroup 2 -0.31455 0.13950 -2.255 0.0262 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7215 on 105 degrees of freedom
## Multiple R-squared: 0.04619, Adjusted R-squared: 0.0371
## F-statistic: 5.084 on 1 and 105 DF, p-value: 0.02622str(model_AbilityIndex_out)## List of 13
## $ coefficients : Named num [1:2] 3.079 -0.315
## ..- attr(*, "names")= chr [1:2] "(Intercept)" "group_fctGroup 2"
## $ residuals : Named num [1:107] 1.1713 -0.3287 0.4213 -0.0787 0.6713 ...
## ..- attr(*, "names")= chr [1:107] "1" "2" "3" "4" ...
## $ effects : Named num [1:107] -30.235 -1.627 0.353 -0.147 0.603 ...
## ..- attr(*, "names")= chr [1:107] "(Intercept)" "group_fctGroup 2" "" "" ...
## $ rank : int 2
## $ fitted.values: Named num [1:107] 3.08 3.08 3.08 3.08 3.08 ...
## ..- attr(*, "names")= chr [1:107] "1" "2" "3" "4" ...
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:107, 1:2] -10.3441 0.0967 0.0967 0.0967 0.0967 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:107] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "(Intercept)" "group_fctGroup 2"
## .. ..- attr(*, "assign")= int [1:2] 0 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ group_fct: chr "contr.treatment"
## ..$ qraux: num [1:2] 1.1 1.09
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 105
## $ contrasts :List of 1
## ..$ group_fct: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ group_fct: chr [1:2] "Group 1" "Group 2"
## $ call : language lm(formula = pasa_ability ~ group_fct, data = df_eq_pasa_out)
## $ terms :Classes 'terms', 'formula' language pasa_ability ~ group_fct
## .. ..- attr(*, "variables")= language list(pasa_ability, group_fct)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "pasa_ability" "group_fct"
## .. .. .. ..$ : chr "group_fct"
## .. ..- attr(*, "term.labels")= chr "group_fct"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(pasa_ability, group_fct)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## .. .. ..- attr(*, "names")= chr [1:2] "pasa_ability" "group_fct"
## $ model :'data.frame': 107 obs. of 2 variables:
## ..$ pasa_ability: num [1:107] 4.25 2.75 3.5 3 3.75 2.25 3 3.5 2.5 2.25 ...
## ..$ group_fct : Factor w/ 2 levels "Group 1","Group 2": 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language pasa_ability ~ group_fct
## .. .. ..- attr(*, "variables")= language list(pasa_ability, group_fct)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "pasa_ability" "group_fct"
## .. .. .. .. ..$ : chr "group_fct"
## .. .. ..- attr(*, "term.labels")= chr "group_fct"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(pasa_ability, group_fct)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## .. .. .. ..- attr(*, "names")= chr [1:2] "pasa_ability" "group_fct"
## - attr(*, "class")= chr "lm"broom::tidy(model_AbilityIndex_out)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.08 0.0982 31.4 3.93e-55
## 2 group_fctGroup 2 -0.315 0.140 -2.25 2.62e- 2tab_model(model_AbilityIndex_out)| pasa ability | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 3.08 | 2.88 – 3.27 | <0.001 |
| group_fct [Group 2] | -0.31 | -0.59 – -0.04 | 0.026 |
| Observations | 107 | ||
| R2 / R2 adjusted | 0.046 / 0.037 | ||
# Ability Index model results
#model_AbilityIndex_out_results <- apa_print(model_AbilityIndex_out)
#Calculate a cohen's D effect size (to compare to multiple R squared below).
library(lsr)
cohensD(pasa_ability ~ group_fct, data = df_eq_pasa)## [1] 0.5107648Stress Index
# Box Plot
pasa_box_stress <- boxplot(pasa_stress~Condition,data=df_eq_pasa, main="Stress Index Scores by Condition",
xlab="Condition", ylab="Stress index") #Boxplot for Ability Index by condition
pasa_box_stress## $stats
## [,1] [,2]
## [1,] -4.500 -2.70
## [2,] -1.400 -0.80
## [3,] 0.175 0.25
## [4,] 1.200 1.30
## [5,] 3.650 4.15
##
## $n
## [1] 56 53
##
## $conf
## [,1] [,2]
## [1,] -0.3739546 -0.2057624
## [2,] 0.7239546 0.7057624
##
## $out
## [1] 4.55
##
## $group
## [1] 2
##
## $names
## [1] "1" "2"Running the regression of Stress Index
model_StressIndex_AsFac <- lm(pasa_stress ~ group_fct, # outcome (Stress Index) ~ predictor (Condition)
data = df_eq_pasa)
summary(model_StressIndex_AsFac)##
## Call:
## lm(formula = pasa_stress ~ group_fct, data = df_eq_pasa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3991 -1.2575 0.0009 1.2509 4.2425
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1009 0.2202 -0.458 0.648
## group_fctGroup 2 0.4084 0.3158 1.294 0.199
##
## Residual standard error: 1.648 on 107 degrees of freedom
## Multiple R-squared: 0.0154, Adjusted R-squared: 0.006194
## F-statistic: 1.673 on 1 and 107 DF, p-value: 0.1986broom::tidy(model_StressIndex_AsFac)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.101 0.220 -0.458 0.648
## 2 group_fctGroup 2 0.408 0.316 1.29 0.199tab_model(model_StressIndex_AsFac)| pasa stress | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -0.10 | -0.54 – 0.34 | 0.648 |
| group_fct [Group 2] | 0.41 | -0.22 – 1.03 | 0.199 |
| Observations | 109 | ||
| R2 / R2 adjusted | 0.015 / 0.006 | ||
# Stress Index model results
#model_StressIndex_AsFac_results <- apa_print(model_StressIndex_AsFac)
#Calculate a cohen's D effect size (to compare to multiple R squared below).
library(lsr)
cohensD(pasa_ability ~ group_fct, data = df_eq_pasa)## [1] 0.5107648Write-Up of Stress Index
A simple linear regression was performed to predict self-perceived stress (threat index - ability index) based on condition. The mean self-perceived stress score for Condition 1 (M = -0.101, SE = 0.235) was not significantly different than Condition 2 (M = 0.308, SE = 0.209) (F(107) = 1.673, p = 0.199), with an R2 of 0.015. Removing outliers from the data did not impact this effect (F(105) = 2.75, p = 0.100), with an R2 of 0.026.
Detecting Outliers
#Obtaining values of outliers
library(dplyr)
library(ggplot2)
is_outlier <- function(x) {
return(x < quantile(x, 0.25) - 1.5 * IQR(x) | x > quantile(x, 0.75) + 1.5 * IQR(x))
}
df_eq_pasa %>%
group_by(Condition) %>%
mutate(outlier = ifelse(is_outlier(pasa_stress), pasa_stress, as.numeric(NA))) %>%
ggplot(., aes(x = factor(Condition), y = pasa_stress)) +
geom_boxplot() +
geom_text(aes(label = outlier), na.rm = TRUE, hjust = -0.3)
# Removing Outliers
df_eq_pasa_stress_out1 <- df_eq_pasa %>%
filter(Condition == 1, pasa_ability > 0)
df_eq_pasa_stress_out2 <- df_eq_pasa %>%
filter(Condition == 2, pasa_ability < 4.375)
df_eq_pasa_stress_out <- rbind(df_eq_pasa_stress_out1, df_eq_pasa_stress_out2)Re-running the regression of Stress Index w/ outliers removed
model_StressIndex_AsFac_out <- lm(pasa_stress ~ group_fct, # outcome (Stress Index) ~ predictor (Condition)
data = df_eq_pasa_stress_out)
summary(model_StressIndex_AsFac_out)##
## Call:
## lm(formula = pasa_stress ~ group_fct, data = df_eq_pasa_stress_out)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.3991 -1.2334 0.0009 1.2166 4.1324
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1009 0.2159 -0.467 0.641
## group_fctGroup 2 0.5185 0.3127 1.658 0.100
##
## Residual standard error: 1.615 on 105 degrees of freedom
## Multiple R-squared: 0.02552, Adjusted R-squared: 0.01624
## F-statistic: 2.75 on 1 and 105 DF, p-value: 0.1002broom::tidy(model_StressIndex_AsFac_out)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.101 0.216 -0.467 0.641
## 2 group_fctGroup 2 0.519 0.313 1.66 0.100tab_model(model_StressIndex_AsFac_out)| pasa stress | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -0.10 | -0.53 – 0.33 | 0.641 |
| group_fct [Group 2] | 0.52 | -0.10 – 1.14 | 0.100 |
| Observations | 107 | ||
| R2 / R2 adjusted | 0.026 / 0.016 | ||
# Stress Index model results
#model_StressIndex_AsFac_results_out <- apa_print(model_StressIndex_AsFac_out)
#Plot model
library(ggplot2)
# get the group means:
with(df_eq_pasa_stress_out, by(pasa_stress, Condition, mean))## Condition: 1
## [1] -0.1008929
## ------------------------------------------------------------
## Condition: 2
## [1] 0.4176471# saving group means and group labels into a new data frame:
stress.means <- c(-0.1009, 0.4176)
stress.groups <- factor(c('Condition 1','Condition 2'))
means.data <- data.frame(stress.means, stress.groups)
# plotting means
ggplot(means.data, aes(y=stress.means, x=stress.groups)) + geom_bar(stat='identity')
# add error bars
standard.error <- function(x){sd(x)/sqrt(length(x))}
with(df_eq_pasa_stress_out, by(pasa_stress, Condition, standard.error))## Condition: 1
## [1] 0.2347577
## ------------------------------------------------------------
## Condition: 2
## [1] 0.2022318means.data$ses <- c(.23, .20)
# add the errorbar geom where the minimum is the mean - 1 standard error and the maximum is the mean + 1 standard error
#ggplot(means.data, aes(y=pasa_stress, x=Condition)) + geom_bar(stat='identity') + geom_errorbar(ymin=stress.means-means.data$ses, ymax=stress.means+means.data$ses) #doesn't work yet
#Calculate a cohen's D effect size (to compare to multiple R squared below).
library(lsr)
cohensD(pasa_ability ~ group_fct, data = df_eq_pasa_stress_out)## [1] 0.6221246#Checking the assumptions of the Stress Index Model
#Assumption of normality met
model_StressIndex_AsFac %>%
augment() %>%
ggplot(aes(x= .resid)) +
geom_density(fill = "purple") +
stat_function(linetype = 2, fun = dnorm, # add normal curve
args = list(mean = mean(augment(model_StressIndex_AsFac)$.resid), # define mean and sd
sd = sd(augment(model_StressIndex_AsFac)$.resid)))+
theme_minimal() #The residuals look fairly normally distributed
model_StressIndex_AsFac %>%
augment() %>%
ggplot() +
stat_qq(aes(sample = .std.resid)) + # Create the qq plot data. "Stat' function does something with a statistical function.
geom_abline() + # diagonal line based on intercept and slope, will default to zero (int) and one(slope).
theme_minimal() 
#Residuals are fairly close to the diagonal here, suggesting that the errors are normally distributed.
#Assumption of homoscedasticity - not met. Not sure if this is relevant for an RCT?
model_StressIndex_AsFac %>%
augment() %>%
ggplot(aes(x = .fitted, y = .resid)) +
geom_point() +
theme_minimal() +
geom_hline(yintercept = 0, color = "red")
# model comparison results
#mod_comp_results <- apa_print(list(model_ThreatIndex_results, model_ChallengeIndex_results, model_AbilityIndex_results, model_ControlIndex_results, model_PrimaryIndex_results, model_SecondIndex_results, model_StressIndex_AsFac_results),#add all models
# setting to 0
# otherwise it creates bootstrapped
# CI for delta R^2 - we will cover bootstrapping
# soon.
#boot_samples = 0)
#mod_comp_results$table %>%
#apa_table(caption = "Model Comparisons",
# note = "* p < 0.05; ** p < 0.01; *** p < 0.001",
# escape = TRUE)