# HW 4 solution R code.

# This is merely the R code;
# Actual HWs should also include clearly stated answers 
# and insightful comments!


################### #1


school.sub.corr <- matrix( c(
1,.44,.41,.29,.33,.25,
.44,1,.35,.35,.32,.33,
.41,.35,1,.16,.19,.18,
.29,.35,.16,1,.59,.47,
.33,.32,.19,.59,1,.46,
.25,.33,.18,.47,.46,1
), nrow=6, ncol=6, byrow=T)


factanal(covmat=school.sub.corr, factors=2, rotation='none', n.obs=220)

factanal(covmat=school.sub.corr, factors=2, rotation='varimax', n.obs=220)

# Communalities:

1-factanal(covmat=school.sub.corr, factors=2, rotation='varimax', n.obs=220)$uniquenesses

# Specific variances:

factanal(covmat=school.sub.corr, factors=2, rotation='varimax', n.obs=220)$uniquenesses


################### #2

# problem 4.7:

pain.corr <- matrix( c(
1,-.04,.61,.45,.03,-.29,-.3,.45,.3,
-.04,1,-.07,-.12,.49,.43,.3,-.31,-.17,
.61,-.07,1,.59,.03,-.13,-.24,.59,.32,
.45,-.12,.59,1,-.08,-.21,-.19,.63,.37,
.03,.49,.03,-.08,1,.47,.41,-.14,-.24,
-.29,.43,-.13,-.21,.47,1,.63,-.13,-.15,
-.3,.3,-.24,-.19,.41,.63,1,-.26,-.29,
.45,-.31,.59,.63,-.14,-.13,-.26,1,.4,
.3,-.17,.32,.37,-.24,-.15,-.29,.4,1
), nrow=9, ncol=9, byrow=T)

factanal(covmat=pain.corr, factors=2, rotation='none', n.obs=123)
factanal(covmat=pain.corr, factors=3, rotation='none', n.obs=123)
factanal(covmat=pain.corr, factors=4, rotation='none', n.obs=123)

# for 4.7(c)

factanal(covmat=pain.corr, factors=3, rotation='varimax', n.obs=123)

#################  #3

food.full <- read.table("http://www.stat.sc.edu/~hitchcock/foodstuffs.txt", header=T)
food.labels <- as.character(food.full[,1])
food.data <- food.full[,-1]


factanal(food.data,factors=1, rotation="varimax")
# One factor is clearly not enough (tiny P-value).

factanal(food.data,factors=2, rotation="varimax")
# Two factors are not quite enough based on the P-value (but it is kind of close).

factanal(food.data,factors=3, rotation="varimax")
# Three factors are too many with only 5 variables.



food.fa.2.v <- factanal(food.data,factors=2, rotation="varimax")

# The communalities:

1-food.fa.2.v$uniquenesses


# Estimating factor scores for the life expectancy data set:

food.fa.2.scores <- factanal(food.data,factors=2, rotation="varimax",scores="regression")$scores


# Creating a data frame with the original data and the scores:

food.df <- data.frame(food.data, food.fa.2.scores)
attach(food.df)



# Looking back at the factor scores:

food.df[,c("Factor1","Factor2")]

cbind(food.labels, food.df[,c("Factor1","Factor2")])

food.labels[order(Factor1)]

food.labels[order(Factor2)]


### Doing 2-D scatterplot of the factor scores:

# Factor 2 vs. Factor 1:

plot(Factor1, Factor2, type='n', xlab='Factor 1 scores', ylab='Factor 2 scores')
text(Factor1, Factor2, labels = food.labels, cex = 0.7 )


#########################################################################################################