Untitled

#TRATAMENTO DE DADOS
dados <- read.delim("Trabalho10_Dados.txt", header= FALSE, sep="")

#escolha dos dados pedidos no enunciado
dados <- dados[,2:13]
names(dados)<- c("Age", "Height", "Sex", "Survival", "Shock.Type", "Systolic.Pressure", 
                 "Mean.Arterial.Pressure", "Heart.Rate", "Diastolic.Pressure", "Mean.Central.Venous.Pressure", 
                 "Body.Surface.Area", "Cardiac.Index")
rows <- nrow(dados)
even_rows <- seq_len(rows) %% 2
dados_inicial <- dados[even_rows == 1,]
rownames(dados_inicial) <- NULL #renumerar linhas
nrows <- nrow(dados_inicial)

#transformar qualitativas
qualitativas <- c("Sex", "Survival", "Shock.Type")
for (i in 1:ncol(dados_inicial)){
  if (names(dados_inicial[i]) %in% qualitativas){
    dados_inicial[,i] <- factor(dados_inicial[,i])
  }
}
#ANALISE DE DADOS
summary(dados_inicial)

#box.plot continuas
dados_cont <- dados_inicial[ ,-c(match(qualitativas, names(dados_inicial)))]
par(mfrow = c(3, ncol(dados_cont)/3), mar=c(2,2,2,2), cex=0.5)
lapply(1:ncol(dados_cont), function(i) boxplot(dados_cont[,i], main=names(dados_cont)[i])) 

#variavel resposta
y <- matrix(dados_inicial$Cardiac.Index)
dados_inicial <- dados_inicial[,-ncol(dados_inicial)]

#correlações
library(corrplot)
a <- cor(data.matrix(dados_cont))
dev.off()
corrplot(a, method = 'color', addCoef.col = 'black', 
         number.cex = 0.45, tl.cex = 0.5, tl.col="black")

#eliminar Mean Arterial
mrl.comp <-lm(y ~ .,data=dados_inicial)
library(car)
library(usdm)
vif(mrl.comp)
dados_inicial <- dados_inicial[, -which(colnames(dados_inicial) %in% c("Mean.Arterial.Pressure"))]
mrl.comp <-lm(y ~ .,data=dados_inicial)
vif(mrl.comp)


#ANALISE PRELIMINAR DO MODELO COMPLETO

##dados de treino e de teste
set.seed(73)             
teste_ind<-sort(sample(nrows,0.2*nrows))

#Treino
dados_treino <-dados_inicial[-teste_ind,]
y_treino <- y[-teste_ind]
summary(dados_treino)

#Teste
dados_teste <- dados_inicial[teste_ind,]
y_teste <- y[teste_ind]

#modelo completo com covariaveis eliminadas e dados de treino
mrl.comp2 <-lm(y ~ .,data=dados_treino)
#plot(mrl.comp2)
summary(mrl.comp2)
extractAIC(mrl.comp2)

library(gvlma)
gvlma(mrl.comp2)

#boxcox para por ys normais ???
library(MASS)
shapiro.test(y_treino)
shapiro.test(y_teste)

b_treino <- boxcox(lm(y_treino ~ 1))
lambda_treino <- b_treino$x[which.max(b_treino$y_treino)] #lambda=0
y_treino <- log(y_treino)

b_teste <- boxcox(lm(y_teste ~ 1))
lambda_teste <- b_teste$x[which.max(b_teste$y_teste)] #lambda=0
y_teste <- log(y_teste)

#teste de ajustamento dos ys
shapiro.test(y_treino) ### p-value = 0.003804 !!!!!!!!!!
shapiro.test(y_teste) ### p-value = 0.9324

#modelo com y normal
mrl.comp3 <-lm(y_treino ~ .,data=dados_treino)
plot(mrl.comp3) #mais proximo de normal

#eliminar outliers - nao tem significativos pelo box plot - eliminar mais tarde
boxplot(y_treino)
boxplot(y_teste)

##capacidade de previsão
library(ggplot2)
mrl.comp3_teste <- lm(y_teste ~ .,data=dados_teste)
predictions_comp3 <- predict(mrl.comp3_teste, newdata = dados_teste)
d2<-data.frame(predictions_comp3, y=y_teste)
ggplot(d2, aes(predictions_comp3, y)) +
  geom_point(shape = 16, size = 3, show.legend = FALSE)  +
  geom_abline(intercept = 0, slope = 1, color = "red")
ri <- y_teste - predictions2
sum(ri^2)/length(y_teste)


#matrizes x 
X <- cbind(1,data.matrix(dados_treino))
n=length(y_treino)

mrl.full <- lm(y_treino~.,data=dados_treino)
summary(mrl.full)
extractAIC(mrl.full)
gvlma(mrl.full)
#anova??

#stepforward
#primeira int
mrl.base <-lm(y_treino~1,
              data=dados_treino)

