25. Dati non gaussiani e categoriali#
Negli esempi precedenti di questa dispensa è stato utilizzato lo stimatore di massima verosimiglianza (ML). Molti dei modelli CFA e SEM riportati nella letteratura di ricerca applicata utilizzano infatti stime di ML. Tuttavia, lo stimatore ML è appropriata solo per dati multivariati normali a livello di scala a intervalli (cioè, quando la distribuzione congiunta delle variabili continue è distribuita normalmente). Quando i dati continui si discostano marcatamente dalla normalità (cioè, mostrano una forte asimmetria o curtosi), o quando alcuni degli indicatori non sono a livello di scala a intervalli (cioè, dati binari, politomici o ordinali), allora è preferibile usare uno stimatore diverso da quello di ML.
25.1. Dati non gaussiani#
La ricerca ha dimostrato che la stima di ML è robusta nel caso di piccole deviazioni nella normalità. Tuttavia, quando la non normalità è più pronunciata, è necessario utilizzare uno stimatore diverso dalla ML per ottenere risultati statistici affidabili (vale a dire, statistiche accurate sulla bontà dell’adattamento ed errori standard accurati delle stime dei parametri). La stima di ML è particolarmente sensibile ad un eccesso di curtosi.
Le conseguenze dell’uso della ML in condizioni di grave non normalità includono
valori eccessivi della statistica \(\chi^2\) del modello;
la sottostima degli indici di bontà dell’adattamento mediante indici quali TLI e CFI;
la sottostima degli errori standard dei parametri.
Questi effetti deleteri sono esacerbati dalla diminuzione della dimensione del campione.
Allo scopo di limitare tali conseguenze indesiderate, nelle condizioni di marcata violazione dell’assunzione di normalità multivariata, vengono usati stimatori diversi dalla ML. I due stimatori più comunemente usati per dati continui non normali sono
ML robusto,
minimi quadrati ponderati (WLS).
L’uso di WLS non è, in generale, raccomandato, a meno che le dimensioni del campione non siano molto grandi. Al contrario, la ricerca ha dimostrato che il metodo ML robusto fornisce uno stimatore adeguato rispetto a diversi livelli di non normalità.
Esaminiamo qui un esempio discusso da Brown [Bro15] nelle sue tabelle 9.5 – 9.7.
d <- readRDS(here::here("data", "brown_table_9_5_data.RDS"))
head(d)
| x1 | x2 | x3 | x4 | x5 | |
|---|---|---|---|---|---|
| <int> | <int> | <int> | <int> | <int> | |
| 1 | 0 | 0 | 0 | 0 | 0 |
| 2 | 0 | 0 | 0 | 0 | 0 |
| 3 | 0 | 0 | 0 | 0 | 0 |
| 4 | 4 | 2 | 2 | 1 | 1 |
| 5 | 1 | 0 | 1 | 6 | 0 |
| 6 | 0 | 0 | 0 | 0 | 0 |
Le statistiche descrittive di questo campione di dati mostrano valori eccessivi di asimmetria e di curtosi.
psych::describe(d)
| vars | n | mean | sd | median | trimmed | mad | min | max | range | skew | kurtosis | se | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| <int> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | <dbl> | |
| x1 | 1 | 870 | 1.4701149 | 2.172832 | 0 | 1.0086207 | 0 | 0 | 8 | 8 | 1.506406 | 1.252591 | 0.07366591 |
| x2 | 2 | 870 | 0.8229885 | 1.601474 | 0 | 0.4152299 | 0 | 0 | 8 | 8 | 2.398394 | 5.670143 | 0.05429505 |
| x3 | 3 | 870 | 1.2655172 | 2.070024 | 0 | 0.7772989 | 0 | 0 | 8 | 8 | 1.797942 | 2.343203 | 0.07018040 |
| x4 | 4 | 870 | 1.0264368 | 1.928047 | 0 | 0.5359195 | 0 | 0 | 8 | 8 | 2.157445 | 3.977564 | 0.06536693 |
| x5 | 5 | 870 | 0.6068966 | 1.519175 | 0 | 0.1839080 | 0 | 0 | 8 | 8 | 3.103965 | 9.373781 | 0.05150485 |
Definiamo un modello ad un fattore e, seguendo @brown2015confirmatory, aggiungiamo una correlazione residua tra gli indicatori X1 e X3:
model <- '
f1 =~ x1 + x2 + x3 + x4 + x5
x1 ~~ x3
'
Procediamo alla stima dei parametri utilizzando uno stimatore di ML robusto. La sintassi lavaan è la seguente:
fit <- cfa(model, data = d, mimic = "MPLUS", estimator = "MLM")
Per esaminare la soluzione ottenuta ci focalizziamo sulla statistica \(\chi^2\) – si consideri la soluzione robusta fornita nell’output.
out = summary(fit)
print(out)
lavaan 0.6.15 ended normally after 28 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 16
Number of observations 870
Model Test User Model:
Standard Scaled
Test Statistic 25.913 10.356
Degrees of freedom 4 4
P-value (Chi-square) 0.000 0.035
Scaling correction factor 2.502
Satorra-Bentler correction (Mplus variant)
Parameter Estimates:
Standard errors Robust.sem
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
f1 =~
x1 1.000
x2 0.703 0.062 11.338 0.000
x3 1.068 0.044 24.304 0.000
x4 0.918 0.063 14.638 0.000
x5 0.748 0.055 13.582 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
.x1 ~~
.x3 0.655 0.143 4.579 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
.x1 1.470 0.074 19.968 0.000
.x2 0.823 0.054 15.166 0.000
.x3 1.266 0.070 18.043 0.000
.x4 1.026 0.065 15.712 0.000
.x5 0.607 0.051 11.790 0.000
f1 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.x1 2.040 0.228 8.952 0.000
.x2 1.241 0.124 10.019 0.000
.x3 1.227 0.169 7.255 0.000
.x4 1.458 0.177 8.233 0.000
.x5 0.807 0.100 8.063 0.000
f1 2.675 0.289 9.273 0.000
Per fare un confronto, adattiamo lo stesso modello ai dati usando lo stimatore di ML.
fit2 <- cfa(model, data = d)
Notiamo come il valore della statistica \(\chi^2\) ora ottenuto sia molto maggiore di quello trovato in precedenza.
out = summary(fit2)
print(out)
lavaan 0.6.15 ended normally after 28 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 11
Number of observations 870
Model Test User Model:
Test statistic 25.913
Degrees of freedom 4
P-value (Chi-square) 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
f1 =~
x1 1.000
x2 0.703 0.035 20.133 0.000
x3 1.068 0.034 31.730 0.000
x4 0.918 0.042 21.775 0.000
x5 0.748 0.033 22.416 0.000
Covariances:
Estimate Std.Err z-value P(>|z|)
.x1 ~~
.x3 0.655 0.091 7.213 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.x1 2.040 0.128 15.897 0.000
.x2 1.241 0.070 17.671 0.000
.x3 1.227 0.095 12.942 0.000
.x4 1.458 0.090 16.135 0.000
.x5 0.807 0.053 15.119 0.000
f1 2.675 0.220 12.154 0.000
25.2. Dati categoriali#
Quando almeno un indicatore è categoriale (cioè binario, politomico o ordinale), il metodo ML ordinario non dovrebbe essere utilizzato per stimare i modelli CFA. Vi sono molte potenziali conseguenze del trattamento delle variabili categoriali come continue in un’analisi CFA, incluso il fatto che può tale scelta può
produrre stime attenuate delle relazioni tra gli indicatori, specialmente quando ci sono effetti pavimento o soffitto;
portare ad individuare “pseudo-fattori” che emergono come artefatti del metodo statistico;
produrre distorsioni negli indici di bontà dell’adattamento e nelle stime degli errori standard;
produrre stime errate dei parametri.
Esistono vari stimatori che possono essere utilizzati con indicatori categoriali; ad esempio, gli stimatori dei minimi quadrati ponderati (WLS), dei minimi quadrati ponderati robusti (WLSMV) e dei minimi quadrati non ponderati (ULS).
25.2.1. Un esempio concreto#
Nell’esempio discusso da @brown2015confirmatory, i ricercatori desiderano verificare un modello uni-fattoriale di dipendenza da alcol in un campione di 750 pazienti ambulatoriali. Gli indicatori di alcolismo sono item binari che riflettono la presenza/assenza di sei criteri diagnostici per l’alcolismo (0 = criterio non soddisfatto, 1 = criterio soddisfatto). I dati sono i seguenti:
d1 <- readRDS(here::here("data", "brown_table_9_9_data.RDS"))
head(d1)
| y1 | y2 | y3 | y4 | y5 | y6 | |
|---|---|---|---|---|---|---|
| <int> | <int> | <int> | <int> | <int> | <int> | |
| 1 | 1 | 1 | 1 | 1 | 1 | 1 |
| 2 | 1 | 1 | 1 | 1 | 1 | 1 |
| 3 | 1 | 1 | 1 | 1 | 1 | 0 |
| 4 | 1 | 1 | 1 | 1 | 1 | 1 |
| 5 | 0 | 0 | 0 | 0 | 0 | 0 |
| 6 | 1 | 1 | 0 | 1 | 1 | 1 |
Il modello viene specificato nel modo seguente:
model1 <- '
etoh =~ y1 + y2 + y3 + y4 + y5 + y6
'
Adattiamo il modello specificando che i dati sono a livello di scala ordinale (stimatore WLSMVS).
fit1 <- cfa(
model1,
data = d1,
ordered = names(d1),
estimator = "WLSMVS",
mimic = "mplus"
)
Esaminiamo la soluzione ottenuta:
out = summary(fit1, fit.measures = TRUE)
print(out)
lavaan 0.6.15 ended normally after 16 iterations
Estimator DWLS
Optimization method NLMINB
Number of model parameters 12
Number of observations 750
Model Test User Model:
Standard Scaled
Test Statistic 5.651 9.540
Degrees of freedom 9 9
P-value (Chi-square) 0.774 0.389
Scaling correction factor 0.592
mean and variance adjusted correction (WLSMV)
Model Test Baseline Model:
Test statistic 1155.845 694.433
Degrees of freedom 15 9
P-value 0.000 0.000
Scaling correction factor 1.664
User Model versus Baseline Model:
Comparative Fit Index (CFI) 1.000 0.999
Tucker-Lewis Index (TLI) 1.005 0.999
Root Mean Square Error of Approximation:
RMSEA 0.000 0.009
90 Percent confidence interval - lower 0.000 0.000
90 Percent confidence interval - upper 0.028 0.051
P-value H_0: RMSEA <= 0.050 0.999 0.944
P-value H_0: RMSEA >= 0.080 0.000 0.000
Standardized Root Mean Square Residual:
SRMR 0.031 0.031
Parameter Estimates:
Standard errors Robust.sem
Information Expected
Information saturated (h1) model Unstructured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
etoh =~
y1 1.000
y2 0.822 0.072 11.392 0.000
y3 0.653 0.092 7.097 0.000
y4 1.031 0.075 13.703 0.000
y5 1.002 0.072 13.861 0.000
y6 0.759 0.076 10.011 0.000
Intercepts:
Estimate Std.Err z-value P(>|z|)
.y1 0.000
.y2 0.000
.y3 0.000
.y4 0.000
.y5 0.000
.y6 0.000
etoh 0.000
Thresholds:
Estimate Std.Err z-value P(>|z|)
y1|t1 -0.759 0.051 -14.890 0.000
y2|t1 -0.398 0.047 -8.437 0.000
y3|t1 -1.244 0.061 -20.278 0.000
y4|t1 -0.795 0.051 -15.436 0.000
y5|t1 -0.384 0.047 -8.148 0.000
y6|t1 -0.818 0.052 -15.775 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.399
.y2 0.594
.y3 0.744
.y4 0.361
.y5 0.397
.y6 0.653
etoh 0.601 0.063 9.596 0.000
Scales y*:
Estimate Std.Err z-value P(>|z|)
y1 1.000
y2 1.000
y3 1.000
y4 1.000
y5 1.000
y6 1.000
Confrontiamo la soluzione ottenuta con lo stimatore WLSMVS con quella ottenuta mediante lo stimatore ML.
fit2 <- cfa(
model1,
data = d1
)
out = summary(fit2, fit.measures = TRUE)
print(out)
lavaan 0.6.15 ended normally after 35 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 12
Number of observations 750
Model Test User Model:
Test statistic 14.182
Degrees of freedom 9
P-value (Chi-square) 0.116
Model Test Baseline Model:
Test statistic 614.305
Degrees of freedom 15
P-value 0.000
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.991
Tucker-Lewis Index (TLI) 0.986
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -2087.600
Loglikelihood unrestricted model (H1) -2080.508
Akaike (AIC) 4199.199
Bayesian (BIC) 4254.640
Sample-size adjusted Bayesian (SABIC) 4216.535
Root Mean Square Error of Approximation:
RMSEA 0.028
90 Percent confidence interval - lower 0.000
90 Percent confidence interval - upper 0.054
P-value H_0: RMSEA <= 0.050 0.914
P-value H_0: RMSEA >= 0.080 0.000
Standardized Root Mean Square Residual:
SRMR 0.021
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Latent Variables:
Estimate Std.Err z-value P(>|z|)
etoh =~
y1 1.000
y2 0.934 0.093 10.057 0.000
y3 0.390 0.055 7.038 0.000
y4 1.008 0.087 11.541 0.000
y5 1.158 0.101 11.468 0.000
y6 0.700 0.077 9.142 0.000
Variances:
Estimate Std.Err z-value P(>|z|)
.y1 0.109 0.007 14.692 0.000
.y2 0.169 0.010 16.781 0.000
.y3 0.085 0.005 18.483 0.000
.y4 0.102 0.007 14.285 0.000
.y5 0.140 0.010 14.506 0.000
.y6 0.132 0.008 17.514 0.000
etoh 0.065 0.009 7.664 0.000
Si noti che la soluzione ottenuta mediante lo stimatore WLSMVS produce indici di bontà di adattamento migliori e errori standard dei parametri più piccoli.