#----------------------------------------------------------------------------------
#                 PRÁTICA 2 - ANÁLISE AMOSTRA DE CRIANÇAS
#----------------------------------------------------------------------------------

# Altere o diretório de trabalho utilizando a função setwd("")

  setwd("C:/A DATA/PSP5103/2021/Praticas")


# Trabalhando com vetores e atribuições 
#----------------------------------------------------------------------------------

# Calculando algumas medidas descritivas para a amostra de crianças
#----------------------------------------------------------------------------------


peso <- c(64, 71, 53, 67, 55, 58, 77, 57, 56, 51, 76, 68)
peso

idade <- c(8, 10, 6, 11, 8, 7, 10, 9, 10, 6, 12, 9)
idade

# Calculando a média

length(peso)
sum(peso) 

peso.m <- sum(peso) / length(peso)   
peso.m

idade.m <- sum(idade) / length(idade)   
idade.m

# Calculando a variância

peso.var <- sum((peso-peso.m)^2) / (length(peso)-1)
peso.var

idade.var <- sum((idade-idade.m)^2) / (length(idade)-1)
idade.var


# Calculando o desvio padrão

peso.dp <- sqrt(peso.var)
peso.dp

idade.dp <- sqrt(idade.var)
idade.dp

# Calculando o coeficiente de correlação entre peso e idade

peso.pad <- (peso-peso.m)/peso.dp
peso.pad

idade.pad <- (idade-idade.m)/idade.dp
idade.pad 

r.peso.idade <- sum(peso.pad*idade.pad)/(length(idade)-1)
r.peso.idade


# Carregando o banco de dados
#----------------------------------------------------------------------------------

# Importe o banco de dados "Amostra de Crianças" e dê o nome de "dados" 
# Utilize o Menu
# File -> Import Dataset ->  from Excel

library(readxl)
dados <- read_excel("C:/A DATA/PSP5103/2021/Dados/Amostra de criancas.xlsx")
View(dados)

names(dados)    # mostra o nome das variáveis no banco de dados 
dim(dados)      # mostra as dimensões (No. de linhas e colunas) do bd
head(dados)
dados


#----------------------------------------------------------------------------------
#                            ANÁLISE DESCRITIVA
#----------------------------------------------------------------------------------


# Produzindo um resumo rápido do conteúdo do banco de dados...

summary(dados)

# ou de uma variável...

summary(dados$peso)


# DESCREVENDO UMA VARIÁVEL DE CADA VEZ
#----------------------------------------------------------------------------------

# Medidas descritivas
#---------------------

# Média, mediana, variância, desvio padrão e quantis

mean(dados$peso)      # média
var(dados$peso)       # variância
sd(dados$peso)        # desvio padrão
quantile(dados$peso,  probs = c(0, 0.25, 0.5, 0.75, 1))  #quartis
quantile(dados$peso,  probs = c(0, 0.333, 0.666, 1))     #tercis


# Medidas descritivas, incluindo o intervalo de confiança

library(epiDisplay)
ci.numeric(dados$peso, alpha = 0.05)   
ci.numeric(dados$idade, alpha = 0.05)   
ci.numeric(dados$altura, alpha = 0.05)   

# Boxplot
#---------

boxplot(dados$peso)
boxplot(dados$peso, ylab="Peso (libras)", col="orange")

# Veja as cores disponíveis no R:

colors()

par(mfrow=c(1,3))

boxplot(dados$peso, ylab="Peso (libras)", col="orange")
boxplot(dados$idade, ylab="Idade (anos)", col="orange")
boxplot(dados$altura, ylab="Altura (pés)", col="orange")

par(mfrow=c(1,1))


# Histograma
#------------

hist(dados$peso)

hist(dados$peso, freq = FALSE)
hist(dados$peso, breaks = 7)
hist(dados$peso, breaks = c(50, 52.5, 55, 57.5, 60, 62.5, 65, 67.5, 70, 72.5, 75, 77.5, 80))
hist(dados$peso, col="lightblue")


