Studio preliminare

Author

Progetto

Vogliamo studiare la relazione tra numero di spriz, tipo di software utilizzato e stress.

Il modello che vogliamo valutare è

\[ stress = \beta_0 + \beta_1 \times spriz + \beta_2 \times software + \beta_3 \times spriz \times software + \epsilon \] Usiamo le seguenti variabili: - stress, scala a 100 punti - spriz, numero di spriz e software utilizzato, R o excel

Design analysis

Ipotizziamo che lo stress di distribuisca normalmente con media 70 e dev.st. 10, quindi graficamente:

Definizione dei parametri

Vogliamo simulare i punteggi standardizzati di stress; per semplicità consideriamo il numero di spriz come variabile centrata sulla media ed il tipo di software codificato con -.5, .5 al posto della classica codifica 0-1.

I parametri del modello per la simulazione saranno i seguenti:

\(\beta_0 = 0\)
\(\beta_1 = -0.3\)
\(\beta_2 = 0.01\)
\(\beta_3 = 0.1\)

Data simulation

Simuliamo 1000 soggetti.

Iniziamo con la variabile spriz:

spriz <- rpois( N, lambda )
spriz_s <- scale( spriz )
plot( table( spriz ), type = "h" )

Poi generiamo la variabile software:

software <- sample( c(-.5,.5), N, replace = TRUE )

e per finire la variabile stress, sulla base del modello:

stress <- b0 + b1*spriz_s + b2*software + b3*spriz_s*software + rnorm( N )

plot( density( stress ) )

Aggreghiamo le variabili in un data frame e le rappresentiamo con grafici univariati e bivariati:

myData <- data.frame(
  stress, spriz = spriz_s, software
)

library( GGally )
ggpairs( myData )

Stima dei parametri dai dati simulati

Per stimare i parametri, utilizziamo un modello lineare bayesiano

library( brms )
fit <- brm( stress ~ spriz*software, data = myData, cores = 4,
            prior = myPrior )

con le seguenti prior:

                 prior     class           coef
1 student_t( 3, 0, 1 )         b               
2                              b       software
3                              b          spriz
4                              b spriz:software
5 student_t(3, 0, 2.5) Intercept               
6 student_t(3, 0, 2.5)     sigma

Visualizziamo le posterior ed i traceplot del modello:

Le stime dei parametri sono le seguenti:

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: stress ~ spriz * software 
   Data: myData (Number of observations: 1000) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept         -0.02      0.03    -0.09     0.04 1.00     4803     3205
spriz             -0.30      0.03    -0.36    -0.24 1.00     5312     3161
software          -0.04      0.06    -0.16     0.09 1.00     4833     2940
spriz:software     0.11      0.07    -0.02     0.24 1.00     4617     3194

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.01      0.02     0.97     1.06 1.00     5775     3337

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).

Ciclo di simulazione

Per determinare quanti soggetti servono per stimare bene questi parametri esploriamo i seguenti livelli di \(n\): 20, 50, 100.

PARTAB <- NULL
for (n in n_vec) {
  for ( b in 1:B) {
    # genero i predittori
    spriz <- rpois( n, lambda )
    spriz_s <- scale( spriz )
    software <- sample( c(-.5,.5), n, replace = TRUE )
    
    # genero la dipendente
    stress <- b0 + b1*spriz_s + b2*software + b3*spriz_s*software + rnorm( n )
    
    # model fit 
    Z <- data.frame(
      stress, spriz = spriz_s, 
      software
    )
    
    if (b == 1) {
      fit <- brm( stress ~ spriz*software, data = Z, cores = 4,
                prior = myPrior )
    } else {
      fit <- update( fit, newdata = Z )
    }
    
    # recupero le stime dei parametri e gli intervalli di credibilità
    PAR <- data.frame( fixef( fit, probs = round(PROBS,3) ) )
    PAR$par <- row.names(PAR)
    PAR$n <- n
    
    PARTAB <- rbind( PARTAB, PAR )
    
  }
}

Vediamo se le stime dei parametri sono coerenti con l’attesa:

Ed infine, le curve di potenza, ottenute contando quante volte l’intervallo di credibilità 89% non contiene lo zero: