# --- Leitura dos dados em .csv --- #

dados <- read.csv("http://www.icmc.usp.br/~cibele/Dados/MBA_CAR_ATTRIB.csv", sep=";")

names(dados)

attach(dados)

# --- Matriz de dados --- #

matrizdados<-na.omit(dados[,3:18])

n <- nrow(matrizdados)

p <- ncol(matrizdados)


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

Xbarra <- mean(matrizdados)

S <- cov(matrizdados)

par(mfrow=c(2,8))
for(i in 1:16)
boxplot(matrizdados[,i], main=names(matrizdados)[i])

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

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

analisefatorial <- factanal(matrizdados, factors=3, 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]


# --- 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)