hist(dados$peso, freq = FALSE, col="lightblue", xlab="Peso (libras)", ylab="Densidade de frequência", main = "")
hist(dados$idade, freq = FALSE,  col="lightblue", xlab="Idade (anos)", ylab="Densidade de frequência", main = "")
hist(dados$altura, freq = FALSE,  col="lightblue", xlab="Altura (pés)", ylab="Densidade de frequência", main = "")


# DESCREVENDO A RELAÇÃO/ASSOCIAÇÃO ENTRE PESO E IDADE
#----------------------------------------------------------------------------------

# Coeficiente de correlação linear de Pearson
#----------------------------------------------

cor(dados$peso, dados$idade)
cor.test(dados$peso, dados$idade, alternative="two.sided", method="pearson")

cor(dados$peso, dados$altura)
cor.test(dados$peso, dados$altura, alternative="two.sided", method="pearson")


#Diagrama de dispersão
#----------------------

par(mfrow=c(1,2))
plot(peso~idade, xlab="Idade (anos)", ylab="Peso (libras)", col="blue", data=dados)
plot(peso~altura, xlab="Altura (pés)", ylab="Peso (libras)", col="blue", data=dados)

par(mfrow=c(1,1))


#----------------------------------------------------------------------------------
#                  MODELO DE REGRESSÃO LINEAR SIMPLES
#----------------------------------------------------------------------------------

# Checando as suposições do modelo
#-----------------------------------

# Gráficos para avaliar normalidade 

hist(dados$peso, freq = FALSE, col="lightblue", xlab="Peso (libras)", 
     ylab="Densidade de frequencia", main="")           

library(car)
qqPlot(dados$peso, dist="norm")   

# Testes de normalidade de Shapiro-Wilk

# Ho: a distribuicao e Normal
# H1: a distribuicao não e Normal


shapiro.test(dados$peso)      # Shapiro Willks - Mais rigoroso                         

ks.test(dados$peso,"pnorm", mean(dados$peso), sd(dados$peso))    # Kolmogorov-Smirnov, para n>50


# Ajustando o modelo 
#--------------------

mod.1 <- lm(peso~idade, data=dados)

summary(mod.1)


# Inserindo a reta ajustada no diagrama de dispersão
#----------------------------------------------------

plot(peso~idade, xlab="Idade (anos)", ylab="Peso (libras)", col="blue", data=dados)
abline(mod.1, col="red")


# Extraindo objetos do modelo
#-------------------------------

objects(mod.1)


# Coeficientes da reta de regressão ajustada
#-------------------------------------------

coefficients(mod.1)  
# ou
coef(mod.1)  
# ou   
mod.1$coefficients    

b0 <- coef(mod.1)[1]       # bo
b1 <- coef(mod.1)[2]       # b1

b0; b1

b0 + b1*8 # Peso previsto para uma criança com 8 anos 

b0 + b1*10 # Peso previsto para uma criança com 10 anos 


# Tabela de ANOVA
#--------------------

anova(mod.1)	         # Cuidado para não digitar Anova (a saída deste comando é um pouco diferente)


# Somas de quadrados
#--------------------

anova(mod.1)[1,4]

SQM <- anova(mod.1)[1,2]
SQR <- anova(mod.1)[2,2]
SQT <- SQM + SQR

SQM; SQR; SQT

SQM/SQT

# Graus de liberdade
#--------------------

gl.m <- anova(mod.1)[1,1] 
gl.r <- anova(mod.1)[2,1]
gl.t <- anova(mod.1)[1,1] + anova(mod.1)[2,1] 

gl.m; gl.r; gl.t 

# Quadrados médios
#--------------------

QMM <- SQM / gl.m
QMR <- SQR / gl.r

QMM; QMR

# Teste F
#--------------------

