# --- Leitura dos matrizdados em .csv --- #

dados <- read.table("http://wiki.icmc.usp.br/images/9/95/BDI.txt")

names(dados)

attach(dados)

matrizdados<-dados[,1:21]


# --- Análise descritiva e exploratória --- #

Xbarra <- mean(matrizdados)

S <- cov(matrizdados)

n <- nrow(matrizdados)

p <- ncol(matrizdados)


par(mfrow=c(4,6))
for(i in 1:21)
boxplot(matrizdados[,i], main=names(matrizdados)[i])


pairs(matrizdados[,1:7])
pairs(matrizdados[,8:14])
pairs(matrizdados[,15:21])

# --- Estimação via componentes principais --- #
# --- Critério: Autovalores maiores do que um --- #

sum(princomp(matrizdados, cor=TRUE)[[1]]>1)

analisefatorial <- factanal(matrizdados, factors=4, rotation="varimax")
analisefatorial

#--- Cargas fatoriais ---#
cargas <- analisefatorial$loadings[]


#--- Comunalidades  ---#
comunalidades <- diag(cargas %*% t(cargas))

#--- Variância especifica ---#
varespecifica <- analisefatorial$uniqueness

#--- Estimativa de R  ---#
Rchapeu <- cargas %*% t(cargas) + varespecifica

#--- Matriz de resíduos  ---#
resíduos <- cor(matrizdados) - Rchapeu



# --- Exercício: como ficaria a análise sem rotação? --- #
# --- Cargas fatoriais --- #

sqrt(eigen(cor(matrizdados))[[1]][1]) * eigen(cor(matrizdados))[[2]][,1]
sqrt(eigen(cor(matrizdados))[[1]][2]) * eigen(cor(matrizdados))[[2]][,2]
sqrt(eigen(cor(matrizdados))[[1]][3]) * eigen(cor(matrizdados))[[2]][,3]
sqrt(eigen(cor(matrizdados))[[1]][4]) * eigen(cor(matrizdados))[[2]][,4]

# --- Critério: Scree plot --- #

 plot(eigen(cor(matrizdados))[[1]], type="b", pch=16, main="Scree plot", ylab="autovalor", xlab="ordem da componente")


# --- Exercício: Os dados tem distribuição normal? --- #
# --- Posso estimar o modelo utilizando o método da máxima verossimilhança?--- #

#install.packages("mvShapiroTest") # Executar se não tiver o pacote

library(mvShapiroTest)

mvShapiro.Test(as.matrix(matrizdados))

# Distâncias dj

par(mfrow=c(1,1))

d<-rep(0,n)

for (i in 1: n)
  d[i] <- as.matrix(matrizdados[i,] - Xbarra) %*% solve(S) %*% as.matrix(t(matrizdados[i,] - Xbarra))

qqplot(d, qchisq(ppoints(n), p), pch=16)
abline(0,1)

# --- Envelope para as distâncias - executar código da função antes --- #

source("http://wiki.icmc.usp.br/images/2/23/Envelope.distancias.txt")
envelope.distancias(d, n, Xbarra, S)




#Exercício: Fazer análises separando por gênero