Introduzione ai modelli di equazioni struttural

27. Introduzione ai modelli di equazioni struttural#

Un modello di equazioni strutturali (SEM) è una classe generale di tecniche multivariate che combina l’analisi fattoriale confermativa e l’analisi dei percorsi. La modellizzazione delle equazioni strutturali comprende due componenti principali: il modello di misura e il modello strutturale. Il “modello di misurazione” descrive le relazioni tra variabili latenti e variabili osservate (ovvero, può esere rappresentato come un modello fattoriale), mentre “il modello strutturale” è un modello di analisi dei percorsi che specifica le relazioni tra le variabili latenti.

Il processo di analisi SEM è composto dai seguenti passaggi e decisioni:

Costruire un diagramma dei percorsi che mostri il modello di misura e strutturale di interesse.
Identificare il livello di misura per ogni item e verificare le ipotesi distributive.
Assicurarsi che la funzione di adattamento scelta sia basata sui tipi di misura (ad esempio, massima verosimiglianza per misure continue, minimi quadrati ponderato per misure ordinali).
Adattare il modello utilizzando la funzione di adattamento appropriata e valutare l’adattamento del modello utilizzando un insieme di indici.
Una volta stabilito un modello plausibile, interpretare i vari parametri a livello di elemento (ad esempio, saturazioni fattoriali, errori standard, valori R-quadrati, termini di errore, ecc.)

suppressPackageStartupMessages({
    library("lavaan")
    library("lavaanExtra")
    library("lavaanPlot")
    library("psych")
    library("dplyr") 
    library("tidyr")
    library("knitr")
    library("mvnormalTest")
    library("semPlot")
    library("DiagrammeRsvg")
    library("rsvg")
})
options(repr.plot.width=6, repr.plot.height=6)
set.seed(42)

In questo esempio useremo il classico dataset di Holzinger e Swineford (1939) che consiste nei punteggi dei test di abilità mentale di bambini di settima e ottava elementare di due scuole diverse (Pasteur e Grant-White). Nel dataset originale, ci sono punteggi per 26 test. Tuttavia, un sottoinsieme più piccolo con 9 variabili è più ampiamente utilizzato nella letteratura (ad esempio, Joreskog del 1969).

x1: Visual perception
x2: Cubes
x3: Lozenges
x4: Paragraph comprehension
x5: Sentence completion
x6: Word meaning
x7: Speeded addition
x8: Speeded counting of dots
x9: Speeded discrimination straight and curved capitals

Leggiamo i dati in R.

data(HolzingerSwineford1939)

Esaminiamo i dati.

glimpse(HolzingerSwineford1939) 

Rows: 301
Columns: 15
$ id     <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, …
$ sex    <int> 1, 2, 2, 1, 2, 2, 1, 2, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 2, 2, 1, …
$ ageyr  <int> 13, 13, 13, 13, 12, 14, 12, 12, 13, 12, 12, 12, 12, 12, 12, 12,…
$ agemo  <int> 1, 7, 1, 2, 2, 1, 1, 2, 0, 5, 2, 11, 7, 8, 6, 1, 11, 5, 8, 3, 1…
$ school <fct> Pasteur, Pasteur, Pasteur, Pasteur, Pasteur, Pasteur, Pasteur, …
$ grade  <int> 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, …
$ x1     <dbl> 3.333333, 5.333333, 4.500000, 5.333333, 4.833333, 5.333333, 2.8…
$ x2     <dbl> 7.75, 5.25, 5.25, 7.75, 4.75, 5.00, 6.00, 6.25, 5.75, 5.25, 5.7…
$ x3     <dbl> 0.375, 2.125, 1.875, 3.000, 0.875, 2.250, 1.000, 1.875, 1.500, …
$ x4     <dbl> 2.333333, 1.666667, 1.000000, 2.666667, 2.666667, 1.000000, 3.3…
$ x5     <dbl> 5.75, 3.00, 1.75, 4.50, 4.00, 3.00, 6.00, 4.25, 5.75, 5.00, 3.5…
$ x6     <dbl> 1.2857143, 1.2857143, 0.4285714, 2.4285714, 2.5714286, 0.857142…
$ x7     <dbl> 3.391304, 3.782609, 3.260870, 3.000000, 3.695652, 4.347826, 4.6…
$ x8     <dbl> 5.75, 6.25, 3.90, 5.30, 6.30, 6.65, 6.20, 5.15, 4.65, 4.55, 5.7…
$ x9     <dbl> 6.361111, 7.916667, 4.416667, 4.861111, 5.916667, 7.500000, 4.8…

27.1. Valutazione delle assunzioni#

