Studio di riferimento

L’esperimento si compone di due fattori sperimentali (match \(\times\) identity) ciascuno con due livelli; i soggetti coinvolti sono 15, ad ognuno dei quali sono stati somministrati 200 trial, 50 per condizione.

Per il nuovo studio si prende in considerazione la sola condizione di matching e si vuole confrontare un campione di controllo con soggetti di un campione target.

Problema: i soggetti target sono pochi (molta memoria).

Obiettivo della power è individuare quanti soggetti servono (per gruppo) tenendo conto che il numero di trial non può essere aumentato ulteriormente.

Gamma distribution

Per modellare i RT utilizziamo la Gamma.

La distribuzione Gamma è definita con \[\begin{eqnarray}\label{eq:dgamma1} f(y|\alpha,\lambda) = \frac{1}{\Gamma(\alpha)} \lambda^{\alpha} y^{\alpha-1} e^{-\lambda y} \end{eqnarray}\] con \(y>0\), in cui \(\alpha>0\) indica il parametro di shape e \(\lambda>0\) il parametro di rate. La distribuzione ha valore atteso \(E(y)=\alpha/\lambda\) e varianza \(\mbox{Var}(Y)=\alpha/ \lambda^2\).

Il modello generativo dei dati è formalizzato nel seguente modo: \[\begin{eqnarray}\label{eq:generativemodel} \begin{array}{cr} RT \sim \text{Gamma}(\alpha, \lambda) & \text{Modello dei dati osservati}\\ \alpha & \text{Parametro di shape}\\ \lambda = \alpha / exp(\mu) & \text{Parametro di rate}\\ \mu = {\mathbf{X}\beta} & \text{Componente lineare} \end{array} \end{eqnarray}\]

Ovvero, i RT si distribuiscono secondo una funzione Gamma di parametri \(\alpha\) (shape) e \(\lambda\) (rate); \(\lambda\) è dato dal rapporto \(\alpha/exp(\mu)\) in cui \(\mu\) è il valore atteso della componente lineare del modello, da cui si ricava: \(E(y)=\alpha/\lambda = \alpha/(\alpha/exp(\mu)) = exp(\mu)\).

In sintesi: dobbiamo individuare i vettori di parametri \(\mathbf{\beta}\) che producono le differenze di RT attese nelle varie condizioni sperimentali.

Stime dei parametri sui dati osservati

Importante Utilizziamo i dati disponibili per valutare: 1) se la Gamma funziona sufficientemente bene; 2) per avere un’idea della variabilità individuale.

> library( brms )
> fitRT <- brm( RT ~ identity + (1|ID),
+       data = subset( X, match == "matching"  & ACCURATEZZA == 1), family = "gamma", seed = SEED )
 Family: gamma 
  Links: mu = log; shape = identity 
Formula: RT ~ identity + (1 | ID) 
   Data: subset(X, match == "matching" & ACCURATEZZA == 1) (Number of observations: 1350) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~ID (Number of levels: 15) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.12      0.03     0.07     0.18 1.00      611     1148

Population-Level Effects: 
           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     -0.16      0.03    -0.22    -0.09 1.00      639      998
identityTU    -0.22      0.01    -0.24    -0.19 1.00     3968     3012

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape    16.50      0.63    15.26    17.73 1.00     3159     2670

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
> pp_check(fitRT)

Effetti ipotizzati

Calcolo dei parametri

Per calcolare i parametri, utilizziamo la parte lineare del modello:

\[\begin{eqnarray} \beta = \mu \mathbf{X}^{-1} \end{eqnarray}\]

con la seguente funzione:

> expected_pars_RT <- function(MU){
+   
+   DESIGN <- expand.grid( 
+     I = c("sconosciuto", "tu"), 
+     G = c("control","target") )
+   DESIGN$MU <- MU
+   MM <- model.matrix( MU ~ I*G, data = DESIGN)
+   
+   ETA <- matrix( log(MU/1000) )
+   DESIGN$BE <- solve( MM, ETA )
+   return(DESIGN)
+   
+ }
> 
> expected_pars_RT( c( 850, 750, 850, 650 ))$BE
                  [,1]
(Intercept) -0.1625189
Itu         -0.1251631
Gtarget      0.0000000
Itu:Gtarget -0.1431008

Per simulare RT coerenti con lo scenario 2 il modello sarà:

\[\begin{eqnarray}\label{eq:exp1simdata1} \begin{array}{l} \mu = (-0.16+\delta_i) -0.13 I + 0 G -0.14 IG\\ \delta \sim Normal( 0, \sigma_{\delta}) \\ \lambda = 16 / exp(\mu) \\ RT \sim Gamma( 16, \lambda) \end{array} \end{eqnarray}\] in cui \(I \in \{0,1\}\) indica le due condizioni sconosciuto/tu, \(G \in \{0,1\}\) i due gruppi – controllo e target – e \(\delta\) rappresenta la variabilità dei soggetti intorno all’intercetta generale; per semplicità la assumiamo costante anche rispetto ai trials.

Come deviazione standard della variabilità individuale \(\sigma_{\delta}\) aumentiamo leggermente il valore stimato empiricamente (0.12) a 0.15. Dal momento che anche la variabilità individuale incide sulla potenza, meno differenze tra soggetti implica più potenza, assumendo un maggior grado di differenza ci poniamo in una posizione più prudente.

Prova di simulazione

Proviamo a generare un campione di 1000 soggetti (500 per gruppo) con 50 trial per condizione.

Stimiamo i parametri:

 Family: gamma 
  Links: mu = log; shape = identity 
Formula: RT ~ I * G + (1 | ID) 
   Data: Z2 (Number of observations: 100000) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~ID (Number of levels: 1000) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.14      0.00     0.14     0.15 1.00      735     1526

Population-Level Effects: 
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept      -0.15      0.01    -0.16    -0.14 1.01      395      697
Itu            -0.12      0.00    -0.13    -0.12 1.00    14533     8278
Gtarget        -0.00      0.01    -0.02     0.01 1.02      360      912
Itu:Gtarget    -0.14      0.00    -0.14    -0.13 1.00    15819     7539

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape    16.06      0.07    15.93    16.20 1.00    18398     7175

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Visualizziamo i RT medi (attesi) per condizione:

Power analysis

Scegliamo 5 numerosità campionarie (per gruppo): \(n = 10, 15, 20, 25, 30\) e, per ciascuna di esse, ripetiamo per 600 volte il seguente algoritmo:

  1. Simuliamo un set di dati con \(n\) soggetti per gruppo e 50 trials per condizione: ogni soggetto ha 100 osservazioni, 50 nella condizione sconosciuto e 50 nella condizione tu.
  2. Stimiamo i parametri del modello utilizzando le seguenti prior
                 prior     class        coef
1   student_t(3, 0, 1) Intercept            
2 student_t(3, 0, .04)         b         Itu
3 student_t(3, 0, .04)         b     Gtarget
4 student_t(3, 0, .04)         b Itu:Gtarget
5   student_t(3, 0, 2)        sd            
6        gamma(17, 17)     shape

Si tratta di prior scettiche, sulla base delle quali, assumiamo che, con una probabilità di 90%, i coefficienti di regressione cadano nell’intervallo \([-0.09, 0.09]\). Questo intervallo, che possiamo definire null interval, definisce in pratica i valori per i quali gli effetti si considerano sostanzialmente nulli. Per l’intercetta invece la prior è più ampia, si assume, con la stessa probabilità, che i valori cadano nell’intervallo \([-2.35, 2.35]\); ricordiamo che non ci interessa fare la power sull’intercetta.

  1. Per ciascun parametro calcoliamo gli intervalli di credibilità al 90% e 80% (Credible Interval, CI) ovvero gli estremi che contengono rispettivamente queste due porzioni di distribuzioni a posteriori. La scelta dei valori 80% e 90% è giustificata dal limite di soggetti che si potranno reclutare e dal fatto che gli effetti attesi sono piuttosto piccoli.

Per calcolare la potenza contiamo quante volte l’intervallo è completamente fuori dal range \([-0.09, 0.09]\). In pratica, la potenza, per ogni parametro e numerosità, viene stimata come rapporto tra il numero di volte in cui l’IC è fuori dall’intervallo nullo ed il totale di iterazioni per numerosità (600). Dal momento che gli IC 80% sono più stretti di quelli 90%, con i primi ci attendiamo un livello di potenza maggiore.

Visualizziamo le medie ottenute nei 12000 campioni simulati:

E quindi le curve di potenza per i \(\beta\) del modello: