library(compareGroups)
library(formattable)
library(tidyverse)
library(epitools)
library(gmodels)
#library(psychometric)

##### Cargar los datos
##### Primero definid el directorio de trabaj0: Session > Set Working Directory
remove(uci)
uci <- readRDS(file='uci.R')
uci[1:20,] %>% formattable()

##### Assignar etiquetas a los códigos de las variable categóricas

uci$SEXO <- factor(uci$SEXO,labels=c('Home','Dona'))
uci$ESTADO <- factor(uci$ESTADO,labels=c('Viu','Mort'))
uci$VME <- factor(uci$VME,labels=c('No','Si'))
uci %>% head() %>% formattable()

##### Descriptiva en función del ESTADO final del paciente

t <- compareGroups(VME~ESTADO,data=uci)
home <- createTable(update(t,subset=SEXO=='Home'),show.ratio = TRUE)
dona <- createTable(update(t,subset=SEXO=='Dona'),show.ratio = TRUE)
home
dona

rbind(home,dona)

compareGroups(SEXO~ESTADO,data=uci) %>% createTable(show.ratio = T)


with(uci,CrossTable(VME,ESTADO,
                    asresid=T,format='SPSS',
                    prop.r=F,
                    prop.t=F,
                    prop.chisq=F,
                    chisq = T))

###### Cáclulo de OR (la referencia es el grupo de vivos)

compareGroups(ESTADO~.-ID,data=uci,ref=1) %>% createTable(show.ratio = T)


##### Gráfica de OR univariantes

DataOR <- matrix(c(1,1.04,1.03,1.06,
                   2,0.72,0.48,1.06,
                   3,0.99,0.98,0.99,
                   4,0.98,0.97,0.99,
                   5,5.54,3.68,8.51,
                   6,1.02,1,1.03,
                   7,1.01,0.98,1.04,
                   8,0.09,0.02,0.41,
                   9,0.97,0.94,0.99),ncol = 4,byrow = T)
DataOR
 
colnames(DataOR) <- c('Group','OR','low','upper')
DataOR <- as.data.frame(DataOR)
DataOR$Group <- factor(DataOR$Group,labels=c('Edat','Sexo(Dona)','TAS','TAD','VME(si)','Estada','Sodio','pH','CO3H'))
DataOR
ggplot(DataOR,aes(x=Group,y=OR))+
  geom_errorbar(aes(ymin=low,ymax=upper),color='blue',size=1,width=.4)+
  geom_point(shape=21,fill='white',size=4)+
  coord_flip()+
  geom_hline(yintercept = 1,color='red',linetype='dashed',size=1)+
  ylim(0,10)

##### Relación entre Hipertension y Estado

df <- uci %>% mutate(TAM = round(TAD+(TAS-TAD)/3,1)) %>% 
  mutate(TAMR = cut(TAM,c(0,50,1000)))
df$TAMR <- factor(df$TAMR,labels=c('Hipotenso','Normal'))
df[1:20,] %>% formattable()

compareGroups(TAMR~ESTADO,data=df,ref=2) %>% createTable(show.ratio = T)

#### Gráfica de barras

t <- with(df,table(ESTADO,TAMR))
tp <- prop.table(t,2) %>% round(2)
tp 
tp <- as.data.frame(tp)
tp

ggplot(tp,aes(x=ESTADO,y=Freq,fill=TAMR))+
  geom_bar(stat='identity',position=position_dodge())+
  geom_text(aes(label=Freq),vjust=1.5,color="white",
            position = position_dodge(0.9), size=5)+
  ylim(0,1)

(643/357)/(311/689)

###### Oddsratio

compareGroups(TAMR~ESTADO,data=df,ref=2) %>% createTable(show.ratio = T)

with(df,CrossTable(ESTADO,TAMR,
                    asresid=T,format='SPSS',
                    prop.r=F,
                    prop.t=F,
                    prop.chisq=F,
                    chisq = T))
###### pH y Co3H

head(df)