Un’analisi SEM inizia con una valutazione sulla distribuzione delle variabili endogene.

temp <- HolzingerSwineford1939 |>
    select(-c(id, sex, school, ageyr, agemo, grade))

describe(temp) |>
    print()

   vars   n mean   sd median trimmed  mad  min   max range  skew kurtosis   se
x1    1 301 4.94 1.17   5.00    4.96 1.24 0.67  8.50  7.83 -0.25     0.31 0.07
x2    2 301 6.09 1.18   6.00    6.02 1.11 2.25  9.25  7.00  0.47     0.33 0.07
x3    3 301 2.25 1.13   2.12    2.20 1.30 0.25  4.50  4.25  0.38    -0.91 0.07
x4    4 301 3.06 1.16   3.00    3.02 0.99 0.00  6.33  6.33  0.27     0.08 0.07
x5    5 301 4.34 1.29   4.50    4.40 1.48 1.00  7.00  6.00 -0.35    -0.55 0.07
x6    6 301 2.19 1.10   2.00    2.09 1.06 0.14  6.14  6.00  0.86     0.82 0.06
x7    7 301 4.19 1.09   4.09    4.16 1.10 1.30  7.43  6.13  0.25    -0.31 0.06
x8    8 301 5.53 1.01   5.50    5.49 0.96 3.05 10.00  6.95  0.53     1.17 0.06
x9    9 301 5.37 1.01   5.42    5.37 0.99 2.78  9.25  6.47  0.20     0.29 0.06

I valori di asimmetria e kurtosi sono accettabili.

Valutiamo l’ipotesi di normalità multivariata. Utilizziamo il pacchetto mvnormalTest per un test sulla normalità univariata (test W di Shapiro-Wilk) e multivariata (test di asimmetria e curtosi multivariata di Mardia).

Iniziamo con la normalità univariata.

mvnout <- mardia(temp)
## Shapiro-Wilk Univariate normality test
mvnout$uv.shapiro

   W      p-value UV.Normality
x1 0.9928 0.1582  Yes         
x2 0.9697 0       No          
x3 0.9523 0       No          
x4 0.9827 0.0011  No          
x5 0.9769 1e-04   No          
x6 0.9538 0       No          
x7 0.9908 0.056   Yes         
x8 0.9807 4e-04   No          
x9 0.9942 0.307   Yes         

mvnout$mv.test |>
    print()

          Test Statistic p-value Result
   Skewness  344.9053       0     NO
   Kurtosis    3.2344  0.0012     NO
MV Normality      <NA>    <NA>     NO

I risultati dei test univariati e multivariati indicano che le misure non provengono da distribuzioni univariate o multivariate normalmente distribuite (i risultati ‘No’ nella tabella). Affronteremo questi problemi nella successiva fase di adattamento del modello ai dati.

27.2. Specificazione del modello#

Definiamo il modello utilizzando la sintassi lavaan.

model <- "

    # [-----Latent variables (measurement model)-----]

    visual =~ x1 + x2 + x3
    textual =~ x4 + x5 + x6
    speed =~ x7 + x8 + x9

    # [-----------Mediations (named paths)-----------]

    speed ~ visual_speed*visual
    textual ~ visual_textual*visual
    visual ~ ageyr_visual*ageyr + grade_visual*grade

    # [---------Regressions (Direct effects)---------]

    speed ~ ageyr + grade
    textual ~ ageyr + grade

    # [------------------Covariances-----------------]

    speed ~~ textual
    ageyr ~~ grade

    # [--------Mediations (indirect effects)---------]

    ageyr_visual_speed := ageyr_visual * visual_speed
    ageyr_visual_textual := ageyr_visual * visual_textual
    grade_visual_speed := grade_visual * visual_speed
    grade_visual_textual := grade_visual * visual_textual

"

Dato che abbiamo riscontrato violazioni dell’ipotesi di normalità multivariata richiesta per i modelli SEM, useremo lo stimatore “MLM” come funzione di adattamento. Questo stimatore utilizza una procedura di massima verosimiglianza e fornisce errori standard robusti e una statistica di test scalata di Satorra-Bentler per affrontare i problemi di violazione della normalità multivariata.

È importante notare che i problemi con i dati non normali conducono ad una sottostima degli errori standard, il che può portare a respingere troppo spesso l’ipotesi nulla che un parametro sia zero e ad un’inflazione della statistica chi-quadrato del modello, portando a respingere troppo spesso il modello.

fit <- sem(model, data=HolzingerSwineford1939, std.lv = TRUE, estimator = "MLM")

27.3. Bontà di adattamento#

