[R] problems trying to reproduce structural equation model using the sem package

[Gustavo Carvalho]

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)