# Ho: Beta1 = 0
# H1: Beta1 <> 0

F1.10 <- QMM/QMR      # Calculando o valor da estatística F
F1.10

pvalor.F <- pf(F1.10, df1=gl.m, df2=gl.r, lower.tail=FALSE)   # Calculando o p-valor
pvalor.F

# Teste T
#--------------------

# Ho: Beta1 = 0
# H1: Beta1 <> 0

idade.m <- mean(dados$idade)
idade.m

idade.sq <- sum((dados$idade-idade.m)^2)  # Soma de quadrados (X - Xb)^2
idade.sq

b1

b1.var = QMR/idade.sq    # Calculando a variância de b1
b1.var

b1.dp <- sqrt(b1.var)    # Calculando o desvio padrão de b1
b1.dp

t10 <-  b1/b1.dp         # Calculando o valor da estatística t
t10

pvalor.t <- 2*pt(t10, 10, lower.tail = FALSE)     # Calculando o p-valor
pvalor.t

# Note que: 
# F1.10 = (t10)^2
# pvalor.F = pvalor.t

F1.10
t10^2

pvalor.F
pvalor.t

# No modelo simples, os testes F e T são equivalentes


# Coeficiente de determinação (R2)
#-------------------------------------------

R2 <- SQM / SQT    # coeficiente de determinação
R2
sqrt(R2)           # coeficiente de correlação


# Valores ajustados
#-------------------

30.571429 + 3.642857*dados$idade   # valores ajustados

b0 + b1*dados$idade                # valores ajustados

fitted(mod.1)                      # Valores ajustados


# Intervalo de confiança para os parâmetros do modelo
#-----------------------------------------------------

t.95 <- qt(0.025, 10, lower.tail = FALSE)
t.95

b1 - t.95*b1.dp
b1 + t.95*b1.dp

confint(mod.1, level=0.95)


# Intervalo de confiança para E(Y) - Intervalo de estimação
#------------------------------------------------------------


# Para idade = 8

xh <- 8

yh <- b0 + b1*xh               # valor ajustado
yh

yhe.var <- QMR*( 1/12 + ((xh-idade.m)^2/idade.sq) )
yhe.var

yhe.dp <- sqrt(yhe.var)
yhe.dp

# Intervalo de estimação

yh - t.95*yhe.dp
yh + t.95*yhe.dp


# Usando a função predict


nd <- data.frame( idade = seq(6,12,0.5) )
nd

pr.co <- predict( mod.1, newdata=nd, interval="conf" )
pr.co


# Intervalo de confiança para Y - Intervalo de predição
#-------------------------------------------------------
# Para idade = 8

xh <- 8

yh <- b0 + b1*xh               # valor ajustado
yh

yhp.var <- QMR*( 1 + 1/12 + ((8-idade.m)^2/idade.sq) )
yhp.var

yhp.dp <- sqrt(yhp.var)
yhp.dp

# Intervalo de predição

yh - t.95*yhp.dp
yh + t.95*yhp.dp


# Usando a função predict

nd <- data.frame( idade = seq(6,12,0.5) )
nd

pr.pr <- predict( mod.1, newdata=nd, interval="pred" )
pr.pr

plot(peso~idade, xlab="Idade (anos)", ylab="Peso (libras)", col="blue", data=dados)
matlines( nd$idade, pr.co, lty=1, lwd=c(5,2,2), col="blue" )
matlines( nd$idade, pr.pr, lty=2, lwd=c(5,2,2), col="blue" )



#----------------------------------------------------------------------------------
#                               ANÁLISE DE RESÍDUOS
#----------------------------------------------------------------------------------

dados$idade
dados$peso

dados$peso - (b0 + b1*dados$idade)   # Resíduos brutos

residuals(mod.1)                     # Resíduos brutos

# Obtendo o resíduo padronizado
#-------------------------------

mean(residuals(mod.1))         # media dos resíduos brutos

