## Example from Venables and Ripley (2002, p. 323)
## Previously from Bartholomew and Knott (1999, p. 68--72)
## Originally from Smith and Stanley (1983)
## Replicated from example(ability.cov)
man <- make_manifest(covmat = ability.cov)
## Not run:
# ## Here is the easy way to set up a SEFA model, which uses pop-up menus
# res <- make_restrictions(manifest = man, factors = 2, model = "SEFA")
# ## End(Not run)
## This is the hard way to set up a restrictions object without pop-up menus
beta <- matrix(NA_real_, nrow = nrow(cormat(man)), ncol = 2)
rownames(beta) <- rownames(cormat(man))
free <- is.na(beta)
beta <- new("parameter.coef.SEFA", x = beta, free = free, num_free = sum(free))
Phi <- diag(2)
free <- lower.tri(Phi)
Phi <- new("parameter.cormat", x = Phi, free = free, num_free = sum(free))
res <- make_restrictions(manifest = man, beta = beta, Phi = Phi)
# This is how to make starting values where Phi is the correlation matrix
# among factors, beta is the matrix of coefficients, and the scales are
# the logarithm of the sample standard deviations. It is also the MLE.
starts <- c( 4.46294498156615e-01, # Phi_{21}
4.67036349420035e-01, # beta_{11}
6.42220238211291e-01, # beta_{21}
8.88564379236454e-01, # beta_{31}
4.77779639176941e-01, # beta_{41}
-7.13405536379741e-02, # beta_{51}
-9.47782525342137e-08, # beta_{61}
4.04993872375487e-01, # beta_{12}
-1.04604290549591e-08, # beta_{22}
-9.44950629176182e-03, # beta_{32}
2.63078925240678e-04, # beta_{42}
9.38038168787216e-01, # beta_{52}
8.43618801925473e-01, # beta_{62}
log(man@sds)) # log manifest standard deviations
sefa <- Factanal(manifest = man, restrictions = res,
# NOTE: Do NOT specify any of the following tiny values in a
# real research situation; it is done here solely for speed
starting.values = starts, pop.size = 2, max.generations = 6,
wait.generations = 1)
nsim <- 101 # number of simulations, also too small for real work
show(sefa)
summary(sefa, nsim = nsim)
model_comparison(sefa, nsim = nsim)
stuff <- list() # output list for various methods
stuff$model.matrix <- model.matrix(sefa) # sample correlation matrix
stuff$fitted <- fitted(sefa, reduced = TRUE) # reduced covariance matrix
stuff$residuals <- residuals(sefa) # difference between model.matrix and fitted
stuff$rstandard <- rstandard(sefa) # normalized residual matrix
stuff$weights <- weights(sefa) # (scaled) approximate weights for residuals
stuff$influence <- influence(sefa) # weights * residuals
stuff$cormat <- cormat(sefa, matrix = "RF") # reference factor correlations
stuff$uniquenesses <- uniquenesses(sefa, standardized = FALSE) # uniquenesses
stuff$FC <- loadings(sefa, matrix = "FC") # factor contribution matrix
stuff$draws <- FA2draws(sefa, nsim = nsim) # draws from sampling distribution
if(require(nFactors)) screeplot(sefa) # Enhanced scree plot
profile(sefa) # profile plots of non-free parameters
pairs(sefa) # Thurstone-style plot
if(require(Rgraphviz)) plot(sefa) # DAG
Run the code above in your browser using DataLab