# Back to the US in 1977
# Do a one-factor analysis and compute factor scores
state.fa1 <- factanal(state.x77,factors=1,scores="regression")
# Have R tell us about it
state.fa1
# Compare the factor loadings here to the first principal component we got
# last time
# Make a map again
# Function is a repeat from last time
# Plot the state abbrevations in position, with scaled sizes
# Linearly scale the sizes from the given minimum to the maximum
# Inputs: vector of raw numbers, minimum size for plot,
# maximum size
# Outputs: Rescaled sizes (invisible)
plot.states_scaled <- function(sizes,min.size=0.4,max.size=2,...) {
out.range = max.size - min.size
in.range = max(sizes)-min(sizes)
scaled.sizes = out.range*((sizes-min(sizes))/in.range)
sizes = scaled.sizes + min.size
plot(state.center,type="n",...)
text(state.center,state.abb,cex=sizes)
invisible(sizes)
}
# Actual map-making
plot.states_scaled(state.fa1$score[,1],min.size=0.3,max.size=1.5,
xlab="longitude",ylab="latitude")
# 3D visualization
# You need the "scatterplot3d" library from CRAN for this
require(scatterplot3d)
# Make a matrix with the x,y,z values (z=factor scores)
state.xyz <- cbind(state.center$x,state.center$y,
state.fa1$scores[,1])
colnames(state.xyz)=c("x","y","z")
# Plot it!
state.3d <- scatterplot3d(state.xyz,type="h",
xlab="longitude",
ylab="latitude",
zlab="factor score",
cex.symbol=0.01,color="grey")
# Add labels
text(state.3d$xyz.convert(state.xyz),state.abb)
# How good is the factor model, though?
# Try models with 1--4 factors and get the p-value for each
pvalues <- sapply(1:4,function(q){factanal(state.x77,factors=q)$PVAL})
# What are these p-values?
signif(pvalues,2)
# Plot vs. the nominal 5% level
plot(1:4,pvalues,xlab="q (number of factors)", ylab="pvalue",
log="y",ylim=c(1e-11,0.04))
abline(h=0.05,lty=2)
# Another example: sleep in mammals
sleep <- read.csv("sleep.txt")
# The fifth column is the sum of the 3rd and 4th columns, which makes everyone
# very unhappy!
sleep <- sleep[,-5]
# There are missing observations, so we need to calculate the covariance
# matrix outside factanal
sleep.cov <- cov(sleep,use="pairwise.complete.obs")
# Try one factor
sleep.fa1 <- factanal(sleep,factors=1,covmat=sleep.cov)
sleep.fa1
# Try to tell a story about this factor based on the loadings
# Notice that it won't give us a p-value because of the missing observations
# Try two factors
sleep.fa2 <- factanal(sleep,factors=2,covmat=sleep.cov)
sleep.fa2
# Notice that the first factor has changed
# Try to tell a story about these two factors