Práctico 3. Análisis de senderos y modelado de ecuaciones estructurales
Introducción
Utilizaremos los datos de Carmona (2006) sobre la evolución fenotípica de la resistencia de Datura stramonium contra su herbívoro folívoro Epitrix sp., los cuales se encuentran contenidos en el archivo datura1.txt. Para realizar el análisis utilizaremos el paquete lavaan de R. Como pre-procesamiento de los datos, vamos a estandarizar todal las variables.
library(lavaan)
datura <- read.table("datura1.txt", header = T)
datos <- data.frame(scale(datura))Hipótesis causal
Utilizaremos el siguiente modelo como hipótesis para explicar la relación entre el daño por herbívoros, el crecimiento, el tiempo hasta el inicio de la floración, y la producción de frutos (adecuación) en Datura stramonium. Cada flecha representa una relación causal: la variable de origen de la flecha es la supuesta causa, la variable que es señalada por la flecha es la variable de respuesta. Las flechas dobles se usan para indicar varianzas o covarianzas.
Modelado de rutas mediante regresiones.
Los coeficientes de ruta son simbolizados con la letra rho (\(\rho\)). Por ejemplo, en la relación daño → crecimiento, el coeficiente de ruta puede representarse como \(\rho_{cre, da}\). Los coeficientes de ruta pueden estimarse a partir de una serie de regresiones múltiples, de donde los coeficientes parciales (i.e., la relación entre dos variables habiendo “controlado” el efecto de las demás) son los coeficientes de ruta.
# sólo daño afecta a crecimiento
m1<- lm(crec ~ dano, data = datos)
summary(m1)##
## Call:
## lm(formula = crec ~ dano, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.60758 -0.56579 -0.09768 0.54944 2.28180
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.13497 0.05980 2.257 0.0248 *
## dano -0.06152 0.05992 -1.027 0.3055
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9698 on 261 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.004022, Adjusted R-squared: 0.0002065
## F-statistic: 1.054 on 1 and 261 DF, p-value: 0.3055# daño y crecimiento afectan a floración
m2 <- lm(flor.dias ~ dano + crec, data = datos)
summary(m2)##
## Call:
## lm(formula = flor.dias ~ dano + crec, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.13365 -0.66284 -0.06507 0.64311 2.97450
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.04909 0.05904 -0.831 0.4065
## dano -0.13624 0.05870 -2.321 0.0211 *
## crec -0.04466 0.06052 -0.738 0.4613
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9483 on 260 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.02159, Adjusted R-squared: 0.01406
## F-statistic: 2.868 on 2 and 260 DF, p-value: 0.0586# finalmente daño, crecimiento y floración afectan a la adecuación
m3 <- lm(frutos ~ dano + crec + flor.dias, data = datos)
summary(m3)##
## Call:
## lm(formula = frutos ~ dano + crec + flor.dias, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6697 -0.5501 -0.1557 0.3681 3.8984
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01466 0.05470 -0.268 0.789
## dano 0.07624 0.05487 1.389 0.166
## crec 0.40442 0.05606 7.215 5.96e-12 ***
## flor.dias -0.39837 0.05738 -6.943 3.08e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8774 on 259 degrees of freedom
## (42 observations deleted due to missingness)
## Multiple R-squared: 0.2958, Adjusted R-squared: 0.2876
## F-statistic: 36.26 on 3 and 259 DF, p-value: < 2.2e-16Con estos datos podemos colocar todos coeficientes de ruta
Coeficientes de efectos
Los coeficientes de efectos sobre una variable dependiente son calculados para estimar el efecto total de una o varias causas. Típicamente son referidas con la letra r. Un coeficiente de efecto tiene 4 componentes: efectos causales (directos e indirectos) y efectos no causales (espurios y “no analizables”).
- Efecto directo: es el coeficiente de ruta.
- Efecto indirecto: es el producto de los coeficientes de ruta involucrados.
- Efectos espurios: ocurren cuando hay una variable causal común para las variables sobre las que se calcula un componente de efecto indirecto.
Cálculo para el modelo de Datura * Daño sobre crecimiento. Directo. -0.062
* Daño sobre floración. Directo + Indirecto. -0.136 + (-0.062 × -0.045) = -0.133
* Crecimiento sobre floración: Directo. -0.045
* Crecimiento sobre adecuación. Directo + Indirecto. 0.404 + (-0.045 × -0.076) = 0.407
* Floración sobre adecuación. Directo. -0.398.
* Daño sobre adecuación. Directo + Indirecto. 0.076 + (-0.062 × 0.404) + (-0.136 × -0.398) = 0.105. Tiene efectos espúreos: -0.062 × -0.045 × -0.398 = -0.001
Cálculo del componente no analizable
\[U = \sqrt{1-R^2}\] Por ejemplo, el efecto de variables desconocidas sobre la adecuación sería \(U_{fru} = \sqrt{1 – 0.296} = 0.839\)
Modelado de ecuaciones estructurales mediante lavaan.
En base al conocimiento del sistema de relaciones causales, pueden proponerse varios modelos alternativos y evaluar cuál de ellos presenta mejor ajuste. En este caso construimos tres modelos para Datura.
MODELO 1
# modelo saturado
sat <- '
frutos ~ crec + flor.dias + dano
flor.dias ~ crec + dano
crec ~ dano'
fit.sat <- sem(model = sat, data = datos, missing = "listwise",
fixed.x = F)
summary(fit.sat, standardized = T, rsq = T, fit.measures = T)## lavaan 0.6-7 ended normally after 12 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free parameters 10
##
## Used Total
## Number of observations 263 305
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 99.018
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1431.267
## Loglikelihood unrestricted model (H1) -1431.267
##
## Akaike (AIC) 2882.534
## Bayesian (BIC) 2918.256
## Sample-size adjusted Bayesian (BIC) 2886.551
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value RMSEA <= 0.05 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 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
## frutos ~
## crec 0.404 0.056 7.270 0.000 0.404 0.377
## flor.dias -0.398 0.057 -6.996 0.000 -0.398 -0.366
## dano 0.076 0.054 1.400 0.161 0.076 0.073
## flor.dias ~
## crec -0.045 0.060 -0.742 0.458 -0.045 -0.045
## dano -0.136 0.058 -2.334 0.020 -0.136 -0.143
## crec ~
## dano -0.062 0.060 -1.031 0.303 -0.062 -0.063
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .frutos 0.758 0.066 11.467 0.000 0.758 0.704
## .flor.dias 0.889 0.078 11.467 0.000 0.889 0.978
## .crec 0.933 0.081 11.467 0.000 0.933 0.996
## dano 0.996 0.087 11.467 0.000 0.996 1.000
##
## R-Square:
## Estimate
## frutos 0.296
## flor.dias 0.022
## crec 0.004MODELO 2
# modelo con efectos indirectos del daño,
# y sin efecto del crecimiento en la floración
exdano <- '
frutos ~ crec + flor.dias
flor.dias ~ dano
crec ~ dano'
fit.exdano <- sem(model=exdano, data=datos, missing="listwise",
fixed.x=FALSE)
summary(fit.exdano, standardized = TRUE, rsq = T, fit.measures = T)## lavaan 0.6-7 ended normally after 12 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free parameters 8
##
## Used Total
## Number of observations 263 305
##
## Model Test User Model:
##
## Test statistic 2.503
## Degrees of freedom 2
## P-value (Chi-square) 0.286
##
## Model Test Baseline Model:
##
## Test statistic 99.018
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.995
## Tucker-Lewis Index (TLI) 0.984
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1432.519
## Loglikelihood unrestricted model (H1) -1431.267
##
## Akaike (AIC) 2881.037
## Bayesian (BIC) 2909.614
## Sample-size adjusted Bayesian (BIC) 2884.251
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.031
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.130
## P-value RMSEA <= 0.05 0.493
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.028
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## frutos ~
## crec 0.399 0.056 7.168 0.000 0.399 0.375
## flor.dias -0.410 0.057 -7.247 0.000 -0.410 -0.379
## flor.dias ~
## dano -0.133 0.058 -2.289 0.022 -0.133 -0.140
## crec ~
## dano -0.062 0.060 -1.031 0.303 -0.062 -0.063
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .frutos 0.764 0.067 11.467 0.000 0.764 0.719
## .flor.dias 0.891 0.078 11.467 0.000 0.891 0.980
## .crec 0.933 0.081 11.467 0.000 0.933 0.996
## dano 0.996 0.087 11.467 0.000 0.996 1.000
##
## R-Square:
## Estimate
## frutos 0.281
## flor.dias 0.020
## crec 0.004MODELO 3
# modelo con efecto del crecimiento en la floración
# y sin efecto del daño en la floración
exdano2<-'
frutos ~ crec + flor.dias
flor.dias ~ crec
crec ~ dano'
fit.exdano2 <- sem(model=exdano2, data=datos, missing="listwise",
fixed.x=FALSE)
summary(fit.exdano2, standardized=TRUE, rsq=T, fit.measures=T)## lavaan 0.6-7 ended normally after 12 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of free parameters 8
##
## Used Total
## Number of observations 263 305
##
## Model Test User Model:
##
## Test statistic 7.345
## Degrees of freedom 2
## P-value (Chi-square) 0.025
##
## Model Test Baseline Model:
##
## Test statistic 99.018
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.943
## Tucker-Lewis Index (TLI) 0.828
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1434.940
## Loglikelihood unrestricted model (H1) -1431.267
##
## Akaike (AIC) 2885.880
## Bayesian (BIC) 2914.457
## Sample-size adjusted Bayesian (BIC) 2889.093
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.101
## 90 Percent confidence interval - lower 0.030
## 90 Percent confidence interval - upper 0.183
## P-value RMSEA <= 0.05 0.102
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.060
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## frutos ~
## crec 0.399 0.056 7.164 0.000 0.399 0.372
## flor.dias -0.410 0.057 -7.243 0.000 -0.410 -0.376
## flor.dias ~
## crec -0.036 0.061 -0.589 0.556 -0.036 -0.036
## crec ~
## dano -0.062 0.060 -1.031 0.303 -0.062 -0.063
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .frutos 0.764 0.067 11.467 0.000 0.764 0.709
## .flor.dias 0.907 0.079 11.467 0.000 0.907 0.999
## .crec 0.933 0.081 11.467 0.000 0.933 0.996
## dano 0.996 0.087 11.467 0.000 0.996 1.000
##
## R-Square:
## Estimate
## frutos 0.291
## flor.dias 0.001
## crec 0.004REPORTE
Construcción de una tabla de resultados Es necesario reportar el \(\chi^2\) del modelo y su significancia. Existe mucha divergencia sobre qué índice de ajuste del modelo reportar, pero es conveniente reportar al menos uno de ellos. Finalmente, pude ser importante reportar algún criterio de información para comparar entre modelos.
| Modelo | \(\chi^2\) | df | P | SRMR | AIC |
|---|---|---|---|---|---|
| sat | 0 | 0 | - | 0 | 2882.53 |
| exdano | 2.503 | 2 | 0.286 | 0.028 | 2881.04 |
| exdano2 | 7.345 | 2 | 0.025 | 0.060 | 2885.88 |
Ejercicios
Utilizaremos los datos del trabajo de Dai & Galloway (2013) que puden descargarse libremente de Dryad. Los autores midieron tres rasgos fenotípicos en plantas de Passiflora incarnata: tamaño de la flor, área del estigma y deflección del estilo. El éxito de apareamiento fue definido como el número de granos de polen deposuitados en el estigma. El éxito reproductivo femenino fue estimado como el número total de semillas producidas por la planta. Finalmente, también se midió el tamaño de las plantas como una medida que puede afectar el éxito reproductivo. La hipótesis principal de este estudio es que los rasgos florales afectan el aparemiento (selección sexual) y a su vez, este influye sobre el éxito reproductivo, como se detalla en este diagrama.
- Reconstruya este análisis de senderos utilizando el set de datos proporcionado.
- Idenifique las variables endógenas y exógenas, las vias que corresponden a efectos y las vías que corresponden a correlaciones.
- ¿Puede proponer otra hipótesis causal para estas variables? Ajústela y compare la adecuación de los dos modelos.
Los datos de este ejercicio fueron obtenidos de:
Dai C & Galloway LF (2013) Sexual selection in a hermaphroditic plant through female reproductive success. Journal of Evolutionary Biology 26: 2622-2632.