{"id":737,"date":"2013-07-31T10:49:20","date_gmt":"2013-07-31T10:49:20","guid":{"rendered":"http:\/\/www.georg-hosoya.de\/wordpress\/?p=737"},"modified":"2014-02-18T15:04:37","modified_gmt":"2014-02-18T15:04:37","slug":"rotational-invariance-in-factor-models-ii","status":"publish","type":"post","link":"https:\/\/www.georg-hosoya.de\/wordpress\/?p=737","title":{"rendered":"Rotational Invariance in factor models II"},"content":{"rendered":"<p>I have been working a bit more on estimating simple factor models with <a href=\"http:\/\/mcmc-jags.sourceforge.net\/\">JAGS<\/a>. I posted the problem from the last post over there at the <a href=\"http:\/\/sourceforge.net\/p\/mcmc-jags\/discussion\/\">open discussion forum<\/a> and had some very helpful input from Vincent Arel-Bundock from Michigan University and one of his co-authors, Walter R. Mebane which ultimately led to the code of the current post.<\/p>\n<p>Here I wanted to know how to estimate a simple factor model for continuous manifest variables and three latent, intercorrelated factors. The code is simple and straight forward. I put a scaled inverse Wishart prior on the variance-covariance matrix of factor loadings, following the approach in Gelman &amp; Hill (2007, p.376). For solving the invariance problem I set the loadings for one item on each factor to 1. In addition, the means of the latent trait distributions (factor values) were all constrained to zero. This is a very traditional approach, that is also implemented in e.g. <a href=\"http:\/\/www.statmodel.com\/\">MPlus<\/a>, the<a href=\"http:\/\/www.stata.com\/\"> STATA<\/a> program <a href=\"http:\/\/www.gllamm.org\/\">GLLAMM<\/a> or the <a href=\"http:\/\/www.r-project.org\/\">R<\/a> package <a href=\"http:\/\/lavaan.ugent.be\/\">lavaan<\/a>.<\/p>\n<p>In the following code I use the Holzinger &amp; Swineford dataset from the lavaan tutorial for comparison. First I estimate the model with lavaan and in a second step with JAGS. Both approaches show very comparable results, but keep in mind that this is somewhat experimental. When comparing coefficients please note that I coded the model in such a way that JAGS returns estimated standard deviations and not variances for the manifest items and latent variables.<\/p>\n<p>Ok, enough talk, here is the code. First the R-script:<\/p>\n<pre lang=\"rsplus\" line=\"1\">\r\n## holzinger.R\r\n# Experimental!\r\n\r\nlibrary(lavaan)\r\n\r\n# From the lavaan tutorial\r\nHS.model <- ' visual  =~ x1 + x2 + x3      \r\n              textual =~ x4 + x5 + x6\r\n              speed   =~ x7 + x8 + x9 '\r\n\r\nfit <- cfa(HS.model, data=HolzingerSwineford1939, meanstructure=TRUE)\r\n\r\nsummary(fit, standardized=TRUE)\r\n\r\nlibrary(R2jags)\r\ndata<-HolzingerSwineford1939\r\n## Transform data to long format----\r\n#  the complicated way\r\n\r\n# Dependent variable\r\ny<-c(data$x1, data$x2, data$x3, data$x4, data$x5, data$x6, data$x7, data$x8, data$x9)\r\n\r\n# Item identifier\r\nitem<-c(rep(1, 301), rep(2, 301), rep(3, 301), rep(4, 301), rep(5, 301), rep(6, 301), rep(7, 301), rep(8, 301), rep(9, 301))\r\n\r\n# Construct identifier (1: visual, 2: textual, 3: speed)\r\nconst<-c(rep(1,3*301 ), rep(2,3*301 ), rep(3,3*301 ))\r\n\r\n# Person id\r\nid<-rep(seq(1:301),9)\r\n## ---------------------------------\r\n# Fitting the model with JAGS\r\nW=(diag(3))\r\n\r\ndat<-list(\"y\", \"item\", \"const\", \"id\", \"W\")\r\nparameters<-(c(\"beta\", \"sigma.y\", \"lam\", \"sigma.eta1\", \"sigma.eta2\", \"sigma.eta3\", \"rho_eta13\", \"rho_eta12\", \"rho_eta23\", \"eta\"))\r\n\r\noutput<-jags.parallel(dat, inits=NULL, parameters,  model.file=\"3fa.txt\",n.chains=2, n.iter=3000, n.burnin=2000)\r\nplot(output)\r\n<\/pre>\n<p>Next the model code for JAGS:<\/p>\n<pre lang=\"rsplus\" line=\"1\">\r\n## 3fa.txt\r\n# Experimental!\r\nmodel{\r\n\r\n  # Likelihood\r\n  for(i in 1:2709) # Loop over observations\r\n  {\r\n      # Linear decomposition could probably done more\r\n      # elegantly with matrix algebra\r\n      mu[i]<-beta[item[i]] +\r\n         lam[const[i], item[i]]*eta[const[i], id[i]]+\r\n         lam[const[i], item[i]]*eta[const[i], id[i]]+\r\n         lam[const[i], item[i]]*eta[const[i], id[i]]\r\n\r\n      # Distributional assumptions for the \r\n      # observed variables\r\n         y[i]~dnorm(mu[i], tau.y[item[i]])\r\n  }\r\n\r\n  # Priors for the items\r\n    for(i in 1:9)\r\n    {\r\n      beta[i]~dnorm(0,1.0E-3) \r\n      tau.y[i] <- pow(sigma.y[i], -2)  \r\n      sigma.y[i] ~ dunif (0, 100)   \r\n    }\r\n\r\n    # Matrix of loadings [construct, item]\r\n    # Constraints and priors \r\n     lam[1,1]<-1\r\n     lam[2,1]<-0\r\n     lam[3,1]<-0\r\n     lam[1,2]~dnorm(0,1.0E-3) \r\n     lam[2,2]<-0\r\n     lam[3,2]<-0\r\n     lam[1,3]~dnorm(0,1.0E-3) \r\n     lam[2,3]<-0\r\n     lam[3,3]<-0\r\n     lam[1,4]<-0\r\n     lam[2,4]<-1\r\n     lam[3,4]<-0\r\n     lam[1,5]<-0\r\n     lam[2,5]~dnorm(0,1.0E-3)\r\n     lam[3,5]<-0\r\n     lam[1,6]<-0\r\n     lam[2,6]~dnorm(0,1.0E-3)\r\n     lam[3,6]<-0\r\n     lam[1,7]<-0\r\n     lam[2,7]<-0\r\n     lam[3,7]<-1\r\n     lam[1,8]<-0\r\n     lam[2,8]<-0\r\n     lam[3,8]~dnorm(0,1.0E-3)\r\n     lam[1,9]<-0\r\n     lam[2,9]<-0\r\n     lam[3,9]~dnorm(0,1.0E-3)\r\n\r\n    # Scaled inverse Wishart prior for\r\n    # variance-covariance matrix of factor values\r\n    # after Gelman and Hill (2007)\r\n\r\n    for(i in 1:301)\r\n    {  \r\n      eta[1,i]<-xi.eta1*eta.raw[i,1]\r\n      eta[2,i]<-xi.eta2*eta.raw[i,2]\r\n      eta[3,i]<-xi.eta3*eta.raw[i,3]\r\n      eta.raw[i,1:3]~dmnorm(eta.raw.hat[i,],tau.eta.raw[,])\r\n      eta.raw.hat[i,1]<-mu.eta1.raw\r\n      eta.raw.hat[i,2]<-mu.eta2.raw\r\n      eta.raw.hat[i,3]<-mu.eta3.raw\r\n    }\r\n      mu.eta1<-xi.eta1*mu.eta1.raw\r\n      mu.eta2<-xi.eta2*mu.eta2.raw\r\n      mu.eta3<-xi.eta3*mu.eta3.raw\r\n      mu.eta1.raw<-0\r\n      mu.eta2.raw<-0\r\n      mu.eta3.raw<-0\r\n      xi.eta1~dunif(0,100)\r\n      xi.eta2~dunif(0,100)\r\n      xi.eta3~dunif(0,100)\r\n      tau.eta.raw[1:3,1:3] ~ dwish(W[,],df)\r\n      df<-4\r\n      sigma.eta.raw[1:3,1:3]<-inverse(tau.eta.raw[,])\r\n      sigma.eta1<-xi.eta1*sqrt(sigma.eta.raw[1,1])\r\n      sigma.eta2<-xi.eta2*sqrt(sigma.eta.raw[2,2])\r\n      sigma.eta3<-xi.eta3*sqrt(sigma.eta.raw[3,3])\r\n      rho_eta12<-sigma.eta.raw[1,2]\/sqrt(sigma.eta.raw[1,1]*sigma.eta.raw[2,2])\r\n      rho_eta13<-sigma.eta.raw[1,3]\/sqrt(sigma.eta.raw[1,1]*sigma.eta.raw[3,3])\r\n      rho_eta23<-sigma.eta.raw[2,3]\/sqrt(sigma.eta.raw[2,2]*sigma.eta.raw[3,3])\r\n}\r\n<\/pre>\n","protected":false},"excerpt":{"rendered":"<p>I have been working a bit more on estimating simple factor models with JAGS. I posted the problem from the last post over there at the open discussion forum and had some very helpful input from Vincent Arel-Bundock from Michigan &hellip; <a href=\"https:\/\/www.georg-hosoya.de\/wordpress\/?p=737\">Continue reading <span class=\"meta-nav\">&rarr;<\/span><\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts\/737"}],"collection":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=737"}],"version-history":[{"count":19,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts\/737\/revisions"}],"predecessor-version":[{"id":840,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=\/wp\/v2\/posts\/737\/revisions\/840"}],"wp:attachment":[{"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=737"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=737"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.georg-hosoya.de\/wordpress\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=737"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}