library(msm)
library(mstate)
library(dynpred)
library(mstate)
library(mstate)
library(dynpred)
library(survival)
library(survminer)
library(foreign)
library(ggplot2)
library(rms)
library(pec)
library(xfun)
library(Hmisc)
library(prodlim)
library(mstate)
library(mstate) 
library(dynpred)
pro<-read.dta("C:/Users/hp/Desktop/working data/merged_all_R18.dta")
head(pro)
pro$sex <- factor(pro$sex,labels=c("Male","Female"))
pro$cd4 <- factor(pro$cd4,labels=c("<=200","201-350", "351-500", ">=501"))
pro$tb <- factor(pro$tb,labels=c("Negative","Positive"))
pro$oi <- factor(pro$oi,labels=c("No OI","Yes"))
tmat1 <- transMat(x = list(c(2, 3, 5, 6), c(4, 5, 6),  c(4, 5, 6), c(5, 6), c(), c()), names = c("1a", "1b", "1c", "1d", "1e", "1f"))
print(tmat1)
pro <- expand.covs(pro, covs=("sex"))
pro<-data.frame(pro)
fix(pro)
pro <- expand.covs(pro, covs="sex")
covs<-c("sex")
pro <- expand.covs(pro, covs, longnames = FALSE)
g$x1 <- factor(tg$x1,labels=c("male","female"))
# expanded covariates
pro<-data.frame(pro)
fix(pro)
pro$sex <- factor(pro$sex,labels=c("male","female"))
pro <- expand.covs(pro, covs="sex")
pro1<-pro[pro$sex=="Male",]
head(pro)
table(pro$sex)
c0 <- coxph(Surv(Tstart,Tstop,status) ~ strata(trans), data=pro, method="breslow")
c1<-coxph(Surv(Tstart,Tstop,status) ~ wt0+ age + tb + sex +  patie_residence + educa_level +   vrlmmm_1000 +   cd4 + strata(trans), data=pro, method = "breslow")
summary(pmodel11)
revpro <- pro[rev(1:nrow(pro)), c("record_id", "to", "Tstop", "status")]
revpro <- revpro[!duplicated(revpro$record_id), ]
tail(revpro)


n <- nrow(pro)
tmat <- attr(pro, "trans")
nd <-pro[rep(1,6),]
nd$trans <- 1:6
nd$strata <- 1:6
attr(nd, "trans") <- tmat
msf0 <- msfit(c0, newdata=nd, variance=FALSE, trans=tmat)
msf0 <- msfit(c1, newdata=nd, variance=FALSE, trans=tmat)

pt0 <- probtrans(msf0, predt=0, variance=FALSE)












###pro<- merge(pro, csl[, c("id", "sex", "age3", "wt0", "cd40_c", "oi_cat", "anemia", "TB_history", "followuptime")], by="id")
class(pro2) <- c("pro2","data.frame")
class(csl) <- c("csl","data.frame")
K <- 4 # number of states in this multi-state model
# Expand covariates
table(pro$sex)
c0 <- coxph(Surv(Tstart,Tstop,status) ~ strata(trans), data=pro, id=record_id, method="breslow")
c1 <- coxph(Surv(start,stop,status) ~ wt0 + age2 +  Oi_cat + TB_history + strata(trans), data=pro2, id=id, method="breslow")
tmat1 <- transMat(x = list(c(2,3,4), c(1,3,4), c(1,2,4), c(1,2,3)), names=c("1", "2", "3", "4"))
tmat1
pmodel11=coxph(Surv(Tstart,Tstop,status) ~ wt0+ age + tb + sex +  patie_residence + educa_level +   vrlmmm_1000 +   cd4 + strata(trans), data=pro, method = "breslow")
summary(pmodel11)
pmodel0=coxph(Surv(followuptime,status)~1,data=pro,x=TRUE,y=TRUE)
perror1=pec(list(Cox=pmodel0),Hist(month,status)~wt0+age, data=pro)
pmodel=coxph(Surv(followuptime,status)~wt0+ age + tb + sex + anemia +  patie_residence + educa_level +   vrlmmm_1000 +   cd4 + strata(trans), sub=c(trans==1),data=pro,x=TRUE,y=TRUE)
summary(pmodel)
perror=pec(list(Cox=pmodel),Hist(month,status)~wt0+ age + tb + sex +patie_residence + educa_level + vrlmmm_1000 + cd4, data=pro)
crps(perror)
plot(perror)
####Estimation
table(pro$status)
fit1=coxph(Surv(start,stop,status==0)~strata(trans),  data=pro, method = "breslow")
summary(fit1)
pro<-read.dta("C:/Users/hp/Desktop/working data/work_23_2.dta")
head(pro)
table(pro$substitute2)
tmat1 <- transMat(x = list(c(2, 3, 4, 5, 6), c(1,3,4, 5, 6),  c(1,2,4, 5, 6), c(1,2,3,5, 6), c(), c()), names = c("1a", "1b", "1c", "1d", "1e", "1f"))
tmat1
pro<-
fit2=coxph(Surv(followuptime, status)~ age + wt0 + strata(reg), sub=c(substitute=="d4T"), data=pro, method = "breslow")
fit1=coxph(Surv(followuptime, status)~ wt0+ age + tb + sex + patie_residence + cd4 + strata(trans),  data=pro, method = "breslow")
c0 <- coxph(Surv(Tstart,Tstop,status) ~ strata(trans), data=pro, id=record_id, method="breslow")
c1 <- coxph(Surv(Tstart,Tstop,status) ~ sex + strata(trans), data=pro, id=record_id, method="breslow")
fit1=coxph(Surv(followuptime,status==0) ~ wt0+ age + tb + sex +  cd4 + strata(trans), data=pro, method = "breslow")
summary(fit1)
fit1=coxph(Surv(followuptime,status) ~ age + sex+ wt0+ cd4 + strata(trans), data=pro, sub=c(substitute2=="NNRTI"), method = "breslow")
fit1=coxph(Surv(followuptime,status==1) ~ wt0+ age + sex +  cd4 + strata(trans), data=pro,  method = "breslow")
fit1=coxph(Surv(followuptime,status) ~ wt0+ age + sex +  cd4 + strata(trans), data=pro, method = "breslow")
fit1=coxph(Surv(followuptime,status) ~ age + sex+ wt0+ cd4 + strata(trans), data=pro, method = "breslow")

head(pro)

nd <- pro[rep(1,6),]
nd$trans <- 1:6
nd$strata <- 1:6
table(pro$reg)
table(pro$oi)
head(pro)
+ age + sex + tb + cd4+
# Prediction
  par(mfrow=c(3,2))
msf0 <- msfit(object=fit1, newdata = nd,  trans=tmat1)
pt0 <- probtrans(msf0, predt=0, method = "greenwood",  direction="forward")
pt0 
plot(pt0, type="single", lwd=2,  from=1, ord=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6)) 
title(main = "From 1")
plot(pt0, type="single", lwd=2,  from=2, ord=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6))
title(main = "From 2")
plot(pt0, type="single", lwd=2,  from=3, ord=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6))
title(main = "From 3")
plot(pt0, type="single", lwd=2,  from=4, ord=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6))
title(main = "From 4")
plot(pt0, type="single", lwd=2,  from=5, ord=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6))
title(main = "From 5")
plot(pt0, type="single", lwd=2,  from=6, ord=c(1,2,3,4,5,6), col=c(1,2,3,4,5,6))
title(main = "From 6")
plot(pt0, type="single", lwd=2,  from=2, ord=c(1,3,2,4, 5,6))
plot(pt0, type="single", lwd=2,  from=3, ord=c(1,3,2,4,5,6))
plot(pt0, type="single", lwd=2,  from=4, ord=c(1,3,2, 4))
pt0 <- probtrans(msf0, predt = 10,direction="forward",method = "greenwood")
plot(pt0, type="single", lwd=2,  from=1, ord=c(1,3,2, 4))
title(main = "From 3")
perror=pec(list(Cox=pmodel),Hist(time,status)~1,data=dat)
brier <- pec(list(pt0),data=pro,formula=Surv(followuptime,status!=0)~strata(trans))
print(predct)
plot(predct)
head(pro)
Models <- list("fit1"=coxph(Surv(followuptime, status)~strata(trans),  data=pro, method = "breslow", x=TRUE,y=TRUE),
               "fit2"=coxph(Surv(followuptime, status)~age3 + wt0 + Oi_cat + cd40_c,data=pro, method = "breslow", x=TRUE,y=TRUE))
PredError <- pec(object=Models, formula=Surv(followuptime,status)~age3 + wt0 + Oi_cat + cd40_c,
                 data=pro,   exact=TRUE, cens.model="marginal",  splitMethod="none", B=0, verbose=TRUE)