ggplot(df,aes(x=PH,y=CO3H))+
  geom_smooth(method='lm',se=F)+
  geom_point()

rv <-with(df,cor(PH,CO3H))
CIr(r=rv, n = length(df$PH), level = .95)


res <- lm(CO3H~PH,data=df)     
summary(res)
confint(res)
anova(res)

cbind(coefficients(res),confint(res)) %>% formattable()

df %>% group_by(SEXO,ESTADO) %>% summarise(r=round(cor(PH,CO3H),2),
                                           low=round(CIr(r=cor(PH,CO3H), n = length(df$PH), level = .95),2)[1],
                                           upper=round(CIr(r=cor(PH,CO3H), n = length(df$PH), level = .95),2)[2]
                                          ) %>%
      formattable()


ggplot(df,aes(x=PH,y=CO3H))+
  geom_smooth(method='lm',se=F)+
  geom_point()+
  facet_wrap(SEXO~ESTADO)

res2 <- lm(CO3H~PH+SEXO*ESTADO,data=df)     
summary(res2)
anova(res2)

confint(res2)



res3 <- lm(CO3H~PH+SEXO,data=df)     
summary(res3)
anova(res3)
confint(res3)

anova(res,res3,res2)

ggplot(df,aes(x=PH,y=CO3H,fill=SEXO))+
  geom_smooth(method='lm',se=F,aes(color=SEXO))+
  geom_point(aes(color=SEXO))

coef <- coefficients(res3)

######## Gráfica del modelo final CO3H vs PH por SEXO

ggplot(df,aes(x=PH,y=CO3H))+
  geom_point(aes(color=SEXO,shape=SEXO),size=2)+
  geom_abline(intercept = -98.65,slope=16.6,color='red')+
  geom_abline(intercept = -98.65-1.2,slope=16.6,color='green')

####### Caracterización de outliers

ggplot(df,aes(x=SEXO,y=CO3H,fill=SEXO))+
  geom_boxplot()

rr <- lm(CO3H~PH*ESTADO,data=df)

summary(rr)
coefficients(rr)

######### Relación entre la estancia y la mortalidad

summary(df$ESTADA)

df <- df %>% mutate(GESTADA=cut(ESTADA,c(0,3,7,14,200)))
table(df$GESTADA) 
df$GESTADA <- factor(df$GESTADA,labels=c('<3 dias','4-7 dias','8-14 dias','+14 dias'))

dt<-df %>% group_by(GESTADA,ESTADO) %>% summarise(n=n()) %>% mutate(freq=n/sum(n))

with(df,CrossTable(GESTADA,ESTADO,
                   asresid=T,format='SPSS',
                   prop.c=F,
                   prop.t=F,
                   prop.chisq=F,
                   chisq = T))

createTable(compareGroups(GESTADA~ESTADO,df),show.p.trend = T)
dt
ggplot(dt,aes(x=GESTADA,y=freq,fill=ESTADO))+
  geom_bar(stat='identity',
           position=position_dodge(),
           color='black')+scale_fill_manual(values=c('brown','orange'))+
  ylim(0,1)

createTable(compareGroups(ESTADO~GESTADA,data=df),show.ratio = T)
createTable(compareGroups(ESTADO~GESTADA,data=df,ref=3),show.ratio = T)

t1 <- createTable(compareGroups(ESTADO~GESTADA,data=df,subset = SEXO== 'Home'),show.ratio = T)
t2 <- createTable(compareGroups(ESTADO~GESTADA,data=df,subset = SEXO== 'Dona'),show.ratio = T)
cbind('Homes'=t1,'Dones'=t2)

###################### Sodio 

t1 <- createTable(compareGroups(ESTADO~SODIO,data=df,subset = SEXO== 'Home'))
t2 <- createTable(compareGroups(ESTADO~SODIO,data=df,subset = SEXO== 'Dona'))
cbind('Homes'=t1,'Dones'=t2)

ggplot(df,aes(x=ESTADO,y=SODIO,fill=SEXO))+geom_boxplot()

