Gustavo Carvalho
2010-Sep-16 13:50 UTC
[R] problems trying to reproduce structural equation model using the sem package
Hello, I've been unsuccessfully trying to reproduce a sem from Grace et al. (2010) published in Ecological Monographs: http://www.esajournals.org/doi/pdf/10.1890/09-0464.1 The model in question is presented in Figure 8, page 81. The errors that I've been getting are: 1. Using a correlation matrix: res.grace <- sem(grace.model, S = grace, N = 190) Warning message: In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars, ?: ?Could not compute QR decomposition of Hessian. Optimization probably did not converge. 2. Using a variances/covariances matrix: res.grace <- sem(grace.model, S = grace.cov, N = 190) Error in solve.default(C) : ?Lapack routine dgesv: system is exactly singular In addition: Warning messages: 1: In log(det(C)) : NaNs produced 2: In sem.default(ram = ram, S = S, N = N, param.names = pars, var.names = vars, ?: ? singular Hessian: model is probably underidentified. (...) So far I've tried: 1. Fixing the variances of the latent variables 2. Allowing the exogenous indicators to covary (fixed.x parameter in sem()) 3. Manually inserting the published parameter estimates during model specification (specify.model()) to see if the starting parameters passed to nlm were the problem 4. Extensively looking for typing mistakes Anyway, there seems to be a problem either with the way I specified the model or with the model itself as it has been published. I can see that the number of degrees of freedom in ?the model that I've specified is 21, as in the published model. Any light you could shed on this would be greatly appreciated. The code to reproduce all steps is presented below. Thank you very much, Gustavo. ########## library(sem) grace <- matrix(ncol = 10, nrow = 10) variables <- c("lightlog", "light", "dstb", "species_count", "masslog", ? ? ? ? ? ? ? "soil_carbon", "soil_organic", "soil_low_flooding", ? ? ? ? ? ? ? "soil_high_flooding", "soil_salinity") rownames(grace) <- colnames(grace) <- variables diag(grace) <- 1 ## Coefficients from the paper. grace.coefficients <- c(0.858, 0.667, -0.251, -0.699, 0.06, 0.012, 0.552, 0.547, ? ? ? ? ? ? ? ? ? ? ? ?0.327, 0.776, -0.404, -0.794, 0.157, 0.120, 0.439, 0.462, ? ? ? ? ? ? ? ? ? ? ? ?0.321, -0.228, -0.686, 0.218, 0.186, 0.249, 0.290, 0.216, ? ? ? ? ? ? ? ? ? ? ? ?0.291, 0.119, 0.132, -0.374, -0.406, -0.292, -0.096, -0.071, ? ? ? ? ? ? ? ? ? ? ? ?-0.426, -0.466, -0.138, 0.973, -0.170, -0.150, 0.249, -0.211, ? ? ? ? ? ? ? ? ? ? ? ?-0.188, 0.244, 0.959, 0.073, 0.052) grace.sds <- c(1.11, 0.285, 3.29, 3.33, 1.44, 0.605, 1.23, 1.33, 1.27, 1.68) grace[lower.tri(grace, diag = F)] <- grace.coefficients grace[upper.tri(grace, diag = F)] <- t(grace)[upper.tri(grace, diag = F)] ## Covariances matrix grace.cov <- outer(grace.sds, grace.sds) * grace ## Specifying the model grace.model <- specify.model() salinity -> soil_salinity, NA, 1 flooding -> soil_high_flooding, NA, 1 flooding -> soil_low_flooding, flooding_low_flooding, NA infertility -> soil_organic, NA, -1 infertility -> soil_carbon, inf_carbon, NA disturbance -> dstb, NA, 1 biomass -> masslog, NA, 1 light -> light_effect, NA, 1 lightlog -> light_effect, lightlog_light_effect, NA richness -> species_count, NA, 1 salinity -> richness, salinity_richness, NA salinity -> light, salinity_light, NA flooding -> richness, flooding_richness, NA flooding -> biomass, flooding_biomass, NA infertility -> richness, infertility_richness, NA disturbance -> biomass, disturbance_biomass, NA disturbance -> light, disturbance_light, NA biomass -> light, biomass_light, NA light_effect -> richness, light_effect_richness, NA salinity <-> flooding, salinity_flooding, NA salinity <-> infertility, salinity_infertility, NA salinity <-> disturbance, salinity_disturbance, NA flooding <-> infertility, flooding_infertility, NA flooding <-> disturbance, flooding_disturbance, NA infertility <-> disturbance, infertility_disturbance, NA biomass <-> biomass, biomass, NA masslog <-> masslog, masslog, NA light <-> lightlog, light_lightlog, NA light_effect <-> light_effect, NA, 0 richness <-> richness, richness, NA species_count <-> species_count, species_count, NA light <-> light, light, NA lightlog <-> lightlog, NA, 1 salinity <-> salinity, salinity, NA disturbance <-> disturbance, disturbance, NA flooding <-> flooding, flooding, NA infertility <-> infertility, infertiliy, NA soil_salinity <-> soil_salinity, soil_salinity, NA soil_high_flooding <-> soil_high_flooding, soil_high, NA soil_low_flooding <-> soil_low_flooding, soil_low, NA soil_organic <-> soil_organic, soil_organic, NA soil_carbon <-> soil_carbon, soil_carbon, NA dstb <-> dstb, dstb, NA res.grace <- sem(grace.model, S = grace.cov, N = 190)