Untitled

dados <- read.delim("Trabalho10_Dados.txt", header= FALSE, sep="", )
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,]

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

X <- cbind(1,data.matrix(dados_inicial[,-ncol(dados_inicial)]))

y <- matrix(dados_inicial$Cardiac.Index)
n=length(y)

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

summary(dados_inicial)

#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 Sistolic e Diastolic Pressure
X <- X[, -which(colnames(X) %in% c("Diastolic.Pressure", "Systolic.Pressure"))]
p=ncol(X)

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


#R^2 nao devia dar tão baixo!!!???
mrl <- lm(formula = y ~ Age + Height + Sex + Survival + Shock.Type + 
            Mean.Arterial.Pressure + Heart.Rate + Mean.Central.Venous.Pressure +
            Body.Surface.Area, dados_inicial)
summary(mrl)

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(MASS)
b <- boxcox(lm(y ~ 1))
lambda <- b$x[which.max(b$y)]
hist(y^lambda)