Untitled

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,]

#tranformar 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])
  }
}

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_inicial))
dev.off()
corrplot(a, method = 'color', addCoef.col = 'black', 
         number.cex = 0.45, tl.cex = 0.5, tl.col="black")

#eliminar Systolic e Diastolic Pressure
dados_inicial <- dados_inicial[, -which(colnames(dados_inicial) %in% c("Diastolic.Pressure", "Systolic.Pressure"))]

#modelo completo com covariaveis eliminadas
mrl.comp <-lm(y ~ .,data=dados_inicial)
plot(mrl.comp)
summary(mrl.comp)
AIC(mrl.comp)

library(gvlma)
gvlma(mrl.comp)

#boxcox para por y normal
library(MASS)
b <- boxcox(lm(y ~ 1))
lambda <- b$x[which.max(b$y)]
y2 <- y^lambda #lambda=0.30303...


#modelo com y normal
mrl.comp2 <-lm(y2 ~ .,data=dados_inicial)
plot(mrl.comp2) #mais proximo de normal


#eliminar outliers - nao tem pelo box plot 
boxplot(y2)

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

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

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

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

#novo modelo completo
mrl.full <- lm(y_treino~.,data=dados_treino)
summary(mrl.full)
AIC(mrl.full)

#stepforward
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)
    #Shock.Type + Sex + Age

#stepbackward
mrl.stepbackward <- step(mrl.full, direction = "backward", trace=FALSE)

summary(mrl.stepbackward)
    #Age + Sex + Shock.Type 


predictions <- predict(mrl.full, data = dados_treino)
#plot(predictions, y_treino, type = "p")
d<-data.frame(predictions, y=y_treino)
ggplot(d, aes(predictions, y, color=predictions)) +
  geom_point(shape = 16, size = 5, show.legend = FALSE) +
  theme_minimal()

###
betahat = solve(t(X) %*% X) %*% t(X) %*% y
betahat

SSE = t(y - X %*% betahat) %*% (y - X %*% betahat)
MSE = SSE/(n-p)

varbetahat = c(MSE) * solve(t(X) %*% X)

fittedval = X %*% betahat
simpleres = y - fittedval
simpleres2 = y^lambda - fittedval^lambda

w = cbind(fittedval, simpleres)
colnames(w) = c("fitted values","residuals")
print(w)



C <- cbind(0,diag(p-1))
m <- 0
SSE_h <- t(y)%*%(diag(n) - X%*%solve(t(X)%*%X)%*%t(X))%*%y
f <- c()
for(i in 2:p-1){ 
  ci <- t(matrix(C[i,])) #linha
  Q_i <- t(ci%*%betahat-m)%*%solve(ci%*%solve(t(X)%*%X)%*%t(ci))%*%(ci%*%betahat-m)
  f_i <- (Q_i/1)/(SSE_h/(n-p))
  f <- c(f,f_i)
}
1-pf(f,1,n-p)   

library(ggplot2)