karina

# Método
###Estudo Transversal
require(RcmdrMisc)
require(dplyr)
library(colorspace, pos=33)

#notifcações match
banco <- read.csv("C:/Users/edson/Downloads/Base_de_dados_match  (1).csv", comment.char="#")
banco$recorrencia2<-0
banco$recorrencia2[banco$recorrencia!="1" & banco$diff<365]<-1
banco$recorrencia2<-as.factor(banco$recorrencia2)
banco$CS_GESTANT<-as.factor(banco$CS_GESTANT)
banco$CS_ESCOL_N<-as.factor(banco$CS_ESCOL_N)
#repetição

## Load survival package
library(survival)

#remoção de duplicatas

dataset    <- banco
dataset$NU_IDADE_N<-dataset$NU_IDADE_N-4000
#dataset$recorrencia2<-0
#dataset$recorrencia2[dataset$recorrencia>2 & dataset$diff<60]<-1
#banco$recorrencia2<-as.factor(banco$recorrencia2)

dataset<-arrange(dataset,desc(dataset$DT_NOTIFIC))

dataset<-dplyr::distinct(dataset,dataset$study_ID,.keep_all=T)




tabela <- within(dataset, {
CS_ESCOL_N <- as.factor(CS_ESCOL_N)
CS_GESTANT <- as.factor(CS_GESTANT)
CS_RACA <- as.factor(CS_RACA)
MES <- as.factor(MES)
PCRUZ <- as.factor(PCRUZ)
RESULT <- as.factor(RESULT)
recorrencia2 <- as.factor(recorrencia2)

})

dataset$TRA_ESQUEM<-as.factor(dataset$TRA_ESQUEM)

#criar variável recorrência

tabela<-subset(dataset, select=c(NU_IDADE_N , CS_SEXO,     CS_GESTANT,  CS_RACA,     CS_ESCOL_N,  RESULT,     
PMM       ,  PCRUZ   ,    TRA_ESQUEM , MES     ,      recorrencia2,
diff,        LVC ,        Tempo_class))
#teste de hiotese em banco bruto as.it is
#require(gtsummary)
#gtsummary::tbl_summary(tabela, by= "CS_SEXO")
#gtsummary::tbl_summary(tabela,by="recorrencia2")

#levando em conta subset de recorrências
#Script de Cox Exaustivo

options(digits=3)





dataset <- within(dataset, {
CS_ESCOL_N <- as.factor(CS_ESCOL_N)
CS_GESTANT <- as.factor(CS_GESTANT)
CS_RACA <- as.factor(CS_RACA)
MES <- as.factor(MES)
PCRUZ <- as.factor(PCRUZ)
RESULT <- as.factor(RESULT)
TRA_ESQUEM <- as.factor(TRA_ESQUEM)
})

dataset$recorrencia2[dataset$recorrencia2==2]<-0

## Create outcome variable, then survival object
dataset <- within(dataset, {
  status.dichotomous <- recorrencia2 != 0
  survival.vector    <- Surv(diff, recorrencia2)
})

dataset$TRA_ESQUEM2<-0
dataset$TRA_ESQUEM2[dataset$TRA_ESQUEM=="4"]<-1
dataset$TRA_ESQUEM2<-as.factor(dataset$TRA_ESQUEM2)

dataset$TRA_ESQUEM3<-0
dataset$TRA_ESQUEM3[dataset$TRA_ESQUEM=="8"]<-1
dataset$TRA_ESQUEM3<-as.factor(dataset$TRA_ESQUEM3)

dataset$TRA_ESQUEM4<-0
dataset$TRA_ESQUEM4[dataset$TRA_ESQUEM=="10"]<-1
dataset$TRA_ESQUEM4<-as.factor(dataset$TRA_ESQUEM4)


dataset$TRA_ESQUEM5<-0
dataset$TRA_ESQUEM5[dataset$TRA_ESQUEM=="3"]<-1
dataset$TRA_ESQUEM5<-as.factor(dataset$TRA_ESQUEM5)


dataset$gestante2<-0
dataset$gestante2[dataset$CS_GESTANT!=1]<-1
dataset$gestante2<-as.factor(dataset$gestante2)

dataset$raca2<-0
dataset$raca2[dataset$CS_RACA!=1]<-1
dataset$raca2<-as.factor(dataset$raca2)

dataset$escola2<-0
dataset$escola2[dataset$CS_ESCOL_N=="7"|dataset$CS_ESCOL_N=="8"|dataset$CS_ESCOL_N=="9"|dataset$CS_ESCOL_N=="10"]<-1
dataset$escola2<-as.factor(dataset$escola2)


## Create vectors for outcome and predictors
outcome    <- c("survival.vector")
#predictors <- c("esquema2","cd8_cd4", "idade","cv","COR","gon","Grupo6")
predictors <- c('NU_IDADE_N' , 'CS_SEXO',  'raca2',     'escola2',  'RESULT',     'PMM'       ,  'PCRUZ'   ,    'TRA_ESQUEM2'   ,'TRA_ESQUEM3','TRA_ESQUEM4','TRA_ESQUEM5'  ,  'LVC')





## The lines below should not need modification.

## Create list of models
list.of.models <- lapply(seq_along((predictors)), function(n) {
  
  left.hand.side  <- outcome
  right.hand.side <- apply(X = combn(predictors, n), MARGIN = 2, paste, collapse = " + ")
  
  paste(left.hand.side, right.hand.side, sep = "  ~  ")
})

## Convert to a 

vector.of.models <- unlist(list.of.models)