In genere, i ricercatori esaminano le statistiche di adattamento del modello prima di procedere all’interpretazione delle stime dei parametri. L’ipotesi nulla in un’analisi SEM è che la matrice di covarianza riprodotta dal modello sia statisticamente la stessa della matrice di covarianza di input. Contrariamente al solito test di ipotesi, desideriamo che le due matrici siano statisticamente le stesse.

Iniziamo a valutare l’adattamento del modello con un test chi-quadrato ottenuto dall’output di lavaan come segue:

fitMeasures(fit, c("chisq.scaled", "df.scaled", "pvalue.scaled")) |>
    print()

 chisq.scaled     df.scaled pvalue.scaled 
      110.247        36.000         0.000 

Il rapporto tra chi-quadrato e gradi di libertà è minore di 4, il che è accettabile, anche se indicativo di un fit non eccellente.

La Root Mean Square Error of Approximation (RMSEA) è una misura della discrepanza tra le matrici di correlazione basate sul modello e quelle osservate. Utilizza il chi-quadrato del modello nel suo calcolo ma apporta correzioni in base alla complessità del modello (correzione per la parsimonia) e ha una distribuzione campionaria nota in modo da poter calcolare gli intervalli di confidenza. Abbiamo ottenuto i valori RMSEA scalati dall’output di lavaan come segue:

fitMeasures(fit, c("rmsea.scaled", "rmsea.ci.lower.scaled", "rmsea.ci.upper.scaled")) |>
    print()

         rmsea.scaled rmsea.ci.lower.scaled rmsea.ci.upper.scaled 
                0.083                 0.066                 0.100 

Sono state proposte varie linee guida per l’interpretazione dell’RMSEA:

RMSEA <= .05 come soglia per un buon adattamento;
RMSEA = .05 - .08 come adattamento ragionevole;
RMSEA >= .10 come adattamento scarso.

Sulla base della stima puntuale RMSEA ottenuta = .083 e dell’intervallo di confidenza al 90% [.066, .100], concludiamo che il modello ha un adattamento appena accettabile.

Per valutare l’adeguatezza del modello, utilizziamo due ulteriori misure di adattamento: l’Indice di Adattamento Comparativo (CFI) e il Residuo Quadratico Medio Radice Standardizzato (srmr). Il CFI è un indice di adattamento incrementale che confronta il modello considerato con un modello di base ristretto. L’srmr, invece, si basa sulle discrepanze tra le covarianze previste dal modello e le covarianze effettive.

Abbiamo ottenuto i valori scalati del CFI e dell’srmr dall’output di lavaan come segue:

fitMeasures(fit, c("cfi.scaled", "srmr")) |>
    print()

cfi.scaled       srmr 
     0.925      0.060 

Per l’interpretazione di queste misure sono state proposte diverse linee guida. In questo esempio, abbiamo utilizzato come valori soglia CFI >= .90 e srmr <= .08. In base a queste soglie, abbiamo concluso che il valore ottenuto di CFI.scaled = .925 e srmr = .060 fornivano ulteriori prove che il nostro modello si adattava ai dati in modo soddisfacente.

Sulla base di questo insieme di misure di adattamento, possiamo concludere che il modello specificato è plausibile.

27.4. Modello di misurazione#

Dato il modello accettabile, passiamo all’esame delle varie stime dei parametri. Concentrandoci prima sul modello di misura, abbiamo ottenuto le stime dall’output di lavaan come segue:

standardizedsolution(fit, type = "std.all", se = TRUE, zstat = TRUE, pvalue = TRUE, ci = TRUE) %>%
    filter(op == "=~") %>%
    select(LV = lhs, Item = rhs, Coefficient = est.std, ci.lower, ci.upper, SE = se, Z = z, "p-value" = pvalue) |>
    print()

       LV Item Coefficient ci.lower ci.upper    SE      Z p.value
visual   x1       0.778    0.644    0.911 0.068 11.402       0
visual   x2       0.425    0.314    0.536 0.057  7.497       0
visual   x3       0.577    0.471    0.683 0.054 10.693       0
textual   x4       0.856    0.809    0.902 0.024 35.915       0
textual   x5       0.858    0.817    0.899 0.021 41.296       0
textual   x6       0.830    0.784    0.876 0.023 35.601       0
 speed   x7       0.599    0.500    0.698 0.050 11.897       0
 speed   x8       0.751    0.665    0.837 0.044 17.081       0
 speed   x9       0.625    0.536    0.714 0.045 13.813       0

