-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbehavior_GLMM_analysis.R
144 lines (116 loc) · 4.67 KB
/
behavior_GLMM_analysis.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
cat("\f")
rm(list=ls())
# Init
graphics.off()
library("MASS")
library(afex)
library(phia)
library(doBy)
library(effsize)
library(lmerTest);
library(dplyr);
library("sjPlot")
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
#Import data
file_name = "data/srm_post_norm_for_glmm.csv" # SRM relationship data of Q2-Q1 normalisation
#read data
data = read.table(file_name, header=TRUE,sep=',')
head(data)
#defining outcomes
relationship <-as.numeric(data$relationship)
#defining predictors (categorical)
participant_condition <- as.factor(data$participant_condition)
participant_condition <- relevel(participant_condition, "U")
other_condition <- as.factor(data$other_condition)
other_condition <- relevel(other_condition, "U")
participant_nb <- as.factor(data$participant_nb)
interacting_partner <- as.factor(data$interacting_partner)
question_content <- as.factor(data$question_content)
dyad <- as.factor(data$dyad)
sex <- as.factor(data$sex)
#create new dataframe
all_data=data.frame( relationship, participant_condition , interacting_partner, dyad, other_condition, participant_nb, question_content, sex )
## -------------------------------------- ##
## -------------------------------------- ##
## GLMM Analysis ------------------------ ##
## -------------------------------------- ##
## -------------------------------------- ##
# Choose the question you want to anlayse
#question = "other_seeing_me_again" # Other seeing me again
question = "good_conversation" # Conversation quality data
#question = "seeing_again" # Romantic attraction data
#question = "smile" # Other smile data
clean_data <- all_data[all_data$question_content == question,]
head(clean_data)
#---- Linear Model, without random factors
nul<- lm( relationship ~ 1 , data= clean_data)
v0 <- lm( relationship ~ participant_condition , data= clean_data)
v1 <- lm( relationship ~ other_condition , data= clean_data)
v2 <- lm( relationship ~ other_condition + participant_condition , data= clean_data)
v3 <- lm( relationship ~ other_condition * participant_condition , data= clean_data)
aov.out = anova(nul,v0)
aov.out
aov.out = anova(nul,v1)
aov.out
aov.out = anova(nul,v2)
aov.out
aov.out = anova(nul,v3)
aov.out
aov.out = anova(v2,v3)
aov.out
## ---- With both participant_nb and interacting_partner as random factors
nul<- lmer( relationship ~ 1 + (1 | participant_nb) + (1 | interacting_partner), data= clean_data, REML = FALSE )
v0 <- lmer( relationship ~ participant_condition + (1 | participant_nb) + (1 | interacting_partner), data= clean_data, REML = FALSE )
v1 <- lmer( relationship ~ other_condition + (1 | participant_nb) + (1 | interacting_partner), data= clean_data, REML = FALSE )
v2 <- lmer( relationship ~ other_condition + participant_condition + (1 | participant_nb) + (1 | interacting_partner), data= clean_data, REML = FALSE )
v3 <- lmer( relationship ~ other_condition * participant_condition + (1 | participant_nb) + (1 | interacting_partner), data= clean_data, REML = FALSE )
aov.out = anova(nul,v0)
aov.out
aov.out = anova(nul,v1)
aov.out
aov.out = anova(nul,v2)
aov.out
aov.out = anova(nul,v3)
aov.out
#Interaction
aov.out = anova(v2,v3)
aov.out
#summary
summary(v3)
#plot
plot_model( v3, type = "pred", terms = c("participant_condition", "other_condition") )
## ---- Corrected post-hocs
library(PairedData)
conds = c("U", "S")
bonferoni_alpha = 0.05/6 # for bonferoni correction
for (pc in conds){
for (oc in conds){
#choose one measure
reduced = clean_data[clean_data$participant_condition==pc,]
reduced = reduced[reduced$other_condition==oc,]
#Choose another one
for (pc2 in conds){
for (oc2 in conds){
reduced2 = clean_data[clean_data$participant_condition==pc2,]
reduced2 = reduced2[reduced2$other_condition==oc2,]
#choose one measure
x = reduced$relationship
y = reduced2$relationship
res = t.test(x, y, paired = TRUE)
if(is.nan(res$p.value)){
next
}
if (res$p.value < bonferoni_alpha)
{
d = round(effectsize::cohens_d(x, y, paired = TRUE), digit=5)
print(paste("pc : ", pc, " oc : ", oc , " vs ", "pc : ", pc2, " oc : ", oc2 ))
print(paste("p = ", round(res$p.value, digit=5) , "--- t = "
, round(res$statistic[[1]], digit=5)
, "d : ", d$Cohens_d ," df = ", res$parameter)
)
print("")
}
}
}
}
}