res <- lm(SODIO~ESTADO*SEXO,data=df)
anova(res)

############## Grupos de nivelde socio

df <- df %>% mutate(GSODIO = cut(SODIO,c(0,135,145,1000)))
table(df$GSODIO)
df$GSODIO <- factor(df$GSODIO,labels=c('Hiponatremia','Normal','Hipernatremia'))

createTable(compareGroups(GSODIO~ESTADO,data=df,ref=2),show.ratio = T)
createTable(compareGroups(ESTADO~GSODIO,data=df,ref=2),show.ratio = T)

df$GSODIO <- relevel(df$GSODIO,ref='Normal')
res <- glm(ESTADO~GSODIO+SEXO,family='binomial',data=df)
summary(res)
exp(confint(res))

#################### REGRESION LOGISTICA

res <- glm(ESTADO~.-ID,family='binomial',data =df)
summary(res)$coefficients %>% formattable()

res <- glm(ESTADO~EDAT+VME+TAMR+GESTADA+GSODIO,family='binomial',data =df)
summary(res)

df$TAMR <- relevel(df$TAMR,ref='Normal')
res <- glm(ESTADO~EDAT+VME+TAMR+GSODIO,family='binomial',data =df)
summary(res)

OR<-summary(res)$coefficients[,1] %>% exp()
r2<-exp(confint(res)) %>% round(2)
res.c <- cbind(OR,r2) %>% round(2)
res.c <- as.data.frame(res.c)
res.c

newdata <- data.frame(EDAT=50,VME='Si',TAMR='Hipotenso',GSODIO='Hiponatremia')
ff1 <- function(x) predict(res, data.frame(EDAT=x,VME='Si',TAMR='Hipotenso',GSODIO='Hiponatremia'), type="response")
ff2 <- function(x) predict(res, data.frame(EDAT=x,VME='No',TAMR='Hipotenso',GSODIO='Hiponatremia'), type="response")

punts <- seq(30,90,5)
d1 <- data.frame(edat=punts,p=ff1(punts))
d2 <- data.frame(edat=punts,p=ff2(punts))

ggplot(d1,aes(x=edat,y=p))+
  geom_line(size=1,color='red')+
  geom_point(shape=21,fill='white',size=4)+
  geom_line(data=d2,aes(x=edat,y=p),size=1,color='blue')+
  geom_point(data=d2,shape=21,fill='white',size=4)
  


ff1 <- function(x) predict(res, data.frame(EDAT=x,VME='No',TAMR='Normal',GSODIO='Hiponatremia'), type="response")
ff2 <- function(x) predict(res, data.frame(EDAT=x,VME='Si',TAMR='Normal',GSODIO='Hiponatremia'), type="response")
ff3 <- function(x) predict(res, data.frame(EDAT=x,VME='No',TAMR='Normal',GSODIO='Normal'), type="response")
ff4 <- function(x) predict(res, data.frame(EDAT=x,VME='Si',TAMR='Normal',GSODIO='Normal'), type="response")

punts <- seq(30,90,5)
d1 <- data.frame(edat=punts,p=ff1(punts))
d2 <- data.frame(edat=punts,p=ff2(punts))
d3 <- data.frame(edat=punts,p=ff3(punts))
d4 <- data.frame(edat=punts,p=ff4(punts))

ggplot(d1,aes(x=edat,y=p))+
  geom_line(size=1,color='red')+
  geom_point(shape=21,fill='white',size=4)+
  geom_line(data=d2,aes(x=edat,y=p),size=1,color='red',linetype='dashed')+
  geom_point(data=d2,shape=21,fill='white',size=4)+
  geom_line(data=d3,aes(x=edat,y=p),size=1,color='blue')+
  geom_point(data=d3,shape=21,fill='white',size=4)+
  geom_line(data=d4,aes(x=edat,y=p),size=1,color='blue',linetype='dashed')+
  geom_point(data=d4,shape=21,fill='white',size=4)