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

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

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

#Teste
dados_teste <- dados_inicial[teste_ind,]

#matrizes x e y
X <- cbind(1,data.matrix(dados_treino[,-ncol(dados_treino)]))
y <- matrix(dados_treino$Cardiac.Index)
n=length(y)

#box.plot continuas
dados_cont <- dados_treino[ ,-c(match(qualitativas, names(dados_treino)))]
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])) 
summary(dados_treino)

#correlação
library(corrplot)
a = cor(X[,2:ncol(X)])
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
X <- X[, -which(colnames(X) %in% c("Diastolic.Pressure", "Systolic.Pressure"))]
p=ncol(X)

#boxcox para por y normal
library(MASS)
b <- boxcox(lm(y ~ 1))
lambda <- b$x[which.max(b$y)]
shapiro.test(y^lambda) #p-value>0.05 logo y está normal 

#modelo completo com covariaveis eliminadas
mrl.comp <- lm(formula = y ~ Age + Height + Sex + Survival + Shock.Type + 
                 Mean.Arterial.Pressure + Heart.Rate + Mean.Central.Venous.Pressure +
                 Body.Surface.Area, dados_treino)
summary(mrl.comp)
plot(mrl.comp)
plot(predict(mrl.comp), y)
AIC(mrl.comp)





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

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)


mrl2 <- lm(formula = y^lambda ~ Age + Height + Sex + Survival + Shock.Type + 
             Mean.Arterial.Pressure + Heart.Rate + Mean.Central.Venous.Pressure +
             Body.Surface.Area, dados_inicial)
summary(mrl2)

residuals(mrl2)
plot(mrl2)

influence(mrl2)

AIC(mrl2)