SEM class
Pablo Salazar
14/12/2023
Modelos estructurales
Un modelo de ecuaciones estructurales (o cualquier modelo matemático) nace a través de una hipótesis. Es la visión teórica de un sistema apoyada en el razonamiento de conocimientos previos. Sin embargo, para comprobar la validez de un modelo no es suficiente con su significancía.
Para comprobar la validez de un modelo es importante tener en cuenta:
¿El modelo representa adecuamente mis datos?
En los modelos estructurales “clásicos” (lavaan), la calidad del modelo en su conjunto se comprueba a través de:
Model FIt Test Statistic Degree of freedom P-value (Chi-square) = ^ 0.05
La hipotesis nula de este análisis establece que la covarianza establecida entre variables en nuestro modelo es indistiguible de la covarianza real encontrada en nuestros datos.
Varianza del modelo = varianza de los datos
¿Todos los componentes importantes de mi modelo representan adecuadamente mis datos?
Cada vez que una variable de respuesta es explicada por una o mas variables explicativas se crea un componente del modelo. Por ejemplo: el ejemplo anterior tendria 3 componentes.
Productividad ~ Fotosintesis Fotosintesis ~ Cierre estomatico Cierre estomatico ~ Sequia
Cada componente cuenta como un análisis de regresion independiente. El summary de Lavaan proporciona un R^2 para cada regresion segun el coeficiente de correlación parcial de las variables y la interrelacion de las variables del modelo.
Todas las relaciones del modelo se analizan simultaneamente para determinar la validez del modelo. Por tanto la validez de un componente puede verse seriamente afectada por otro componente.
Por ejemplo, si la ejecución de este modelo nos idican que la fotosíntesis no se puede explicar por el cierre estomatico:
Es posible que esto se deba a que el modelo asume que no existe una relación entre la Sequia y la fotosíntesis, y esto interefiere con la relación entre el cierre estomático y la fotosíntesis.
Es importante recordar que el número de variables no define el número de componentes del modelo. Por ejemplo, el siguiente modelo solo tiene 1 componente y se escribiría como:
Productividad ~ sequia + cierre_estomatico + Fotosintesis
Por tanto, al analizar la validez de un modelo, es posible que un componente o alguna relación dentro de un componente no represente adecuadamente los datos (no sea significativo).
¿Existen componentes importantes que no son parte de mi modelo?
En Lavaan, un modelo exitoso puede omitir relaciones importantes y seguir siendo un modelo exitoso. Sin embargo, un modelo fallido puede volverse exitoso si incluyes relaciones importantes que antes no tenias en cuenta.
Pueden revisar estas relaciones usando la función modindices, aunque siempre es importante que exista una base teorica para establecer estas nuevas relaciones.
¿Puedes identificar cuántos componentes tiene el modelo anterior?
¿Puedes identificar cuáles son?
Piecewise model
Los modelos Piecewise (piezas sabias) y los modelos clásicos (en Lavaan) tienen muchas diferencias entre si, como veras en la tabla siguiente. Pero quiza lo mas importante que debes saber es que cada componente (pieza) del modelo se analiza de manera independiente y tiene un efecto limitado sobre el siguiente.
| Lavaan | Piecewise |
|---|---|
|
|
|
|
|
|
|
|
|
|
¿Y si hacemos un modelooooo?
snowman_anatomy
Primero vamos a crear una base de datos a partir de variables creadas de manera aleatoria.
dat <- data.frame(Productividad = runif(50),
fotosintesis = runif(50),
cierre_estomatico = runif(50),
Sequia = runif(50))
summary(dat)## Productividad fotosintesis cierre_estomatico Sequia
## Min. :0.006809 Min. :0.002448 Min. :0.02212 Min. :0.0138
## 1st Qu.:0.189668 1st Qu.:0.245594 1st Qu.:0.25444 1st Qu.:0.2681
## Median :0.456925 Median :0.519435 Median :0.50412 Median :0.4640
## Mean :0.489459 Mean :0.514315 Mean :0.48481 Mean :0.4689
## 3rd Qu.:0.760275 3rd Qu.:0.821213 3rd Qu.:0.70232 3rd Qu.:0.6660
## Max. :0.998982 Max. :0.990174 Max. :0.97398 Max. :0.9527Ahora podemos intenar crear nuestro primer modelo. Piecewise usa la funcion “psem” para construir el modelo, el cual se construye a partir de la modelación de cada componente. Piecewise acepta varios algoritmos de modelación. En este caso usaremos la función “lm” (linear model).
model <- psem(
lm(Productividad ~ fotosintesis, dat), # Componente 1
lm(fotosintesis ~ cierre_estomatico, dat), # Componente 2
lm(cierre_estomatico ~ Sequia, dat) # Componente 3
)
summary(model, .progressBar = F)##
## Structural Equation Model of model
##
## Call:
## Productividad ~ fotosintesis
## fotosintesis ~ cierre_estomatico
## cierre_estomatico ~ Sequia
##
## AIC BIC
## 26.146 43.354
##
## ---
## Tests of directed separation:
##
## Independ.Claim Test.Type DF Crit.Value P.Value
## fotosintesis ~ Sequia + ... coef 47 0.8184 0.4172
## Productividad ~ Sequia + ... coef 47 0.7399 0.4630
## Productividad ~ cierre_estomatico + ... coef 46 -1.7423 0.0881
##
## Global goodness-of-fit:
##
## Fisher's C = 8.146 with P-value = 0.228 and on 6 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value
## Productividad fotosintesis -0.2418 0.1493 48 -1.6194 0.1119
## fotosintesis cierre_estomatico 0.1350 0.1567 48 0.8618 0.3931
## cierre_estomatico Sequia -0.1926 0.1458 48 -1.3208 0.1928
## Std.Estimate
## -0.2276
## 0.1234
## -0.1873
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method R.squared
## Productividad none 0.05
## fotosintesis none 0.02
## cierre_estomatico none 0.04El summary de Piecewise se compone de 4 grandes partes. Primero encontramos la sección “Call:”
Aquí aparece la estructura de nuestro modelo, el akaike information criterion (AIC), y el bayesian information criterion (BIC). Esto es muy util cuando queremos comparar modelos exitosos y tenemos que elegir cual explica mejor nuestros datos
La sección “Tests of directed separation:”
Un modelo Piecewise no puede medir la calidad del modelo en su conjunto porque en la práctica solo esta superponiendo una conjunto de modelos (componentes) uno sobre otro. Para determinar si el modelo es exitoso o no, Piecewise analiza las relaciones de las variables que no has tenido en cuenta. Si esas relaciones son significativas, significa que nuestros datos pueden ser explicados con relaciones adicionales y nuestro modelo debe ser corregido. Si esas relaciones no son significativas, significa que nuestro nuestros datos solo pueden modelarse con la relaciones que hemos establecido.
A partir de la significancia de la relaciones que NO has tenido en cuenta, Piecewise construye el estimador global Fisher’s C
\[C = -2\sum_{i=1}^{k}ln(p_i)\]
Donde k es el número de relaciones que no has tenido en cuenta y p es su significancia. Según Bill Shippley, el estimador C sigue una distribición chi cuadrado con 2k grados de libertad. Por tanto la calidad del modelo se mide comparando el valor del indicador con la tabla de distribuciones chi cuadrado dentro del nivel apropiado de grados de libertad.
El valor AIC tambien se calcula usando C, pero no vamos a entrar en tantos detalles.
La sección “Coefficients:”
Es quizá la parte mas tradicional de los resultados. Aquí se muestra la significancia de cada una de las relaciones que has establecido entre la variable predictora y la variable de respuesta. Los valores estimados estandarizados se usan cuando quieres comparar el peso de las variables predictoras.
La última sección “Individual R-squared:”
Aquí se muestra los coeficientes de correlación de cada componente, Piecewise es especialmente útil cuando quieres hacer análisis anidados. En ese caso no existe realmente un R^2 tradicional, y se utiliza metodos como. R^2 marginal o compuesto.
El paquete piecewise usa un motor grafico extraido de DiagrammeR.
plot(model) En caso quieras comparar tus resultados con los que aparecen en esta web, podemos fijar la producción de valores aleatorios en un número fijo:
set.seed(123)
dat <- data.frame(Productividad = runif(50),
fotosintesis = runif(50),
cierre_estomatico = runif(50),
Sequia = runif(50))
summary(dat)## Productividad fotosintesis cierre_estomatico Sequia
## Min. :0.02461 Min. :0.0006248 Min. :0.01047 Min. :0.1111
## 1st Qu.:0.27137 1st Qu.:0.2259940 1st Qu.:0.25960 1st Qu.:0.3296
## Median :0.50295 Median :0.4453582 Median :0.50487 Median :0.4999
## Mean :0.52009 Mean :0.4770270 Mean :0.51539 Mean :0.5131
## 3rd Qu.:0.78084 3rd Qu.:0.7425265 3rd Qu.:0.73584 3rd Qu.:0.6569
## Max. :0.99427 Max. :0.9849570 Max. :0.98422 Max. :0.9856Cualquier persona que ejecute estas lineas de codigo deberian tener el mismo resultado. Ahora tu modelo deberia ser igual al siguiente:
model <- psem(
lm(Productividad ~ fotosintesis, dat), # Componente 1
lm(fotosintesis ~ cierre_estomatico, dat), # Componente 2
lm(cierre_estomatico ~ Sequia, dat) # Componente 3
)
summary(model, .progressBar = F)##
## Structural Equation Model of model
##
## Call:
## Productividad ~ fotosintesis
## fotosintesis ~ cierre_estomatico
## cierre_estomatico ~ Sequia
##
## AIC BIC
## 24.267 41.475
##
## ---
## Tests of directed separation:
##
## Independ.Claim Test.Type DF Crit.Value P.Value
## fotosintesis ~ Sequia + ... coef 47 -1.6317 0.1094
## Productividad ~ Sequia + ... coef 47 -0.8225 0.4150
## Productividad ~ cierre_estomatico + ... coef 46 0.0513 0.9593
##
## Global goodness-of-fit:
##
## Fisher's C = 6.267 with P-value = 0.394 and on 6 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value
## Productividad fotosintesis 0.0834 0.1531 48 0.5451 0.5882
## fotosintesis cierre_estomatico 0.1574 0.1339 48 1.1750 0.2458
## cierre_estomatico Sequia -0.1148 0.1817 48 -0.6318 0.5305
## Std.Estimate
## 0.0784
## 0.1672
## -0.0908
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method R.squared
## Productividad none 0.01
## fotosintesis none 0.03
## cierre_estomatico none 0.01¿Qué pasa cuando el modelo no es exitoso?
Al igual que los modelos creados con Lavaan, en piecewise empezamos un proceso de “podado” o ajuste de las variables. En el caso de Lavaan, lo ideal seria ejecutar modindices y revisar que relaciones causales (~) o covarianzas (~~) hemos dejado de lado y deberian incluirse. En el caso de Piecewise usamos los resultados del tests of directed separation. Aqui aparecen las relaciones significativas que no hemos incluido pero explican la variabilidad de nuestros datos. Analizamos si existe un fundamento teorico para incluirlos en el modelo, e intentamos crear un nuevo modelo que incluya dichas relaciones.
Nuestro modelo final sera exitoso solo cuando:
- El coeficience de Ficher C sea mayor de 0.05
- La significacion de todos los coeficientes del modelo sea menor de 0.05
- La significacion de todos los coeficientes del test de directed separation sea mayor de 0.05
Es facil confundirse al principio, pero la experiencia hace al maestro.
¿Qué tal si intentamos crear un modelo Piecewise con los datos de Grace Keeley?
En la clase anterior exploraron los datos de Grace keeley para construir un modelo que explica la diversidad a partir de factores abioticos. La base de datos de keeley esta incluida dentro del paquete Piecewise, asi que pueden copiar el siguiente código para que aparezca el dataframe dentro el ecosistema de R.
data(keeley)
keeley
names(keeley)En teoria, en la clase anterior tendrian que haber construido el siguiente modelo usando Lavaan. ¿Pueden construirlo ahora usando Piecewise?
Análisis de datos anidados con Piecewise model
Una de las grandes ventajas de usar Lavaan es la posibilidad de crear variables latentes. De hecho el ejemplo de Keeley tiene mas de una variable latente que se suele omitir en los ejemplos por condiciones practicas. En piecewise no es posible crear variables latentes pero permite tener en cuenta otro tipo de diseños experimentales: los datos anidados.
En un experimento X es bastante común que los datos se agrupen de manera espacial, temporal o geográfica. Estas formas de agrupación suelen estar ligadas a variaciones en la distribución del error de las variables, y se suelen considerar como “random effects”.
Por ejemplo, cuando estas midiendo las caracteristicas de la vegetación dentro de una serie de parcelas dispuestas a lo largo de un gradiente ambiental. Las características de la vegetación cambiarán según las condiciones ambientales, pero pertenecer a una parcela determinada, incluso si las condiciones ambientales son las mismas, puede alterar la varibilidad de tus datos.
Otro ejemplo se da cuando un equipo de personas estan tomando datos en un experimento. Imagina un equipo compuesto por Rafael y Lucia. Ambos estan a cargo de un experimento de invernadero donde buscan medir el efecto de la fertilización sobre el crecimiento. Rafa toma datos los lunes, miercoles y viernes. Lucia toma datos los martes, jueves y sabado. Sin embargo, al analizar los datos te das cuenta que según la persona que haya tomado los datos, la relación de las variables es diferente. Esta claro que no trabajan igual, pero es dificil saber quien lo hace mejor o peor porque eso implicaría que el investigador favorezca la hipótesis nula o la hipótesis alternativa de su experimento. Para solucinar esto podemos anidar los datos segun la persona que los tomó.
La inclusión de diseños anidados en Piecewise es mucho mas realista y práctica en comparación a Lavaan. Nos permite modelar nuestros datos de la misma manera en la que diseñamos nuestros experimentos.
Como ejemplo real, os muestro el analisis de los datos de Shipley (2009) que tambien estan incluidos dentro del paquete Piecewise:
Aqui se intentó modelar la supervivencia de los arboles a partir de su crecimiento, fecha de siembra, dias acumulados de calor, y latitud.
Al analizarlo con lavan podemos ver que:
data("shipley")
shipleylibrary(lavaan)## This is lavaan 0.6-8
## lavaan is FREE software! Please report any bugs.shipley_model <- '
DD ~ lat
Date ~ DD
Growth ~ Date
Live ~ Growth
'
shipley_sem <- sem(shipley_model, shipley)
summary(shipley_sem, standardize = T, rsq = T)## lavaan 0.6-8 ended normally after 27 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 8
##
## Used Total
## Number of observations 1431 1900
##
## Model Test User Model:
##
## Test statistic 38.433
## Degrees of freedom 6
## P-value (Chi-square) 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## DD ~
## lat -0.860 0.023 -37.923 0.000 -0.860 -0.708
## Date ~
## DD -0.517 0.016 -32.525 0.000 -0.517 -0.652
## Growth ~
## Date 0.173 0.020 8.508 0.000 0.173 0.219
## Live ~
## Growth 0.006 0.001 9.854 0.000 0.006 0.252
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .DD 52.628 1.967 26.749 0.000 52.628 0.499
## .Date 38.080 1.424 26.749 0.000 38.080 0.575
## .Growth 38.981 1.457 26.749 0.000 38.981 0.952
## .Live 0.025 0.001 26.749 0.000 0.025 0.936
##
## R-Square:
## Estimate
## DD 0.501
## Date 0.425
## Growth 0.048
## Live 0.064Podemos ver que el modelo no es significativo, y los R^2 de los componentes tampoco son muy buenos, aunque las regresiones si lo son. Ahora vamos a intentar ejecutar el mismo modelo con Piecewise, considerando que los arboles no estan dispuestos al azar sino anidados en poblaciones.
library(nlme)
library(lme4)## Loading required package: Matrix##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmListshipley_psem <- psem(
lme(DD ~ lat, random = ~ 1 | site / tree, na.action = na.omit,
data = shipley),
lme(Date ~ DD, random = ~ 1 | site / tree, na.action = na.omit,
data = shipley),
lme(Growth ~ Date, random = ~ 1 | site / tree, na.action = na.omit,
data = shipley),
glmer(Live ~ Growth + (1 | site) + (1 | tree),
family = binomial(link = "logit"), data = shipley)
)
summary(shipley_psem, .progressBar = FALSE)## boundary (singular) fit: see ?isSingular## boundary (singular) fit: see ?isSingular##
## Structural Equation Model of shipley_psem
##
## Call:
## DD ~ lat
## Date ~ DD
## Growth ~ Date
## Live ~ Growth
##
## AIC BIC
## 49.536 149.592
##
## ---
## Tests of directed separation:
##
## Independ.Claim Test.Type DF Crit.Value P.Value
## Date ~ lat + ... coef 18 -0.0798 0.9373
## Growth ~ lat + ... coef 18 -0.8929 0.3837
## Live ~ lat + ... coef 1431 1.0280 0.3039
## Growth ~ DD + ... coef 1329 -0.2967 0.7667
## Live ~ DD + ... coef 1431 1.0046 0.3151
## Live ~ Date + ... coef 1431 -1.5617 0.1184
##
## Global goodness-of-fit:
##
## Fisher's C = 11.536 with P-value = 0.484 and on 12 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## DD lat -0.8355 0.1194 18 -6.9960 0 -0.6877 ***
## Date DD -0.4976 0.0049 1330 -100.8757 0 -0.6281 ***
## Growth Date 0.3007 0.0266 1330 11.2917 0 0.3824 ***
## Live Growth 0.3479 0.0584 1431 5.9552 0 - ***
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method Marginal Conditional
## DD none 0.49 0.70
## Date none 0.41 0.98
## Growth none 0.11 0.84
## Live delta 0.16 0.18Las estrategias de análisis de los modelos estructurales cambian con el tiempo. Lo mas importante es el proceso de construcción. Concebir una idea que sea única, útil y facil de entender. Vamos a crear un matriz aleatoria de variables, asignarles nombres de procesos ecologicos, y tratar de establecer relaciones entre ellas.
library(corrplot)## Warning: package 'corrplot' was built under R version 4.1.3## corrplot 0.92 loadedset.seed(123)
db <- data.frame(replicate(3,sample(10:100,100,rep=TRUE)))
rd <- data.frame(sample(-5:40,100,rep=TRUE),
sample(-3:30,100,rep=TRUE),
sample(-10:20,100,rep=TRUE),
sample(15:40,100,rep=TRUE),
sample(4:20,100,rep=TRUE),
sample(-10:0,100,rep=TRUE),
sample(10:40,100,rep=TRUE),
sample(8:35,100,rep=TRUE),
sample(-7:20,100,rep=TRUE))
db2 <- cbind(db + 2*rd[1:3], db + 4*rd[4:6], db*-1 + 3*rd[7:9], db*-1+ rd[1:3] + rd[4:6], db + 5*rd[1:3])
db2 <- db2[,sample(ncol(db2))]
colnames(db2) <- c("primary_production", "respiration", "livestock", "carbon_emissions", "Soil_N",
"Seed_dispersal", "forest_regeneration", "Litter_decomposition", "climate_sesonality", "Microbial_activity",
"evapotranspiration", "precipitation", "Biodiversity", "lumber_production", "pollination")
m<- cor(db2)
corrplot(m, diag = FALSE)