18.6 Path analysis con lavaan | Appunti di Costruzione e validazione di strumenti di misura dell’efficacia dell’intervento psicologico in neuropsicologia

18.6 Path analysis con `lavaan`

Usiamo lavaan per svolgere l’analisi statistica descritta nell’esempio precedente.

Esercizio 18.1 La matrice di correlazioni di partenza è:

lower <- "
  1
  .70 1
  .65 .50 1
"

Converto tali dati in una matrice simmetrica.

dat.cov <- getCov(lower, names = c("y", "x1", "x2"))
dat.cov
#>       y  x1   x2
#> y  1.00 0.7 0.65
#> x1 0.70 1.0 0.50
#> x2 0.65 0.5 1.00

Specifico il modello con la sintassi di lavaan.

mr_model <- "y ~ x1 + x2"

Adatto il modello ai dati.

fit <- sem(mr_model, sample.cov = dat.cov, sample.nobs = 100)

Esamino i risultati.

summary(fit, standardized = TRUE)
#> lavaan 0.6.15 ended normally after 1 iteration
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                         3
#> 
#>   Number of observations                           100
#> 
#> Model Test User Model:
#>                                                       
#>   Test statistic                                 0.000
#>   Degrees of freedom                                 0
#> 
#> Parameter Estimates:
#> 
#>   Standard errors                             Standard
#>   Information                                 Expected
#>   Information saturated (h1) model          Structured
#> 
#> Regressions:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   y ~                                                                   
#>     x1                0.500    0.072    6.934    0.000    0.500    0.500
#>     x2                0.400    0.072    5.547    0.000    0.400    0.400
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>    .y                 0.386    0.055    7.071    0.000    0.386    0.390

Con la funzione sem() del pacchetto lavaan abbiamo dunque replicato i risultati ottenuti in precedenza. Il valore \(0.386\) rappresenta la quota di varianza della \(y\) non spiegata dalle variabili esogene.

Esercizio 18.2 Ripeto ora la stessa procedura simulando 100 osservazioni su tre variabili, in maniera tale da avere a disposizione i dati grezzi. Imponendo un effetto causale diretto delle variabili \(x_1\) e \(x_2\) sulla \(y\) e una correlazione \(> 0\) tra le variabili \(x_1\) e \(x_2\), otteniamo i dati seguenti.

set.seed(3)
n <- 100
x1 <- rnorm(n, 100, 9)
x2 <- x1 + rnorm(n, 0, 10)
cor(x1, x2)
#> [1] 0.5346271
y <- 10 + 3 * x1 + 1.5 * x2 + rnorm(n, 0, 15)
dd <- data.frame(y, x1, x2)
print(cor(dd), 3)
#>        y    x1    x2
#> y  1.000 0.831 0.786
#> x1 0.831 1.000 0.535
#> x2 0.786 0.535 1.000

Per fare un esempio, immaginiamo che queste correlazioni siano state ricavate da qualche fonte e dunque le leggiamo nella memoria di lavoro di \(\mathsf{R}\) nel modo seguente.

lower <- "
 1
 .831 1
 .786 .535 1
"
dat.cov <- getCov(lower, names = c("y", "x1", "x2"))
dat.cov
#>        y    x1    x2
#> y  1.000 0.831 0.786
#> x1 0.831 1.000 0.535
#> x2 0.786 0.535 1.000

Data una matrice di correlazioni e data la specificazione delle relazioni tra le variabili, la funzione sem() contenuta nel pacchetto lavaan consente di stimare i coefficienti di percorso.

mr.model <- "y ~ x1 + x2"
fit <- sem(mr.model, sample.cov = dat.cov, sample.nobs = 100)

Esaminiamo i risultati con la funzione summary().

summary(fit, standardized = TRUE)
#> lavaan 0.6.15 ended normally after 1 iteration
#> 
#>   Estimator                                         ML
#>   Optimization method                           NLMINB
#>   Number of model parameters                         3
#> 
#>   Number of observations                           100
#> 
#> Model Test User Model:
#>                                                       
#>   Test statistic                                 0.000
#>   Degrees of freedom                                 0
#> 
#> Parameter Estimates:
#> 
#>   Standard errors                             Standard
#>   Information                                 Expected
#>   Information saturated (h1) model          Structured
#> 
#> Regressions:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>   y ~                                                                   
#>     x1                0.575    0.045   12.710    0.000    0.575    0.575
#>     x2                0.478    0.045   10.571    0.000    0.478    0.478
#> 
#> Variances:
#>                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
#>    .y                 0.145    0.020    7.071    0.000    0.145    0.146

Ora calcoliamo i coefficienti del modello di regressione multipla usando i dati grezzi (standardizzati). Si noti che i risultati ottenuti da lavaan sono identici a quelli prodotti dal modello di regressione multipla – nel caso presente, infatti, il modello statistico esaminato da sem() non era altro che il modello di regressione multipla.

summary(lm(scale(y) ~ scale(x1) + scale(x2)))
#> 
#> Call:
#> lm(formula = scale(y) ~ scale(x1) + scale(x2))
#> 
#> Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.87200 -0.24623  0.01628  0.25714  0.90382 
#> 
#> Coefficients:
#>               Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) -7.472e-16  3.867e-02     0.0        1    
#> scale(x1)    5.748e-01  4.599e-02    12.5   <2e-16 ***
#> scale(x2)    4.785e-01  4.599e-02    10.4   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3867 on 97 degrees of freedom
#> Multiple R-squared:  0.8535, Adjusted R-squared:  0.8505 
#> F-statistic: 282.6 on 2 and 97 DF,  p-value: < 2.2e-16

La quota di varianza non spiegata della variabile endogena è:

1 - .145
#> [1] 0.855