forked from Louisenyholm/Assignment-2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathA2_P3_LanguageASD_power_instructions.Rmd
176 lines (123 loc) · 8.33 KB
/
A2_P3_LanguageASD_power_instructions.Rmd
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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
---
title: "Assignment 1 - Language Development in ASD - Power and simulations"
author: "[YOUR NAME]"
date: "[DATE]"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Welcome to the third exciting part of the Language Development in ASD exercise
In this part of the assignment, we try to figure out how a new study should be planned (i.e. how many participants?) in order to have enough power to replicate the findings (ensuring our sample size is adequate, our alpha at 0.05 and our beta at 0.8):
1- if we trust the estimates of the current study. Report the power analysis and comment on what you can (or cannot) use its estimates for.
2- if we are skeptical of the current study. Report the power analysis and comment on what you can (or cannot) use its estimates for.
3- if we only have access to 30 participants. Identify the power for each relevant effect and discuss whether it's worth to run the study and why
The list above is also what you should discuss in your code-less report.
## Learning objectives
- Learn how to calculate statistical power
- Critically appraise how to apply frequentist statistical power
### Exercise 1
How much power does your study have (if your model estimates are quite right)?
- Load your dataset (both training and testing), fit your favorite model, assess power for your effects of interest (probably your interactions).
- Report the power analysis and comment on what you can (or cannot) use its estimates for.
- Test how many participants you would have to have to replicate the findings (assuming the findings are correct)
N.B. Remember that main effects are tricky once you have interactions in the model (same for 2-way interactions w 3-way interactions in the model). If you want to test the power of main effects, run a model excluding the interactions.
N.B. Check this paper: https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504
You will be using:
- powerSim() to calculate power
- powerCurve() to estimate the needed number of participants
- extend() to simulate more participants
```{r}
setwd("C:/Users/louis/OneDrive - Aarhus universitet/AU Onedrive - RIGTIG/- 3. Semester/Experimental Methods III/Classes/Assignment 2/Assignment-2")
# Installing and loading packages
# install.packages("simr", dependencies = T)
# library(pacman)
# p_load(tidyverse, lme4)
# library(simr)
library(githubinstall)
githubinstall("simr", lib = .libPaths())
print(pacman::p_path() == .libPaths())
pacman::p_load(tidyverse, lmerTest, simr)
# Loading data
df_train <- read_csv("train.csv")
df_test <- read_csv("test.csv")
# Making the test set an extension of the training set regarding IDs
df_test$Child.ID <- 66 + df_test$Child.ID
# Merging training and test data
df <- merge(df_train, df_test, all = T)
# Removing NAs
df <- df[-which(is.na(df$CHI_MLU)),]
# We assume a normal distribution, because we have seen that human perform normally-distributed in similar intelligent tasks.
# Fitting a model
model0 <- lmer(CHI_MLU ~ Visit + Diagnosis + verbalIQ1 + (1 + Visit|Child.ID), df, REML = F, control = lmerControl(optimizer = "nloptwrap", calc.derivs = F))
model <- lmer(CHI_MLU ~ Visit * Diagnosis * verbalIQ1 + (1 + Visit|Child.ID), df, REML = F, control = lmerControl(optimizer = "nloptwrap", calc.derivs = F))
# Checking the main effect estimates
fixef(model0)["Visit"] # main effect of 0.23
fixef(model0)["DiagnosisTD"] # main effect of 0.17
fixef(model0)["verbalIQ1"] # main effect of 0.08
# Checking the interaction effect estimate
fixef(model)["Visit:DiagnosisTD:verbalIQ1"] # interaction effect of -0.02
# Running the power analysis
set.seed(123) # Setting a fixed "random"
# On the main effects (model without interactions)
powerSim(model0, test = fixed("Visit"), nsim = 100) # 100% power for detecting effect of visit (nsim = 100)
powerSim(model0, test = fixed("Diagnosis"), nsim = 100) # 52% power for detecting effect of diagnosis (nsim = 100)
powerSim(model0, test = fixed("verbalIQ1"), nsim = 100) # 100% power for detecting effect of verbal IQ (nsim = 100)
# On the interaction effect
powerSim(model, test = fixed("Visit:Diagnosis:verbalIQ1"), nsim=1000) # 98.30% power for detecting the three-way interaction effect (nsim = 1000)
# Report the power analysis and comment on what you can (or cannot) use its estimates for.
# The interaction of Visit * Diagnosis * VerbalIQ1
# Calculating power based on 1000 simulations gives a beta/power value of 0.97. This means, that an interaction effect of this size (-0.02) is easily detected with a sample like this. It might even be an overpowered study.
# COMMENT ON MAIN EFFECTS
# Correspondingly, the main effects of visit and verbal IQ have a great power (100%). The effect of diagnosis is however just slightly above chance (0.52) for detecting an effect, when it is present.
# Test how many participants you would have to have to replicate the findings (assuming the findings are correct)
# We assume that we do not need to extend our simulation in most of following assessments, since the powers of the most effects are already above 0.80. It is rather a question about not overpowering the study. We do not need this many participants to detect these effects - neither for two of the main effects nor interaction effect.
# However, for the main effect of diagnosis, an extension - to simulate more participants - might be needed.
model_extended <- extend(model0, along = "Child.ID", n = 200)
# Main effects
pc_Visit <- powerCurve(model0, test = fixed("Visit"), along = "Child.ID", nsim = 100)
print(pc_Visit)
plot(pc_Visit)
# In order to replicate the findings of a main effect of Visit, one would need around 11 participants (resulting in a beta-value/power of .90)
pc_Diagnosis <- powerCurve(model_extended, test = fixed("Diagnosis"), along = "Child.ID", nsim = 100)
print(pc_Diagnosis)
plot(pc_Diagnosis)
# In order to replicate the findings of a main effect of Diagnosis, one would need around 156 participants (resulting in a beta-value/power of .80)
pc_VerbalIQ <- powerCurve(model0, test = fixed("verbalIQ1"), along = "Child.ID", nsim = 100)
print(pc_VerbalIQ)
plot(pc_VerbalIQ)
# In order to replicate the findings of a main effect of verbal IQ, one would need around 11 participants (resulting in a beta-value/power of .88)
pc_inter <- powerCurve(model, test = fixed("Visit:Diagnosis:verbalIQ1"), along = "Child.ID", nsim = 100)
print(pc_inter)
plot(pc_inter)
# In order to replicate the findings of a three-way interaction effect, one would need around 34 participants (resulting in a beta-value/power of .84)
```
### Exercise 2
How would you perform a more conservative power analysis?
- Identify and justify a minimum effect size for each of your relevant effects
- take the model from exercise 1 and replace the effects with the minimum effect size that you'd accept.
- assess the power curve by Child.ID, identifying an ideal number of participants to estimate each effect
- if your power estimates do not reach an acceptable threshold simulate additional participants and repeat the previous analysis
- Report the power analysis and comment on what you can (or cannot) use its estimates for.
```{r}
# - Identify and justify a minimum effect size for each of your relevant effects
# JUSTIFY A MINIMUM EFFECT SIZE
# For visit
ggplot(df, aes(Visit, CHI_MLU, colour = Diagnosis))+
geom_point()+
geom_smooth(method = lm)
# - Take the model from exercise 1 and replace the effects with the minimum effect size that you'd accept.
# Specifying a minimum effect size for the three-way interaction effect
fixef(model)["Visit:DiagnosisTD:verbalIQ1"] <- 0.015
# - Assess the power curve by Child.ID, identifying an ideal number of participants to estimate each effect
pc_inter2 <- powerCurve(model, test = fixed("Visit:Diagnosis:verbalIQ1"), along = "Child.ID", nsim = 100)
print(pc_inter2)
plot(pc_inter2)
# In order to find a minimum effect size of xxx of the three-way interaction effect, one would need xx participants (resulting in a beta-value/power of .xx)
# - if your power estimates do not reach an acceptable threshold simulate additional participants and repeat the previous analysis
# - Report the power analysis and comment on what you can (or cannot) use its estimates for.
```
### Exercise 3
Assume you have only the resources to collect 30 kids (15 with ASD and 15 TDs). Identify the power for each relevant effect and discuss whether it's worth to run the study and why
```{r}
```