Setting iniziale e visualizzazione effetto
Stabilisco i parametri del modello
intercept = 0
multitaskingEffect = -.28
dysfluencyEffect = .39
interactionEffect = .20
sigma = 1
Calcolo i livelli medi nelle varie condizioni per visualizzarli
# visualize effect
mean_MT0_DF0 = intercept
mean_MT0_DF1 = intercept + dysfluencyEffect
mean_MT1_DF0 = intercept + multitaskingEffect
mean_MT1_DF1 = intercept + dysfluencyEffect + multitaskingEffect + interactionEffect
Visualizzo con grafico ggplot
dataplot = data.frame(estimate = c(mean_MT0_DF0, mean_MT0_DF1, mean_MT1_DF0, mean_MT1_DF1),
multitasking = c("MT0","MT0","MT1","MT1"),
dysfluency = c("DF0","DF1","DF0","DF1"))
dataplot$dysfluency = factor(dataplot$dysfluency,levels=c("DF1","DF0"))
ggplot(data=dataplot,aes(y=estimate,x=multitasking,group=dysfluency,color=dysfluency,shape=dysfluency)) +
geom_hline(yintercept=0, linetype=2, size=0.8)+
geom_point(size=7) +
geom_line(size=1.5) +
theme(text=element_text(size=27)) +
scale_y_continuous(breaks=seq(-2,2,.1)) +
ylab("z-score (1 = SD)")
Simulo un singolo (piccolo) dataset
N = 16
n = N/4
id = 1:N
MT = c( rep(0,n), rep(0,n), rep(1,n), rep(1,n) )
DY = c( rep(0,n), rep(1,n), rep(0,n), rep(1,n) )
id; MT; DY
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## [1] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
## [1] 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1
# determino lo score "atteso" per ogni soggetto
score_atteso = intercept +
MT * multitaskingEffect +
DY * dysfluencyEffect +
MT * DY * interactionEffect
cbind(id, MT, DY, score_atteso)
## id MT DY score_atteso
## [1,] 1 0 0 0.00
## [2,] 2 0 0 0.00
## [3,] 3 0 0 0.00
## [4,] 4 0 0 0.00
## [5,] 5 0 1 0.39
## [6,] 6 0 1 0.39
## [7,] 7 0 1 0.39
## [8,] 8 0 1 0.39
## [9,] 9 1 0 -0.28
## [10,] 10 1 0 -0.28
## [11,] 11 1 0 -0.28
## [12,] 12 1 0 -0.28
## [13,] 13 1 1 0.31
## [14,] 14 1 1 0.31
## [15,] 15 1 1 0.31
## [16,] 16 1 1 0.31
# determino e aggiungo la varianza residua
residual = rnorm(N,0,sigma)
score = score_atteso + residual
# metto tutto nel dataset
df = data.frame(id=id, multitasking=MT, dysfluency=DY, score=score)
df
## id multitasking dysfluency score
## 1 1 0 0 -0.14180297
## 2 2 0 0 0.85175640
## 3 3 0 0 -0.33599124
## 4 4 0 0 -0.41089111
## 5 5 0 1 0.68914635
## 6 6 0 1 -1.21302280
## 7 7 0 1 1.07210925
## 8 8 0 1 0.68658147
## 9 9 1 0 -0.25401097
## 10 10 1 0 -0.67662425
## 11 11 1 0 0.69632704
## 12 12 1 0 0.02401117
## 13 13 1 1 0.06071607
## 14 14 1 1 1.20459748
## 15 15 1 1 -0.18771538
## 16 16 1 1 0.54085832
anche se inutile qui, fitto e guardo il modello
fit1 = lm(score ~ multitasking * dysfluency, data=df)
summary(fit1)
##
## Call:
## lm(formula = score ~ multitasking * dysfluency, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.52173 -0.35834 -0.02799 0.47256 0.86099
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.009232 0.363356 -0.025 0.980
## multitasking -0.043342 0.513863 -0.084 0.934
## dysfluency 0.317936 0.513863 0.619 0.548
## multitasking:dysfluency 0.139253 0.726713 0.192 0.851
##
## Residual standard error: 0.7267 on 12 degrees of freedom
## Multiple R-squared: 0.0895, Adjusted R-squared: -0.1381
## F-statistic: 0.3932 on 3 and 12 DF, p-value: 0.7602
Simulo 2000 volte
ora faccio la power analysis (con N tot = 320) e vedo anche le stime
ottenute
# decido il numero di iterazioni
niter = 2000
# stabilisco un N un po' decente
N = 320
n = round(N/4)
# inizializzo (per praticità) il vettore dei pvalue, dei CI e degli effetti stimati
pvalue = rep(NA,niter)
ci_width = rep(NA,niter)
estimated_interactionEffect = rep(NA,niter)
# lancio il ciclo for
for(i in 1:niter){
# simulo i punteggi
id = 1:N
MT = c(rep(0,n),rep(0,n),rep(1,n),rep(1,n))
DY = c(rep(0,n),rep(1,n),rep(0,n),rep(1,n))
score_atteso = intercept +
MT * multitaskingEffect +
DY * dysfluencyEffect +
MT * DY * interactionEffect
residual = rnorm(N,0,sigma)
score = score_atteso + residual
# metto tutto nel dataset
df = data.frame(id=id, multitasking=MT, dysfluency=DY, score=score)
# calcolo il p-value dall'anova
fit1 = lm(score ~ multitasking * dysfluency, data=df)
fit0 = lm(score ~ multitasking + dysfluency, data=df)
pvalue[i] = anova(fit1,fit0)$"Pr(>F)"[2]
# estraggo il coefficiente di interazione stimato e lo salvo nel vettore
estimated_interactionEffect[i] = fit1$coefficients["multitasking:dysfluency"]
# estraggo i CI e li salvo nel vettore
ci = confint(fit1)
lower = ci["multitasking:dysfluency",1]
upper = ci["multitasking:dysfluency",2]
ci_width[i] = (upper-lower)
}
Calcolo il power
mean(pvalue < .05, na.rm=T)
## [1] 0.1575
Determino il minimo effetto significativo
minsign = min( estimated_interactionEffect [pvalue < .05 & estimated_interactionEffect > 0], na.rm=T )
minsign
## [1] 0.3878323
Guardo la distribuzione degli effetti stimati e dove si colloca
quello minimo significativo (linea tratteggiata rossa) e quello “vero” (linea blu)
hist(estimated_interactionEffect, breaks=30)
abline(v=minsign, col="red",lwd=4,lty=2)
abline(v=interactionEffect, col="blue",lwd=4,lty=1)
Do anche un’occhiata alla distribuzione dei CI
hist(ci_width/2, breaks=20)
abline(v=mean(ci_width/2),lwd=8,lty=1)
