SEM
From InterSciWiki
- lavaan GitHub -- see also SEM 2sls
- OpenMx Advanced Structural Equation Modeling -
- OpenMx Wiki
setwd('/Users/drwhite/Documents/3sls/sccs/')
load("comb.Rdata", .GlobalEnv)
load("HiGod.Rdata", .GlobalEnv)
load(url('http://intersci.ss.uci.edu/wiki/R/HiGod.Rdata'), .GlobalEnv)
HiGod$AnimXbwealth
HiGod$HiGod4
WX<-read.csv(file = "WYWX2.csv")
save(WX, file = "WYWX2.Rdata")
load("WX.Rdata", .GlobalEnv)
WXAnimXbwealth<-WX$AnimXbwealth
WXSuperjhWriting<-WX$SuperjhWriting
WXNo_rain_Dry<-WX$No_rain_Dry
#load("WYWX.Rdata", .GlobalEnv) #shows variable names:
#WX.AnimXbwealth WX.ecorich WX.FxCmtyWages WX.SuperjhWriting WX.No_rain_Dry WX.Missions
#load(url('http://intersci.ss.uci.edu/wiki/R/WYWX.Rdata'), .GlobalEnv)
#WXAnimXbwealth<-WYWX$WX.AnimXbwealth
#WXSuperjhWriting<-WYWX$WX.SuperjhWriting
#WXAnimXbwealth
library(sem)
#load("HiGodScal.Rdata", .GlobalEnv)
#load("WYhat.Rdata", .GlobalEnv) #from the two_stage_ols.R run with save to WYhat.Rdata
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting, instruments=~ WXAnimXbwealth, WXSuperjhWriting,WXNo_rain_Dry,data=HiGod)
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting, instruments=~ WXAnimXbwealth, WXSuperjhWriting,WXNo_rain_Drydata=HiGod))
#WX.AnimXbwealth WX.ecorich WX.FxCmtyWages WX.SuperjhWriting WX.No_rain_Dry WX.Missions
Error in chol.default(XtZ %*% invZtZ %*% t(XtZ)) :
table(HiGod$AnimXbwealth,round(WXAnimXbwealth,1))
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting, instruments=~ WXAnimXbwealth, WXSuperjhWriting,WXNo_rain_Dry,data=HiGod)
Error in model.frame.default(formula = HiGod4 ~ AnimXbwealth + SuperjhWriting + : attempt to apply non-function
table(HiGod$AnimXbwealth,WXAnimXbwealth) #test that WX var differs from var
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.5 2.6
0 2 23 4 7 17 4 1 1 2 4 7 8 4 2 3 0 3 3 3 0 0 0 0 2 2
1 2 0 0 1 0 0 0 0 1 0 1 2 0 0 0 0 0 1 0 2 0 0 0 0 0
2 0 0 1 0 0 0 0 0 0 2 5 1 0 0 1 1 0 0 0 0 0 0 0 0 2
3 0 1 0 0 0 0 0 0 2 1 1 0 0 0 0 0 0 2 0 0 1 1 0 0 0
4 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 1
5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
8 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1
2.7 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9
0 1 1 1 1 0 1 0 4 0 0 1 2 1
1 1 2 0 0 0 0 0 0 0 0 0 0 0
2 3 0 0 0 0 0 0 0 0 0 0 0 0
3 1 2 1 2 0 0 0 0 0 0 0 0 0
4 2 0 1 1 2 0 0 1 1 1 0 0 0
5 0 0 0 0 0 0 1 0 0 2 0 0 0
6 0 0 0 1 0 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0 0 0 0 0 0
9 0 0 0 1 1 0 0 0 0 0 0 0 0
10 0 0 0 0 0 0 0 1 1 0 0 0 0
the leading minor of order 3 is not positive definite
- Yves_Rosseel Lavaan tsls new - qgraph
- Structural equation modeling and Structural equation modeling
Contents |
Lavaan references
Hu, Li‐tze, and Peter M. Bentler. 1999. Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling: A Multidisciplinary Journal 6(1): 1-55.
MCMC
Wikipedia:OpenBUGS OpenBUGS is a computer software for the Bayesian analysis of complex statistical models using Markov chain Monte Carlo (MCMC) methods. OpenBUGS is the open source variant of WinBUGS (Bayesian inference Using Gibbs Sampling). It runs under Windows and Linux, as well as from inside the R statistical package. Versions from v3.0.7 onwards have been designed to be at least as efficient and reliable as WinBUGS over a range of test applications.
Wall, Melanie M., and Xuan Liu. 2009. Spatial latent class analysis model for spatially distributed multivariate binary data· Journal: Computational Statistics & Data Analysis archive Volume 53 Issue 8:3057-3069. Authors: at University of Minnesota, Division of Biostatistics, United States
- Abstract. A spatial latent class analysis model that extends the classic latent class analysis model by adding spatial structure to the latent class distribution through the use of the multinomial probit model is introduced. Linear combinations of independent Gaussian spatial processes are used to develop multivariate spatial processes that are underlying the categorical latent classes. This allows the latent class membership to be correlated across spatially distributed sites and it allows correlation between the probabilities of particular types of classes at any one site. The number of latent classes is assumed to be fixed but is chosen by model comparison via cross-validation. An application of the spatial latent class analysis model is shown using soil pollution samples where 8 heavy metals were measured to be above or below government pollution limits across a 25 square kilometer region. Estimation is performed within a Bayesian framework using MCMC and is implemented using the OpenBUGS software.
Xuan Liu, Melanie M. Wall and James S. Hodges. 2005. Generalized spatial structural equation models. Biostatistics Volume6, Issue4Pp. 539-557. - Oxford JournalsLife Sciences & Mathematics & Physical Sciences.
- Abstract. It is common in public health research to have high-dimensional, multivariate, spatially referenced data representing summaries of geographic regions. Often, it is desirable to examine relationships among these variables both within and across regions. An existing modeling technique called spatial factor analysis has been used and assumes that a common spatial factor underlies all the variables and causes them to be related to one another. An extension of this technique considers that there may be more than one underlying factor, and that relationships among the underlying latent variables are of primary interest. However, due to the complicated nature of the covariance structure of this type of data, existing methods are not satisfactory. We thus propose a generalized spatial structural equation model. In the first level of the model, we assume that the observed variables are related to particular underlying factors. In the second level of the model, we use the structural equation method to model the relationship among the underlying factors and use parametric spatial distributions on the covariance structure of the underlying factors. We apply the model to county-level cancer mortality and census summary data for Minnesota, including socioeconomic status and access to public utilities.
- Key words: Bayesian Factor analysis Latent Linear model of coregionalization (LMC) Markov chain Monte Carlo (MCMC) Spatial Structural equation models (SEM)
- THIS ONE U Minnesota Biochemistry
- NOT THIS XUAN LIU Xuan Liu Biochem UC Riverside.
SEM & autocorrelation?
Course in SEM that includes autocorrelation.
- 2-level modeling Structural Equation Modeling: Foundations and Extensions (Advanced Quantitative Techniques in the Social Sciences) by David W. Kaplan (Jul 23, 2008)
http://www.guilford.com/cgi-bin/cartscript.cgi?page=pr/kline.htm&dir=research/MSS_series&cart_id=
Kline
http://www.guilford.com/cgi-bin/cartscript.cgi?page=pr/kline.htm&dir=research/MSS_series&cart_id=
Build data for Fox SEM of HiGods following model of Klein
- SEM - John Fox - you need at least as many instrumental variables as coefficients to estimate sem equations - JF
- Lisrel
Klein
library(sem) data(Klein) Klein$P.lag <- c(NA, Klein$P[-22]) Klein$X.lag <- c(NA, Klein$X[-22]) Klein.eqn1 <- tsls(C ~ P + P.lag + I(Wp + Wg), instruments=~G + T, data=Klein)
summary(tsls(C ~ P + P.lag + I(Wp + Wg), instruments=~1 + G + T + Wg + I(Year - 1931) + K.lag + P.lag + X.lag, data=Klein))
#It looks like these two variables are being lagged so they are not essential
summary(tsls(I ~ P + P.lag + K.lag, instruments=~1 + G + T + Wg + I(Year - 1931) + K.lag + P.lag + X.lag, data=Klein))
HiG
setwd('/Users/drwhite/Documents/3sls/sccs/')
library(sem)
load("HiGodScal.Rdata", .GlobalEnv)
load(url('http://intersci.ss.uci.edu/wiki/R/HiGod.Rdata'), .GlobalEnv)
load("WYhat.Rdata", .GlobalEnv) #from the two_stage_ols.R run with save to WYhat.Rdata
HiGodScal$HiGod <- c(NA, HiGodScal$No_rain[-186]) #works, HiGod$HiGod <- c(NA, No_rain[-186]) #does not
HiGodScal$WY=WYhat
#WY=WYhat
save(HiGodScal, file="HiGodScal.Rdata", .GlobalEnv)
load("HiGodScal.Rdata", .GlobalEnv)
save(HiGodScal, file="HiGodScal.Rdata", .GlobalEnv)
HiG=HiGodScal
save(HiG, file="HiG.Rdata", .GlobalEnv)
save(HiGodScal, file="HiGodScal.Rdata", .GlobalEnv)
data(HiG) #NOT WORKING
read.csv(file="FxCmtyWage.csv",header=TRUE) Printable_Moral_Gods_Model_&_Data#Retrieve_the_imputed_FxCmtyWages_and_create_a_missing_data_dummy_for_FxCmtyWages
load("WYWX.Rdata", .GlobalEnv)
#Doug: You can look inside an object like WYWX using the str() command: Str(WYWX)
#if WYWX is an lm object (i.e., WYWX<-lm(WY~WX) ), then you can access your original variable WX.milk using: WYWX$model$WX.milk
somevar = read.csv(file='WYWX.csv', header=TRUE)
HG <- tsls(HiGod4 ~ No_rain_Dry + AnimXbwealth + SuperjhWriting + Missions, instruments=~ WY, data=HiG)
Error in chol.default(XtZ %*% invZtZ %*% t(XtZ)) : is not positive definite ERROR the leading minor of order 3
HG <- tsls(HiGod4 ~ No_rain_Dry + AnimXbwealth, instruments=~ WY, data=HiG) #ERROR leading minor of order 3
HG <- tsls(HiGod4 ~ No_rain_Dry, instruments=~ WY, data=HiG) # THIS WORKS but no second variable works!
HG <- tsls(HiGod4 ~ No_rain_Dry + SuperjhWriting, instruments=~ WY, data=HiG) #no second variable works!
summary(tsls(HiGod4 ~ No_rain_Dry, instruments=~WY, data=HiG))
summary(tsls(HiGod4 ~ No_rain_Dry + AnimXbwealth + SuperjhWriting + Missions, instruments=WY, data=HiG))
Model Formula: HiGod4 ~ No_rain_Dry
Instruments: ~WY
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-4.85 -4.85 -4.00 0.00 4.09 12.20
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.734e-15 0.4761 -7.843e-15 1.0000
No_rain_Dry -6.028e+00 39.8043 -1.514e-01 0.8798
HG <- tsls(HiGod4 ~ Missions, instruments=~ WY, data=HiG) # THIS WORKS but no second variable works!
summary(tsls(HiGod4 ~ Missions, instruments=~WY, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.488e-15 0.1159 2.147e-14 1.0000
Missions -1.001e+00 1.6081 -6.222e-01 0.5346
tsls = 2SLS
- I had commas here!!! Not pluses
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting, instruments=~ WY, data=HiG) # DOESNT WORK leading minor of order 2 summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting, instruments=~WY, data=HiG)) Error in chol.default(XtZ %*% invZtZ %*% t(XtZ)) : the leading minor of order 2 is not positive definite
- Dear Doug,
- Do you have at least as many instrumental variables as coefficients to estimate?
- John
- w_inv<-solve(ww) SCCS_R_package#neff_almost_working As per invertibiliy
- try multiple IVs
Logdate=as.numeric(Logdate)
Hig$Logdate=as.numeric(Logdate)
HG<-tsls(HiGod4~No_rain_Dry, instruments=~WY, data=HiG)
summary(tsls(HiGod4~No_rain_Dry, instruments=~WY+Logdate+Missions, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.366e-15 0.5519 -7.912e-15 1.0000
No_rain_Dry -7.068e+00 14.9915 -4.715e-01 0.6378
HG<-tsls(HiGod4~AnimXbwealth, instruments=~WY, data=HiG)
summary(tsls(HiGod4~AnimXbwealth, instruments=~WY+Logdate+Missions, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.259e-16 0.0712 1.019e-14 1.0000
AnimXbwealth 6.921e-01 1.7527 3.949e-01 0.6934
HG<-tsls(HiGod4~FxCmtyWages, instruments=~WY, data=HiG)
summary(tsls(HiGod4~FxCmtyWages, instruments=~WY+Logdate+Missions, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03865 0.1156 -0.3344 0.73870
FxCmtyWages 0.81813 0.3591 2.2784 0.02464
John: HG<-tsls(HiGod4~FxCmtyWages, instruments=~WY, data=HiG)
summary(tsls(HiGod4~FxCmtyWages, instruments=~WY, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03865 0.2305 -0.1676 0.8672
FxCmtyWages 2.35395 4.7289 0.4978 0.6196 <-- non.signif give 1 IV
HG<-tsls(HiGod4~FxCmtyWages, instruments=~WY+Logdate+Missions, data=HiG)
summary(tsls(HiGod4~FxCmtyWages, instruments=~WY+Logdate+Missions, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03865 0.1156 -0.3344 0.73870
FxCmtyWages 0.81813 0.3591 2.2784 0.02464 <-- signif given 3 IVs
lm(HiGod4 ~ FxCmtyWages+Logdate+Missions + WY)
summary(lm(HiGod4 ~ FxCmtyWages+Logdate+Missions + WY))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.47991 0.52300 0.918 0.361
FxCmtyWages 0.11408 0.09962 1.145 0.255 <--looks like my WY isn't right
Logdate 0.37461 0.45693 0.820 0.414
Missions 0.48978 0.52265 0.937 0.351
WY -0.40699 0.37495 -1.085 0.280 <--looks like my WY extract isn't right <-- DOES IT CORRESPOND TO 2SLS-IS??
2SLS-IS
coef xrange dx dy/dx Fstat ddf sig.value
(Intercept) 0.027 NA NA NA 0.188 1056800.384 0.66468665 distance 1.148 2.000 2.295 0.899 66.098 204522.625 0.00000000 FxCmtyWages 0.072 2.914 0.211 0.083 1.128 1282.762 0.28841394 <-- agrees Logdate 0.149 12.591 1.871 0.733 5.367 96155.897 0.02052832 <-- does not Missions 0.155 2.042 0.317 0.124 5.536 48594.932 0.01863034 <-- does not
HG<-tsls(HiGod4~FxCmtyWages, instruments=~WY, data=HiG)
summary(tsls(HiGod4~FxCmtyWages, instruments=~WY, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.03865 0.2305 -0.1676 0.8672
FxCmtyWages 2.35395 4.7289 0.4978 0.6196 <--
lm(HiGod4 ~ FxCmtyWages + WY)
summary(lm(HiGod4 ~ FxCmtyWages + WY))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.52936 0.51767 1.023 0.309
FxCmtyWages 0.14298 0.09695 1.475 0.143
WY -0.41801 0.37430 -1.117 0.267 <--
HG <- tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages, instruments=~WY + Logdate, data=HiG) summary(tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages, instruments=~WY + Logdate, data=HiG)) #N.Sig. p(F)=.28
HG<-tsls(HiGod4~No_rain_Dry+FxCmtyWages+AnimXbwealth, instruments=~WY+Logdate+Missions, data=HiG)
summary(tsls(HiGod4~No_rain_Dry+FxCmtyWages+AnimXbwealth, instruments=~WY+Logdate+Missions, data=HiG)) #p(F)=.28
HG<-tsls(HiGod4~No_rain_Dry+FxCmtyWages+SuperjhWriting, instruments=~WY+Logdate+Missions, data=HiG)
summary(tsls(HiGod4~No_rain_Dry+FxCmtyWages+SuperjhWriting, instruments=~WY+Logdate+Missions, data=HiG)) #n.sig
HG <- tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages, instruments=~WY + Logdate, data=HiG)
summary(tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages, instruments=~WY + Logdate, data=HiG)) #N.Sig. p(F)=.28
HG <- tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages, instruments=~WY + ww %*% Logdate, data=HiG)
summary(tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages, instruments=~WY + ww %*% Logdate, data=HiG)) #N.Sig. p(N)=.28
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + Missions, instruments=~WY + Logdate, data=HiG) # WORKS
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + Missions, instruments=~WY + Logdate, data=HiG)) #but P=1.0
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + FxCmtyWages, instruments=~WY + Missions, data=HiG) #WORKS
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + FxCmtyWages, instruments=~WY + Missions, data=HiG)) #but P=1.0
HG <- tsls(HiGod4 ~ AnimXbwealth + FxCmtyWages, instruments=~WY + Logdate, data=HiG)
summary(tsls(HiGod4 ~ AnimXbwealth + FxCmtyWages, instruments=~WY + Logdate, data=HiG)) #N.Sig. p(F)=.44
HG <- tsls(HiGod4 ~ No_rain_Dry + SuperjhWriting, instruments=~WY + ww %*% Logdate, data=HiG)
summary(tsls(HiGod4 ~ No_rain_Dry + SuperjhWriting, instruments=~WY + ww %*% Logdate, data=HiG)) #N.Sig. p(S)=.68
###MORE VARIABLES THAN INSTRUMENTS: DOESNT WORK #non-invertible
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + Missions, instruments=~WY + ww %*% Logdate, data=HiG) #order 4
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + Missions, instruments=~WY + ww %*% Logdate, data=HiG)) #error
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + FxCmtyWages + Missions, instruments=~WY + Logdate, data=HiG) #
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + FxCmtyWages+ Missions, instruments=~WY + Logdate, data=HiG))
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + FxCmtyWages, instruments=~WY + Logdate, data=HiG) #
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + FxCmtyWages, instruments=~WY + Logdate, data=HiG))
HG <- tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + Missions, instruments=~WY, data=HiG) # DOESNT WORK non-invertible
summary(tsls(HiGod4 ~ AnimXbwealth + SuperjhWriting + Missions, instruments=~WY, data=HiG))
HG <- tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages+ SuperjhWriting, instruments=~WY + ww %*% Logdate, data=HiG)
summary(tsls(HiGod4 ~ No_rain_Dry + FxCmtyWages+ SuperjhWriting, instruments=~WY + ww %*% Logdate, data=HiG)) #non-inv
lm() Linear models THESE ARE *NOT* tsls
lm(HiGod4 ~ FxCmtyWages + WY)
summary(lm(HiGod4 ~ FxCmtyWages + WY))
lm(HiGod4 ~ FxCmtyWages + WY)
summary(lm(HiGod4 ~ FxCmtyWages + WY))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.52936 0.51767 1.023 0.309 FxCmtyWages 0.14298 0.09695 1.475 0.143 WY -0.41801 0.37430 -1.117 0.267
lm(HiGod4 ~ No_rain_Dry + AnimXbwealth + FxCmtyWages+ SuperjhWriting + WY)
summary(lm(HiGod4 ~ No_rain_Dry + AnimXbwealth + FxCmtyWages+ SuperjhWriting + WY))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.30687 0.44387 0.691 0.49086 No_rain_Dry 0.23509 0.09277 2.534 0.01274 * AnimXbwealth 0.24216 0.07961 3.042 0.00296 ** FxCmtyWages 0.07455 0.08989 0.829 0.40880 SuperjhWriting 0.38068 0.09003 4.228 5.01e-05 *** WY -0.56699 0.30080 -1.885 0.06217 .
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.3661 1.243e-15 -2.944e+14 0
Missions -0.8334 1.314e-15 -6.340e+14 0
summary(lm(HiGod4 ~ Missions + SuperjhWriting + WY))
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.34086 0.37783 0.902 0.368
Missions 0.29578 0.39176 0.755 0.451
SuperjhWriting 0.34265 0.07157 4.788 3.48e-06 ***
WY -0.25155 0.27147 -0.927 0.355
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.36465 0.35088 1.039 0.300
AnimXbwealth 0.35101 0.06457 5.436 1.74e-07 ***
SuperjhWriting 0.31173 0.06459 4.827 2.93e-06 ***
WY -0.26652 0.25217 -1.057 0.292
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.37712 0.35082 1.075 0.284 AnimXbwealth 0.37722 0.06868 5.493 1.32e-07 *** SuperjhWriting 0.32952 0.06648 4.956 1.64e-06 *** Missions -0.43178 0.38706 -1.116 0.266 WY -0.27209 0.25205 -1.079 0.282
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.004754 0.063990 0.074 0.941
AnimXbwealth 0.376124 0.068699 5.475 1.44e-07 ***
SuperjhWriting 0.330724 0.066503 4.973 1.52e-06 ***
Missions -0.423511 0.387163 -1.094 0.275
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.70788 0.42297 1.674 0.0972 .
AnimXbwealth 0.36085 0.08096 4.457 2.07e-05 ***
SuperjhWriting 0.46194 0.08933 5.171 1.10e-06 ***
Missions -0.74544 0.45999 -1.621 0.1081
FxCmtyWages 0.03265 0.08873 0.368 0.7136
WY -0.55730 0.30596 -1.821 0.0714
HG <- tsls(HiGod4 ~ SuperjhWriting, instruments=~ WY, data=HiG) # THIS WORKS but no second variable works!
summary(tsls(HiGod4 ~ SuperjhWriting, instruments=~WY, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.641e-16 0.2341 -3.264e-15 1.0000 SuperjhWriting 3.403e+00 11.0488 3.080e-01 0.7585
HG <- tsls(HiGod4 ~ AnimXbwealth, instruments=~ WY, data=HiG) # THIS WORKS but no second variable works!
summary(tsls(HiGod4 ~ AnimXbwealth, instruments=~WY, data=HiG))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.870e-16 0.9548 -4.053e-16 1.0000
AnimXbwealth -1.256e+01 166.3591 -7.551e-02 0.9399
HighGods (Moral)
Lavaan manual Page 35 Yves_Rosseel Lavaan setwd('/Users/drwhite/Documents/3sls/sccs/') #load("HiG.Rdata") load(url('http://intersci.ss.uci.edu/wiki/R/source/HiG.Rdata')) library(lavaan) model <- ' #first try # latent variable definitions Malth =~ FxCmtyWages + AnimXbwealth + SuperjhWriting # regressions HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY +Malth # residual correlations #Logdate ~~ HiGod4 #Logdate ~~ FxCmtyWages #Logdate ~~ AnimXbwealth #Missions ~~ HiGod4 #Missions ~~ FxCmtyWages '
fit <- sem(model, data=HiG) #lavaan WARNING: lhs and rhs are not BOTH in ov.names or lv.names. #lhs = Logdate #BLOCKED #rhs = HiGod4 #BLOCKED summary(fit, fit.measures=TRUE)
Lavaan (0.4-9) converged normally after 36 iterations
Used Total Number of observations 112 186
Estimator ML Minimum Function Chi-square 59.109 Degrees of freedom 7 P-value 0.000
Chi-square test baseline model:
Minimum Function Chi-square 121.337 Degrees of freedom 14 P-value 0.000
Full model versus baseline model:
Comparative Fit Index (CFI) 0.515 Tucker-Lewis Index (TLI) 0.029
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -780.349 Loglikelihood unrestricted model (H1) -750.794
Number of free parameters 11 Akaike (AIC) 1582.697 Bayesian (BIC) 1612.601 Sample-size adjusted Bayesian (BIC) 1577.837
Root Mean Square Error of Approximation:
RMSEA 0.258 90 Percent Confidence Interval 0.200 0.320 P-value RMSEA <= 0.05 0.000
Standardized Root Mean Square Residual:
SRMR 0.151
Parameter estimates:
Information Expected Standard Errors Standard
Estimate Std.err Z-value P(>|z|)
Latent variables:
Malth =~ FxCmtyWages 1.000 AnimXbwealth 0.530 0.308 1.721 0.085 SuperjhWritin -2.019 2.625 -0.769 0.442
Regressions:
HiGod4 ~ No_rain_Dry 0.308 0.082 3.753 0.000 FxCmtyWages 0.232 0.083 2.800 0.005 AnimXbwealth 0.272 0.076 3.555 0.000 WY -0.492 0.318 -1.546 0.122
Variances:
FxCmtyWages 1.179 0.302 3.897 0.000 AnimXbwealth 1.214 0.178 6.834 0.000 SuperjhWritin 1.751 1.079 1.623 0.105 HiGod4 0.754 0.101 7.483 0.000 Malth -0.188 0.244 -0.770 0.441 #instruments=~WY+Logdate+Missions
2nd model
The TLI is greater (improved by the Latent Variable 'Malth'), as incremental fit, bigger is better see: TLI versus CFI
setwd('/Users/drwhite/Documents/3sls/sccs/')
load(url('http://intersci.ss.uci.edu/wiki/R/source/HiG.Rdata'))
library(lavaan)
model <- '
# latent variable definitions
Malth =~ FxCmtyWages + AnimXbwealth + SuperjhWriting
# regressions
HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth #adding Mathusian event raises the TFI measure
# residual correlations
#Logdate ~~ HiGod4
#Logdate ~~ FxCmtyWages
#Logdate ~~ AnimXbwealth
#Missions ~~ HiGod4
#Missions ~~ FxCmtyWages
'
fit <- sem(model, data=HiG) summary(fit, fit.measures=TRUE) Lavaan (0.4-9) converged normally after 46 iterations
Used Total Number of observations 112 186
Estimator ML Minimum Function Chi-square 40.805 Degrees of freedom 6 P-value 0.000
Chi-square test baseline model:
Minimum Function Chi-square 121.337 Degrees of freedom 14 P-value 0.000
Full model versus baseline model:
Comparative Fit Index (CFI) 0.676 Tucker-Lewis Index (TLI) 0.243
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -771.197 Loglikelihood unrestricted model (H1) -750.794
Number of free parameters 12 Akaike (AIC) 1566.394 Bayesian (BIC) 1599.016 Sample-size adjusted Bayesian (BIC) 1561.091
Root Mean Square Error of Approximation:
RMSEA 0.228 90 Percent Confidence Interval 0.165 0.296 P-value RMSEA <= 0.05 0.000
Standardized Root Mean Square Residual:
SRMR 0.134
Parameter estimates:
Information Expected Standard Errors Standard
Estimate Std.err Z-value P(>|z|)
Latent variables:
Malth =~ FxCmtyWages 1.000 AnimXbwealth 0.530 0.308 1.721 0.085 SuperjhWritin -2.019 2.625 -0.769 0.442
Regressions:
HiGod4 ~ No_rain_Dry 0.209 0.076 2.767 0.006 FxCmtyWages 0.346 0.167 2.071 0.038 AnimXbwealth 0.382 0.105 3.638 0.000 WY -0.518 0.293 -1.767 0.077 Malth 0.651 0.284 2.294 0.022
Variances:
FxCmtyWages 1.179 0.302 3.897 0.000 AnimXbwealth 1.214 0.178 6.834 0.000 SuperjhWritin 1.751 1.079 1.623 0.105 HiGod4 0.861 0.176 4.904 0.000 Malth -0.188 0.244 -0.770 0.441
Bollen 1989
Page 35
library(lavaan) model <- ' # latent variable definitions ind60 =~ x1 + x2 + x3 dem60 =~ y1 + a*y2 + b*y3 + c*y4 dem65 =~ y5 + a*y6 + b*y7 + c*y8 # regressions dem60 ~ ind60 dem65 ~ ind60 + dem65 # residual correlations y1 ~~ y5 y2 ~~ y4 + y6 y3 ~~ y7 y4 ~~ y8 y6 ~~ y8'
fit <- sem(model, data=PoliticalDemocracy) summary(fit, fit.measures=TRUE)
Error in solve.default(E) : system is computationally singular: reciprocal condition number = 3.06231e-19
[lavaan message:] could not compute standard errors!
You can still request a summary of the fit to inspect the current estimates of the parameters.
Error discussion
data() in its most common usage just tells R to load a stored data object from an attached package. if you do help("data") you will see there is another application if package=NULL, so that it will look in a folder data of the current working directory. that might be simplest for you.
do not confuse it with the parameter name data= in functions -- that generally just allows you to specify a data frame that gives bindings for variables.
the solution for you is to have in your working .GlobalEnv a data frame of any name, say "foo", and then you can use data=foo as one of the parameter settings for your sem call
please do not ask me more questions about this. this is for R-help
#this works: lm(HiGod$HiGod4 ~ HiGod$No_rain_Dry + HiGod$AnimXbwealth + HiGod$SuperjhWriting + HiGod$Missions + HiGod$WY)
lm(HiG$HiGod4 ~ HiG$No_rain_Dry + HiG$AnimXbwealth + HiG$SuperjhWriting + HiG$Missions, instruments=HiG$WY) #works
SKIP THIS BELOW INSTRUMENT (AFTER IMPUTATION) IS lm(WY~WX) Instrument = Wdd %*% No_rain_Dry + Wdd %*% AnimXbwealth + Wdd %*% SuperjhWriting + Wdd %*% Missions
3rd model
For TLI, as incremental fit, bigger is better, see: TLI versus CFI
- Chi-square test baseline model P-value .000.
- Standardized Root Mean Square Residual (SRMR) .087: The SRMR is an absolute measure of fit and a value of zero indicates perfect fit. It is an absolute measure of fit and is defined as the standardized difference between the observed correlation and the predicted correlation. It is a positively biased measure and an absolute measure of fit. The bias is greater for small N and for low df studies. This measure tends to be smaller as sample size increases and as the number of parameters in the model increases. The SRMR has no penalty for model complexity. A value less than .08 is generally considered a good fit (Hu & Bentler, 1999). Li‐tze Hu & Peter M. Bentler. 1999. Cutoff criteria for fit indexes in covariance structure analysis: Conventional criteria versus new alternatives. Structural Equation Modeling: A Multidisciplinary Journal 6(1): 1-55.
- Root Mean Square Error of Approximation (RMSEA) .156, at the threshold of significance: P-value=.011. Make sure RMSEA is no smaller than 0.158 for the null model.
- Comparative Fit Index (CFI) .878++ (RMSEA=.156): The CFI should not be computed if the RMSEA of the null model is less than 0.158 or otherwise one will obtain too small a value of the CFI. If the CFI is less than one, then the CFI is always greater than the TLI. CFI pays a penalty of one for every parameter estimated.
- Tucker-Lewis Index (TLI) .587 on 6 variables: If you run the CFA on just the 5 indicators, you might have a nice TLI of .95. However, if you add in the 7 experimental variables, your TLI might sink below .90 because now the null model will not be so "bad" because you now have added to the model 7 variables who have zero correlations with each other.
Lavaan (0.4-9) converged normally after 46 iterations
setwd('/Users/drwhite/Documents/3sls/sccs/')
load(url('http://intersci.ss.uci.edu/wiki/R/source/HiG.Rdata'))
library(lavaan)
Program script (lavaan)
setwd('/Users/drwhite/Documents/3sls/sccs/')
load(url('http://intersci.ss.uci.edu/wiki/R/source/HiG.Rdata'))
library(lavaan)
model <- '
# latent variable definitions
Malth =~ FxCmtyWages + SuperjhWriting
# regressions
HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth
'
fit <- sem(model, data=HiG)
summary(fit, fit.measures=TRUE)
Results = Lavaan (0.4-9) converged normally after 22 iterations
HiGod4 Malth =~ FxCmtyWages + SuperjhWriting Malth
.288 .602 .616 .389 Latent variable variances
1.000 .975 Latent variable coefficients
HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth HiGod4= dependent variable
Coef .209 -.352 .234 -.518 1.427 HiGod4 regression coefficients
OLDER VERSION OBSOLETE Results = using only the 112 cases completely coded: SHOULD IMPROVE WITH MISSING DATA IMPUTATION (MI) HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth coef .209 -.352 .234 -.518 1.427 Malth =~ FxCmtyWages + SuperjhWriting # + AnimXbwealth HiGod4 variances .206 .616 .389 .288 model <- ' # latent variable definitions Malth =~ FxCmtyWages + SuperjhWriting # + AnimXbwealth # regressions HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth ' #adding Mathusian raises TFI fit <- sem(model, data=HiG) summary(fit, fit.measures=TRUE) Lavaan (0.4-9) converged normally after 22 iterations
Used Total Number of observations 112 186
Estimator ML Minimum Function Chi-square 18.665 Degrees of freedom 5 P-value 0.002
Chi-square test baseline model:
Minimum Function Chi-square 91.445 Degrees of freedom 12 P-value 0.000
Full model versus baseline model:
Comparative Fit Index (CFI) 0.828 Tucker-Lewis Index (TLI) 0.587
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -760.127 Loglikelihood unrestricted model (H1) -750.794
Number of free parameters 10 Akaike (AIC) 1540.254 Bayesian (BIC) 1567.439 Sample-size adjusted Bayesian (BIC) 1535.835
Root Mean Square Error of Approximation:
RMSEA 0.156 90 Percent Confidence Interval 0.085 0.235 P-value RMSEA <= 0.05 0.011
Standardized Root Mean Square Residual:
SRMR 0.087
Parameter estimates:
Information Expected Standard Errors Standard
Estimate Std.err Z-value P(>|z|)
Latent variables:
Malth =~ FxCmtyWages 1.000 SuperjhWritin 0.975
Regressions:
HiGod4 ~ No_rain_Dry 0.209 FxCmtyWages -0.352 AnimXbwealth 0.234 WY -0.518 Malth 1.427
Variances:
FxCmtyWages 0.602 SuperjhWritin 0.616 HiGod4 0.288 Malth 0.389
4th model
Wlink=matrix(data = 0, nrow = 186, ncol = 186)
for(i in 1:185) {
Wlink[i,i+1]=1
Wlink[i+1,i]=1}
for(i in 1:184) {
Wlink[i,i+2]=.7
Wlink[i+2,i]=.7}
for(i in 1:183) {
Wlink[i,i+3]=.3
Wlink[i+3,i]=.3}
#for(i in 1:182) {
#Wlink[i,i+4]=.2
#Wlink[i+4,i]=.2}
Wlink=Wlink/rowSums(Wlink)
#Try to copy and past this script into R and then check the first and last elements of the matrix:
Wlink[1:3,1:3]
Wlink[184:186,184:186]
sum(Wlink[55,1:186]) #try changing 55 to any other number between 1 and 186: rows sum to 1.
rowSums(Wlink)
Wdate=Wlink %*% HiG$Logdate
HiG$Wdata = Wdate
setwd('/Users/drwhite/Documents/3sls/sccs/')
load(url('http://intersci.ss.uci.edu/wiki/R/source/HiG.Rdata'))
library(lavaan)
Results = using only the 112 cases completely coded: SHOULD IMPROVE WITH MISSING DATA IMPUTATION (MI)
HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth
coef .209 -.352 .234 -.518 1.427
Malth =~ FxCmtyWages + SuperjhWriting # + AnimXbwealth HiGod4
variances .206 .616 .389 .288
model <- '
# latent variable definitions
Malth =~ FxCmtyWages + SuperjhWriting + Wdata # + AnimXbwealth
# regressions
HiGod4 ~ No_rain_Dry+FxCmtyWages+AnimXbwealth + WY + Malth
' #adding Mathusian raises TFI
fit <- sem(model, data=HiG)
summary(fit, fit.measures=TRUE)
Lavaan (0.4-9) converged normally after 31 iterations
Used Total Number of observations 112 186
Estimator ML Minimum Function Chi-square 22.699 Degrees of freedom 10 P-value 0.012
Chi-square test baseline model:
Minimum Function Chi-square 111.642 Degrees of freedom 18 P-value 0.000
Full model versus baseline model:
Comparative Fit Index (CFI) 0.864 Tucker-Lewis Index (TLI) 0.756
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -821.153 Loglikelihood unrestricted model (H1) -809.803
Number of free parameters 12 Akaike (AIC) 1666.306 Bayesian (BIC) 1698.927 Sample-size adjusted Bayesian (BIC) 1661.003
Root Mean Square Error of Approximation:
RMSEA 0.106 90 Percent Confidence Interval 0.048 0.165 P-value RMSEA <= 0.05 0.056
Standardized Root Mean Square Residual:
SRMR 0.084
Parameter estimates:
Information Expected Standard Errors Standard
Estimate Std.err Z-value P(>|z|)
Latent variables:
Malth =~ FxCmtyWages 1.000 SuperjhWritin 1.244 0.416 2.989 0.003 Wdata -0.358 0.126 -2.849 0.004
Regressions:
HiGod4 ~ No_rain_Dry 0.186 0.084 2.210 0.027 FxCmtyWages -0.189 0.174 -1.087 0.277 AnimXbwealth 0.244 0.079 3.099 0.002 WY -0.445 0.288 -1.542 0.123 Malth 1.316 0.323 4.073 0.000
Variances:
FxCmtyWages 0.694 0.141 4.908 0.000 SuperjhWritin 0.526 0.143 3.688 0.000 Wdata 0.163 0.024 6.772 0.000 HiGod4 0.413 0.155 2.659 0.008 Malth 0.297 0.143 2.084 0.037