mrl.stepforward <- step(mrl.base,
                        scope = list(upper = formula(mrl.full), 
                                     lower = formula(mrl.base)),
                        direction = "forward", trace="FALSE")

summary(mrl.stepforward)
extractAIC(mrl.stepforward)
formula(mrl.stepforward)
anova(mrl.stepforward,mrl.full)

#segunda int
mrl.int2 <- lm(y_treino ~ (Shock.Type + Sex + Diastolic.Pressure + Systolic.Pressure + 
                             Age + Body.Surface.Area)^2, 
               data = dados_treino)

mrl.stepforward2 <- step(mrl.stepforward,
                         scope = list(upper = formula(mrl.int2), 
                                      lower = formula(mrl.stepforward)),
                         direction = "forward", trace="FALSE")
summary(mrl.stepforward2)
extractAIC(mrl.stepforward2)
formula(mrl.stepforward2)
anova(mrl.stepforward2,mrl.int2)

#terceira int
mrl.int3 <- update(mrl.stepforward2, .~. + Sex:Body.Surface.Area:Diastolic.Pressure)

mrl.stepforward3 <- step(mrl.stepforward2,
                         scope = list(upper = formula(mrl.int3), 
                                      lower = formula(mrl.stepforward2)),
                         direction = "forward", trace="FALSE")
summary(mrl.stepforward3)
extractAIC(mrl.stepforward3)
formula(mrl.stepforward3)
anova(mrl.stepforward3,mrl.int3)

#bestsubset com r^2 ajustado
library(leaps)
bestsubsets <- regsubsets(y_treino~., data = dados_treino, nvmax = ncol(dados_treino))  
#plot(bestsubsets, scale = "adjr2")
bs_summary <- summary(bestsubsets) 
ind_r2 <- which.max(bs_summary$adjr2)
coef.r <- coef(bestsubsets, ind_r2) 
#Age, Sex, Shock.Type, Systolic.Pressure, Diastolic.Pressure, Body.Surface.Area

mrl.best.r<- lm(y_treino ~ Age + Sex + Shock.Type + Systolic.Pressure +
                  Diastolic.Pressure + Body.Surface.Area, data = dados_treino)

summary(mrl.best.r)
extractAIC(mrl.best.r)
formula(mrl.best.r)
anova(mrl.best.r, mrl.full)

#com interações 2 a 2

bestsubsets2 <- regsubsets(y_treino ~ (Age + Sex + Shock.Type + Systolic.Pressure +
                                         Diastolic.Pressure + Body.Surface.Area)^2, 
                           data = dados_treino, nvmax = 21, force.in = c(1,3,5,6,8,10))

plot(bestsubsets2, scale = "adjr2")
bs_summary2 <- summary(bestsubsets2) 
ind_r22 <- which.max(bs_summary2$adjr2)
coef.r2 <- coef(bestsubsets2, ind_r22) 


#Age, Sex, Shock.Type, Systolic.Pressure, Diastolic.Pressure, Body.Surface.Area,
#Age:Shock.Type, Sex:Systolic.Pressure, Sex:Body.Surface.Area, Shock.Type:Systolic.Pressure,
#Shock.Type:Diastolic.Pressure, Diastolic.Pressure:Body.Surface.Area


mrl.best.r2<- update(mrl.best.r, .~. + Age:Shock.Type + Sex:Systolic.Pressure + 
                       Sex:Body.Surface.Area + Shock.Type:Systolic.Pressure  + 
                       Shock.Type:Diastolic.Pressure + Diastolic.Pressure:Body.Surface.Area)

summary(mrl.best.r2)
extractAIC(mrl.best.r2)
formula(mrl.best.r2)
anova(mrl.best.r2,lm(y_treino ~ (Age + Sex + Shock.Type + Systolic.Pressure +  Diastolic.Pressure + Body.Surface.Area)^2,  data = dados_treino))
#não há interações 3 a 3

bestsubsets <- regsubsets(y_treino~., data = dados_treino, nvmax = ncol(dados_treino))  
#plot(bestsubsets, scale = "adjr2")
bs_summary <- summary(bestsubsets) 
ind_r2 <- which.max(bs_summary$adjr2)
coef.r <- coef(bestsubsets, ind_r2) 
#Age, Sex, Shock.Type, Systolic.Pressure, Diastolic.Pressure, Body.Surface.Area

mrl.best.r<- lm(y_treino ~ Age + Sex + Shock.Type + Systolic.Pressure +
                  Diastolic.Pressure + Body.Surface.Area, data = dados_treino)

###############
bestsubsets_cp <- regsubsets(y_treino~., data = dados_treino, nvmax = ncol(dados_treino))  
bs_summary <- summary(bestsubsets_cp) 
ind_cp <- which.max(bs_summary$cp)
coef.cp <- coef(bestsubsets_cp, ind_cp) 

mrl.best.cp <- lm(y_treino ~ Age + Sex + Shock.Type + Systolic.Pressure +
                  Diastolic.Pressure + Body.Surface.Area, data = dados_treino)