cfa-example.rnw

Reproducible example
====================

This example goes with a question posted to the lavaan google group. This code depends on a dataframe called mydata. It is a multivariate ordinal dataset simulated from the the MultiOrd package. The simulation takes a long time to run, so this file just brings in the 9 variable dataframe (n=500) dataframe from a gist. The gist is basically a dput of the dataframe that is too long for this document.

Full R code for this example [here](https://gist.github.com/ericpgreen/7092883).

First, load packages.

```{r setup, message=FALSE, highlight=TRUE}
  library(devtools) # to run script from gist and bring in data
  library(lavaan)
  library(semPlot)
  library(knitr)
  library(ggplot2)
  library(reshape2)
  library(psych)
```

Then read in the data.

```{r import, message=FALSE, highlight=TRUE}
# read in data from gist
# gist is at https://gist.github.com/ericpgreen/7091485
# this creates data frame mydata
  gist <- "https://gist.github.com/ericpgreen/7091485/raw/f4daec526bd69557874035b3c175b39cf6395408/simord.R"
  source_url(gist, sha1="da165a61f147592e6a25cf2f0dcaa85027605290")
  # if you get an error "SHA-1 hash of downloaded file does not match expected value"
  # it means the gist has been modified
  # run source_url(gist) to get new sha1 and replace string in quotes
# list of variables
  items <- c("v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8", "v9")
# head
  head(mydata)
```

Run some descriptives and plot histograms.

```{r des, message=FALSE, highlight=TRUE}
# alpha
  alpha(mydata)
# cor
  cor(mydata)
# histograms
  mydata.melt <- melt(mydata)
  ggplot(mydata.melt, aes(x = value)) + 
      facet_wrap(~variable,scales = "free_x") +
      geom_density() +
      xlim(1, 4)
```

Model A
-------
Model A is one first-order factor.

```{r a, message=FALSE, highlight=TRUE}
  mod.A <- '
  f0 =~ v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8 + v9
  '
  fit.mod.A <- cfa(mod.A, data = mydata, ordered=items)
  semPaths(fit.mod.A, "std", title = FALSE, 
           nCharNodes=0, edge.label.cex=0.6)
```

Model B
-------
Model B is 3 first order correlated factors (f1, f2, f3).

```{r b, message=FALSE, highlight=TRUE}
  mod.B <- '
  f1 =~ v1 + v2 + v3
  f2 =~ v4 + v5 + v6
  f3 =~ v7 + v8 + v9
  '
  fit.mod.B <- cfa(mod.B, data = mydata, ordered=items)
  semPaths(fit.mod.B, "std", title = FALSE,
           nCharNodes=0, edge.label.cex=0.6)
```

Model C
-------
Model C is 3 first order correlated factors (f1, f2, f3) with two cross-loading items on f1 and f2.

```{r c, message=FALSE, highlight=TRUE}
  mod.C <- '
  f1 =~ v1 + v2 + v3 + v4
  f2 =~ v4 + v5 + v6 + v7
  f3 =~ v7 + v8 + v9
  '
  fit.mod.C <- cfa(mod.C, data = mydata, ordered=items)
  semPaths(fit.mod.C, "std", title = FALSE,
           nCharNodes=0, edge.label.cex=0.6)
```

Model D
-------
Model D is a hierarchical model with 1 second-order factor (g) and 3 first-order uncorrelated factors (f1, f2, f3).

```{r d, message=FALSE, highlight=TRUE}
  mod.D <- '
  g =~ f1 + f2 + f3
  f1 =~ v1 + v2 + v3
  f2 =~ v4 + v5 + v6
  f3 =~ v7 + v8 + v9
  f1 ~~ 0*f2 + 0*f3
  f2 ~~ 0*f3
  '
  fit.mod.D <- cfa(mod.D, data = mydata, ordered=items, std.lv=TRUE)
  semPaths(fit.mod.D, "std", title = FALSE,
           nCharNodes=0, edge.label.cex=0.6)
```

Model E
-------
Model E is a bifactor model with 1 general factor (g) and 3 first-order uncorrelated factors (f1, f2, f3).

```{r e, message=FALSE, highlight=TRUE}
  mod.E <- '
  g =~ v1 + v2 + v3 + v4 + v5 + v6 + v7 + v8 + v9
  f1 =~ v1 + v2 + v3
  f2 =~ v4 + v5 + v6
  f3 =~ v7 + v8 + v9
  '
  fit.mod.E <- cfa(mod.E, data = mydata, ordered=items, orthogonal=TRUE, std.lv=TRUE)
  semPaths(fit.mod.E, "std", title = FALSE,
           nCharNodes=0, edge.label.cex=0.6, rotate=2)
  # rotate does not seem to be working
  # know how to flip so g is on left and f1, f2, f3 on right?
```