(residuals(mod.1) - mean(residuals(mod.1))) / sqrt(QMR)    # Resíduos semi-estudentizados

rstandard(mod.1)      # Resíduos estudentizados


# Gráficos importantes
#----------------------


# Diagrama de dispersão dos resíduos x variável explicativa

plot(dados$idade, rstandard(mod.1) , xlab="Idade (anos)", ylab="Resíduo Padronizado")
abline(h=0, col="blue")
abline(h=-1.96, col="blue", lty=3)
abline(h= 1.96, col="blue", lty=3)

help(par)

# ylim=c(-3, 3))


# Diagrama de dispersão dos resíduos x valores ajustados

plot(fitted(mod.1), rstandard(mod.1), xlab="Valores ajustados", ylab="Resíduo Padronizado")
abline(h=0, col="blue")
abline(h=-1.96, col="blue", lty=3)
abline(h= 1.96, col="blue", lty=3)


# Box-plot, histograma e qq-plot dos resíduos

boxplot(rstandard(mod.1), ylab="Resíduo Padronizado")

hist(rstandard(mod.1), scale="density",  xlab="Resíduo Padronizado", ylab="Densidade de frequência", main="")

qqPlot(rstandard(mod.1), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Padronizado")


# Default do R

par(mfrow=c(2,2))
plot(mod.1)

par(mfrow=c(1,1))


# DESCREVENDO A RELACAO/ASSOCIACAO ENTRE PESO E ALTURA
#----------------------------------------------------------------------------------

cor(dados$peso, dados$altura)
cor.test(dados$peso, dados$altura, alternative="two.sided", method="pearson")

plot(peso~altura, xlab="Altura (pés)", ylab="Peso (libras)", col="blue", data=dados)

mod.2 <- lm(peso~altura, data=dados)
summary(mod.2)
anova(mod.2)
confint(mod.2)

plot(peso~altura, xlab="Altura (pés)", ylab="Peso (libras)", col="blue", data=dados)
abline(mod.2, col="red")

# Intervalo de confiança para E(Y)
#---------------------------------

summary(dados$altura)

nd <- data.frame( altura = 42:62 )
nd

pr.co <- predict( mod.2, newdata=nd, interval="conf" )
pr.co


# Intervalo de predicao para Y
#---------------------------------

pr.pr <- predict( mod.2, newdata=nd, interval="pred" )
pr.pr

plot(peso~altura, xlab="Altura (pes)", ylab="Peso (libras)", col="blue", data=dados)
matlines( nd$altura, pr.co, lty=1, lwd=c(5,2,2), col="blue" )
matlines( nd$altura, pr.pr, lty=2, lwd=c(5,2,2), col="blue" )


# Análise de resíduos
#---------------------------------

# Diagrama de dispersão dos resíduos x variável explicativa

plot(dados$altura, rstandard(mod.2) , xlab="Altura (pés)", ylab="Resíduo Padronizado")
abline(h=0, col="blue")
abline(h=-1.96, col="blue", lty=3)
abline(h= 1.96, col="blue", lty=3)

# Diagrama de dispersão dos resíduos x valores ajustados

plot(fitted(mod.1), rstandard(mod.2), xlab="Valores ajustados", ylab="Resíduo Padronizado")
abline(h=0, col="blue")
abline(h=-1.96, col="blue", lty=3)
abline(h= 1.96, col="blue", lty=3)


# Box-plot, histograma e qq-plot dos resíduos

boxplot(rstandard(mod.2), ylab="Resíduo Padronizado")

hist(rstandard(mod.2), scale="density",  xlab="Resíduo Padronizado", ylab="Densidade de frequência", main="")

qqPlot(rstandard(mod.2), dist="norm", xlab="Percentis da N(0,1)", ylab="Resíduo Padronizado")


# Default do R

par(mfrow=c(2,2))
plot(mod.2)

par(mfrow=c(1,1))