print(PredError, followuptimes=seq(0,30,80))
summary(PredError)
plot(PredError,xlim=c(0, 90))

brier <- pec(list(nullmodel,pmodel1, pmodel2),data=pro, formula=Surv(followuptime,status!=0)~age + wt0 )
print(brier)
plot(brier)
plot(brier,xlim=c(30, 100))

nullmodel=prodlim(Hist(followuptime,status)~1,data=pro)
pmodel1=coxph(Surv(followuptime,status)~age3 + wt0 + Oi_cat + cd40_c,data=pro,x=TRUE,y=TRUE)
pmodel2=coxph(Surv(followuptime,status)~age3 ,data=pro,x=TRUE,y=TRUE)
perror=pec(list(Cox1=pmodel1,Cox2=pmodel2),Hist(followuptime,status)~age3 + wt0 + Oi_cat + cd40_c,data=pro,reference=TRUE)
R2(perror, times=seq(0,1,.1),reference=1)
plot(perror,xlim=c(0, 90))

# assessing parametric survival models in learn/validation setting
learndat <- read.dta("C:/Users/hp/Desktop/data22/data22/learningdata.dta")
testdat <- read.dta("C:/Users/hp/Desktop/data22/data22/testdata.dta")
library(rms)
f1 <- psm(Surv(followuptime,status)~age3 + wt0 + Oi_cat + cd40_c,data=learndat)
f2 <- psm(Surv(followuptime,status)~age3,data=learndat)
pf <- pec(list(f1,f2),formula=Surv(followuptime,status)~ age3 + wt0 + Oi_cat + cd40_c,data=testdat,maxtime=84)
plot(pf)
summary(pf) 
pf <- pec(list(f1,f2),formula=Surv(followuptime,status)~ age3 + wt0 + Oi_cat + cd40_c,data=testdat,maxtime=84)

Models <- list("fit1"=coxph(Surv(followuptime,status)~age3 + wt0 + Oi_cat + cd40_c,data=learndat, method = "breslow", x=TRUE,y=TRUE),
               "fit2"=coxph(Surv(followuptime,status)~1,data=learndat, method = "breslow", x=TRUE,y=TRUE))
PredError <- pec(object=Models, formula=Surv(followuptime,status)~ age3 + wt0 + Oi_cat + cd40_c,data=testdat,maxtime=84,
                 exact=TRUE, cens.model="marginal",  splitMethod="none", B=0, verbose=TRUE)
print(PredError, followuptimes=seq(0,30,80))
summary(PredError)
plot(PredError,xlim=c(0, 90))

##################################################################
library(riskRegression)
        data(Melanoma)
 Melanoma[1:5,1:5]

 Melanoma$id<-1:nrow(Melanoma)
  set.seed(123)
  ind.split<-sample(1:nrow(Melanoma), round(nrow(Melanoma)*4/5), replace = F)
  dftrain<-Melanoma[ind.split,]
  dftest<-Melanoma[-ind.split,]
  csc <- CSC(Hist(time, status)~age+thick+ulcer,  data=dftrain)
  fgr <- FGR(Hist(time,status)~age+thick+ulcer,  data=dftrain,cause=1)
  fgr.full<-FGR(Hist(time,status)~age+thick+ici+epicel+ulcer+invasion+sex, data=dftrain,cause=1)
  pred.csc<-predict(csc, newdata = dftest[1:5,], time = 1000:4000, cause = 1)
  autoplot(pred.csc)
predictRisk(csc,newdata = dftest[20,],times = 1500,cause = 1)
 
  predictRisk(fgr,newdata = dftest[20,],times = 1500)
   library(survival)
  score<-Score(list("Cause-specific Cox"=csc,"Fine-Gray"=fgr,  "Full"=fgr.full),
                 formula = Hist(time,status)~1,
                 data=dftest,times = seq(1900,4000,100),
                 plots = "calibration",
                 summary = "risks")
  dev.new(width=5,height=4)
  plotCalibration(score,times = 2000)
  
  ######################################################
  pro$id<-1:nrow(pro)
  set.seed(123)
  ind.split<-sample(1:nrow(pro),  round(nrow(pro)*4/5), replace = F)
  dftrain<-pro[ind.split,]
  dftest<-pro[-ind.split,]
  head(dftrain)