Questo output presenta i coefficienti standardizzati (saturazioni fattoriali) per gli item sulle variabili latenti (LV), gli intervalli di confidenza (ci.lower, ci.upper), gli errori standard (SE), i valori Z (test di Wald) e i valori p che testano l’ipotesi nulla che un coefficiente = 0. Le saturazioni fattoriali variavano da .42 a .86, indicando che l’entità delle relazioni tra gli item e i fattori era adeguata (sebbene non ci siano soglie rigide per i carichi accettabili). Si noti che gli errori standard sono robusti, il che significa che sono corretti per le influenze della non-normalità. Si noti anche che tutti i coefficienti sono statisticamente significativi, il che significa che l’ipotesi nulla che un coefficiente = 0 è respinta.

È anche utile esaminare i valori R2, ovvero le saturazioni fattoriali standardizzate al quadrato degli item. Nel framework SEM, qualsiasi variabile che ha una freccia rivolta verso di essa è una variabile endogena e avrà un valore R2 associato ad essa. I valori R2 mostrati di seguito per ogni item indicano la proporzione di varianza di quell’item spiegata dalla variabile latente corrispondente. Più alta è la percentuale di varianza di un item spiegata dal fattore, maggiore è il contributo dell’item per la misurazione del fattore. Abbiamo ottenuto i coefficienti R2 dall’output di lavaan come segue:

parameterEstimates(fit, standardized = TRUE, rsquare = TRUE) %>%
  filter(op == "r2") %>%
  select(Item = rhs, R2 = est) |>
  print()

      Item    R2
     x1 0.605
     x2 0.180
     x3 0.333
     x4 0.732
     x5 0.736
     x6 0.689
     x7 0.359
     x8 0.563
     x9 0.390
visual 0.091
textual 0.327
 speed 0.312

I valori R2 per gli item variano da 0.18 a 0.74 e suggeriscono che gli item hanno una relazione da piccola a sostanziale con una variabile latente. Si noti che non esiste un limite rigido per i valori R2 accettabili, ma valori oltre .50 sono desiderabili.

27.5. Modello strutturale#

Ora ci concentriamo sul modello strutturale. Otteniamo le stime dall’output di lavaan come segue:

standardizedsolution(fit, type = "std.all", se = TRUE, zstat = TRUE, pvalue = TRUE, ci = TRUE) %>%
  filter(op == "~") %>%
  select(LV = lhs, Item = rhs, Coefficient = est.std, ci.lower, ci.upper, SE = se, Z = z, "p-value" = pvalue) |>
  print()

       LV   Item Coefficient ci.lower ci.upper    SE      Z p.value
 speed visual       0.380    0.224    0.535 0.079  4.792   0.000
textual visual       0.366    0.216    0.516 0.077  4.776   0.000
visual  ageyr      -0.220   -0.377   -0.063 0.080 -2.744   0.006
visual  grade       0.348    0.199    0.498 0.076  4.558   0.000
 speed  ageyr       0.124   -0.018    0.266 0.073  1.711   0.087
 speed  grade       0.271    0.111    0.431 0.082  3.325   0.001
textual  ageyr      -0.385   -0.491   -0.280 0.054 -7.147   0.000
textual  grade       0.323    0.190    0.456 0.068  4.743   0.000

Questo output presenta coefficienti di regressione standardizzati che rappresentano le relazioni tra variabili esogene e le variabili endogene, intervalli di confidenza (ci.lower, ci.upper), errori standard (SE), valori Z (test di Wald), valori p che testano l’ipotesi nulla che un coefficiente = 0. Un coefficiente di regressione rappresenta la forza della relazione tra una variabile esogena e una variabile endogena e il segno rappresenta la direzione della relazione.

27.6. Diagramma di percorso#

Per generare il diagramma di percorso del modello adattato usiamo la funzione nice_lavaanPlot del pacchetto lavaanExtra.

nice_lavaanPlot(fit, stand = TRUE, stars= FALSE)

27.7. Conclusioni#

È possibile interpretare i risultati ottenuti come segue. Il modello complessivo sembra plausibile sulla base di vari indici di adattamento (anche con un numero ridotto di variabili misurate). Le variabili latenti ipotizzate (visual, speed, textual) risultano adeguatamente misurate dagli indicatori associati considerando le saturazioni fattoriali (da piccole ad elevate) e i valori R2. Le relazioni strutturali indicano che il grado scolastico influenza in maniera positiva i punteggi delle variabili latenti textual e speed, mentre l’età ha un effetto negativo su textual, ma un effetto positivo su speed. La variabile latente visual ha un effetto positivo sia su textual che su speed. Nel complesso, il modello spiega circa un terzo della varianza di textual e speed, ma solo il 9% della varianza della abilità cognitiva visual.