## Fit coxph to all models
#list.of.fits <- lapply(vector.of.models, #function(x) {
#  
#  formula    <- as.formula(x)
#  fit        <- coxph(formula, data = #dataset, id=study_ID)
#  suma<- summary(fit)
#  result.AIC <- extractAIC(fit)
#  
#  data.frame(num.predictors = #result.AIC[1],
#             AIC            = #result.AIC[2],
#             model          = x,
#             suma= #suma[["coefficients"]])
#})


## Collapse to a data frame
#result <- do.call(rbind, list.of.fits)
#result$ciup<-exp(result$suma.coef + #2*result$suma.se.coef.)
#result$cidown<-exp(result$suma.coef - #2*result$suma.se.coef.)
#result$var<- row.names(result)


## Sort and print
#library(doBy)
#a<-orderBy(~ AIC + model, result)
#a<-as.data.frame(a, digits=2)
#result$pvalorsig<-0
#result$pvalorsig[result$suma.Pr...z..<0.05]<-1

#f<-result %>%
#group_by(model) %>%
#summarise(somamodelo= sum(pvalorsig))
#a<-merge(a,f)

#a$dif<-a$num.predictors- a$somamodelo




#filtrar modelos com todas as var em todos os fatores < 0.05
#filtrar modelos com mais de um preditor, vai ser todos uni se binária ou numérica todas as var em todos os fatores < 0.05
#modelomulticomtudosig<-subset(a,a$dif<1 & a$num.predictors>1)
##filtrar modelos com mais de um preditor, vai ser todos uni se binária ou numérica todas -1 as var em todos os fatores < 0.05, para pegar multinomial
#modelomulticomtudosigmenos1<-subset(a,a$dif<2 & a$num.predictors>1)
#modelomulticomtudosigmenos3<-subset(a,a$dif<2 & a$num.predictors>4)
#modelo_uni<-subset(a,a$num.predictors==1)

#htmlTable::htmlTable(modelo_uni)

#htmlTable::htmlTable(modelomulticomtudosig)

#htmlTable::htmlTable(modelomulticomtudosig)

#dataset

table(dataset$LVC)

## 
##   LVC Ñ_LVC 
##   635  5599

tabela$TRA_ESQUEM2<-0
tabela$TRA_ESQUEM2[tabela$TRA_ESQUEM=="4"]<-1
tabela$TRA_ESQUEM2<-as.factor(tabela$TRA_ESQUEM2)

tabela$TRA_ESQUEM3<-0
tabela$TRA_ESQUEM3[tabela$TRA_ESQUEM=="8"]<-1
tabela$TRA_ESQUEM3<-as.factor(tabela$TRA_ESQUEM3)

tabela$TRA_ESQUEM4<-0
tabela$TRA_ESQUEM4[tabela$TRA_ESQUEM=="10"]<-1
tabela$TRA_ESQUEM4<-as.factor(tabela$TRA_ESQUEM4)


tabela$TRA_ESQUEM5<-0
tabela$TRA_ESQUEM5[tabela$TRA_ESQUEM=="3"]<-1
tabela$TRA_ESQUEM5<-as.factor(tabela$TRA_ESQUEM5)


#a<-dataset$diff[dataset$recorrencia3==1]
#edu::resumosnumericos(a)
#a<-as.data.frame(a)
#edu::resumosnumericos(a)
#a<-dataset$diff[dataset$recorrencia2==1]
#a<-as.data.frame(a)
#table(dataset$recorrencia3)
#table(dataset$recorrencia3)
#a<-dataset$diff[dataset$recorrencia2==1]
#edu::resumosnumericos(a)
#a<-dataset$diff[dataset$recorrencia2==1]
#a<-as.data.frame(a)
#edu::resumosnumericos(a)

#modelos univariado logistico
 require(survminer)
#edu::modelo_uni(tabela)
#modelo multi logistico logit
sjPlot::tab_model(glm(recorrencia2 ~ CS_SEXO+PCRUZ+LVC+TRA_ESQUEM3+TRA_ESQUEM4,data=dataset,family =binomial(logit)))

	recorrencia2
Predictors	Odds Ratios	CI	p
(Intercept)	0.37	0.27 – 0.50	<0.001
CS_SEXO [M]	1.27	1.03 – 1.57	0.029
PCRUZ [2]	1.23	0.88 – 1.70	0.220
PCRUZ [3]	1.33	0.99 – 1.80	0.061
PCRUZ [4]	1.61	1.26 – 2.06	<0.001
PCRUZ [5]	1.55	1.10 – 2.18	0.011
PCRUZ [6]	1.38	0.79 – 2.29	0.235
LVC [Ñ_LVC]	0.14	0.11 – 0.16	<0.001
TRA_ESQUEM3 [1]	8.19	1.07 – 49.79	0.022
TRA_ESQUEM4 [1]	7.11	3.32 – 15.19	<0.001
Observations	6234
R2 Tjur	0.100

#Modelo gee
#sjPlot::tab_model(gee::gee(recorrencia2 ~ CS_SEXO +PCRUZ+LVC+TRA_ESQUEM3+TRA_ESQUEM4,data=dataset,id=study_ID))

#modelo de cox
sjPlot::tab_model(coxph(Surv(diff, recorrencia2) ~ escola2 + TRA_ESQUEM2 + TRA_ESQUEM3,data=dataset,id=study_ID))

	Surv(diff, recorrencia2)
Predictors	Estimates	CI	p
escola2 [1]	1.23	1.06 – 1.44	0.008
TRA_ESQUEM2 [1]	0.50	0.26 – 0.97	0.041
TRA_ESQUEM3 [1]	4.21	1.05 – 16.94	0.043
Observations	6234
R2 Nagelkerke	0.021

.Survfit <- survfit(Surv(diff, recorrencia2) ~ escola2 + TRA_ESQUEM2 + TRA_ESQUEM3 , error="greenwood",data=dataset,id=study_ID)

plot(.Survfit)

#3 é 8
#4 é 10