R/Markov_noipd.R
In HealthEconomicsHackathon/ceatutorial: What the Package Does (One Line, Title Case)
Defines functions
Markov_noipd
##########################################################################################
##========================================================================================
## START OF MARKOV_NOIPD FUNCTION
##========================================================================================
##########################################################################################
##========================================================================================
## THIS FUNCTION RETURNS THE STATE OCCUPANCY PROBABILITIES (PROBABILITIES OVER TIME OF
## BEING IN EACH STATE) BASED ON THE MODEL COEFFICIENTS OR MORTALITY RATES GIVEN AS INPUTS
## MODIFIBLE ARGUMENTS:
## ntrans NUMBER OF TRANSITIONS IN THE MODELLING
## type SOURCE OF INPUT. EITHER "parametric model" OR "mortality hazard".
## REQUIRED FOR EACH TRANSITION
## morhaz APPLICABLE IF TYPE="mortality hazard".
## ncovs NUMBER OF COVARIATE PARAMETERS IN THE MODEL (EXCLUDING INTERCEPT) FOR EACH
## TRANSITION I.E. ONE FOR EACH BINARY AND CONTINUOUS VARIABLE AND K-1 FOR
## EACH CATEGORICAL VARIABLE, WHERE K= NO OF CATEGORIES
## covs VARIABLE NAMES FOR THE COVARIATES
## covcoeff COEFFICIENTS RELATING TO EACH OF THE COVARIATES FOR EACH OF THE TRANSITIONS
## coveval VALUE AT WHICH TO EVALUATE EACH COVARIATE
## dist DISTRIBUTION TO USE FOR EACH TRANSITION. OPTIONS ARE:
## wei FOR WEIBULL, exp FOR EXPONENTIAL, gom FOR GOMPERTZ,
## logl FOR LOGLOGISTIC,logn FOR LOGNORMAL AND gam FOR GENERALISED GAMMA
## logscale COEFFICIENT FOR LOG(SCALE) FOR EACH OF THE TRANSITIONS - APPLICABLE TO
## WEIBULL, EXPONENTIAL, LOGLOGISTIC, LOGNORMAL
## logshape COEFFICIENT FOR LOG(SHAPE) FOR EACH OF THE TRANSITIONS - APPLICABLE TO
## WEIBULL, LOGLOGISTIC,LOGNORMAL
## rate COEFFICIENT FOR RATE - APPLICABLE TO GOMPERTZ
## loglevel COEFFICIENT FOR LOG(LEVEL) - APPLICABLE TO GOMPERTZ
## mu COEFFICIENT FOR MU - APPLICABLE TO GENERALISED GAMMA
## sigma COEFFICIENT FOR SIGMA - APPLICABLE TO GENERALISED GAMMA
## Q COEFFICIENT FOR Q - APPLICABLE TO GENERALISED GAMMA
## alt1end THIS SHOULD BE SET TO TRUE IF AN ADJUSTMENT IS REQUIRED TO TIMEPOINTS AT
## THE END OF SEQUENCE (SEE README FILE). OTHERWISE IT SHOULD BE FALSE
## alt2beg THIS SHOULD BE SET TO TRUE IF AN ADJUSTMENT IS REQUIRED TO TIMEPOINTS AT
## THE BEGINNING OF THE SEQUENCE (SEE README FILE).OTHERWISE IT SHOULD BE
## FALSE
## ttend APPLICABLE IF alt1end and alt2beg are FALSE. TIMEPOINT AT WHICH TO END
## EXTRAPOLATION
## ttinc APPLICABLE IF alt1end and alt2beg are FALSE. GAP BETWEEN TIMEPOINTS
## tt1end APPLICABLE IF alt1end = TRUE. THE ENDPOINT OF TIME BEFORE ADJUSTMENT
## TO THE LAST PART OF THE TIME SEQUENCE IS REQUIRED
## tt1inc APPLICABLE IF alt1end = TRUE. THE GAP BETWEEN TIMEPOINTS TO USE
## BEFORE ADJUSTMENT TO THE LAST PART OF THE TIME SEQUENCE
## tt1end2 APPLICABLE IF alt1end = TRUE. THE TIME HORIZON I.E. THE ENDPOINT OF
## THE ADJUSTMENT TO THE LAST PART OF THE TIME SEQUENCE
## tt1inc2 APPLICABLE IF alt1end = TRUE. THE GAP BETWEEN TIMEPOINTS TO USE FOR
## THE ADJUSTMENT TO THE LAST PART OF THE TIME SEQUENCE (WHICH SHOULD BE
## SMALLER THAN tt1inc)
## tt2beg APPLICABLE IF alt2beg = TRUE. THE TIMEPOINTS TO BE ADDED TO THE
## BEGINNING OF THE TIME SEQUENCE
## tt2start APPLICABLE IF alt2beg = TRUE. POINT AT WHICH TO START THE NEXT SECTION OF
## THE TIME SEQUENCE (AFTER THE ADJUSTMENT OF ADDITIONAL TIMEPOINTS TO THE
## BEGINNING OF THE SEQUENCE)
## tt2inc APPLICABLE IF alt2beg = TRUE and alt1end = FALSE. THE GAP BETWEEN
## TIMEPOINTS TO USE FOR THE LAST PART OF THE TIME SEQUENCE
## tt2end APPLICABLE IF alt2beg = TRUE and alt1end = FALSE. THE TIME HORIZON
## I.E. THE FINAL POINT OF THE TIME SEQUENCE
## data DATASET TO USE FOR MODELLING
## trans TRANSITION MATRIX
##========================================================================================
## NB. THIS FUNCTION USES THE SAME PARAMETERISATIONS AS phreg AND aftreg IN THE eha PACKAGE
## AND flexsurvreg IN THE flexsurv PACKAGE
Markov_noipd<-function(ntrans=3,
type=c("parametric model","parametric model", "mortality hazard"),
morhaz=c(0,0, 0.05),
ncovs=c(1,1,2),
covs=rbind("covariate1", "covariate1",c("covariate1", "covariate2")),
covcoeff=rbind(-0.537, -0.402, c(0.487, -0.687)),
coveval=rbind(0,0,c(0,1)),
dist=cbind("wei", "wei", "wei"),
logscale=c(1.289, 3.997, 0.174 ),
logshape=c(0.441, -0.236, -0.524 ),
rate=c(0.474, -0.487),
loglevel=c( -2.187, -2.825),
mu=c(1.300, 3.160, 1.370),
sigma=c(0.623, 5.780, 0.806),
Q=c(1.050, -3.280, 1.160 ),
alt1end=FALSE, alt2beg=FALSE,
ttend=15, ttinc=1/12,
tt1end=118/12, tt1inc=1/12,tt1end2=15, tt1inc2=0.007,
tt2beg=c(0,0.00001,0.005), tt2start=1/12, tt2inc=1/12,tt2end=15,
data=msmcancer, trans=tmat){
#### create the timepoints
if(alt1end==FALSE & alt2beg==FALSE){
tt<-seq(0,ttend,ttinc)
}
if(alt1end==TRUE & alt2beg==FALSE){
tt<-c(seq(0,tt1end,tt1inc), seq(tt1end+tt1inc2, tt1end2,tt1inc2))
}
if(alt2beg==TRUE & alt1end==FALSE){
tt<-c(tt2beg, seq(tt2start, tt2end,tt2inc))
}
if(alt2beg==TRUE & alt1end==TRUE){
tt<-c(tt2beg, seq(tt2start, tt1end,tt1inc), seq(tt1end+tt1inc2, tt1end2,tt1inc2))
}
#### set up required lists
covcoeffs<-vector("list", ntrans)
lp<-vector("list", ntrans)
x<-vector("list", ntrans)
cumHaz<-vector("list", ntrans)
gamma<-vector("list", ntrans)
mu2<-vector("list", ntrans)
z<-vector("list", ntrans)
u<-vector("list", ntrans)
for (i in 1:ntrans) {
#### calculations based on mortality hazards
if (type[i]=="mortality hazard"){
cumHaz[[i]] <-morhaz[i]*tt
}
#### calculations based on coefficients from each transition
x[[i]]<-coveval[i,1:ncovs[i]]
covcoeffs[[i]]<-covcoeff[i,1:ncovs[i]]
lp[[i]]<-sum(covcoeffs[[i]]* x[[i]] )
#### cumulative hazards for each transition - distribution specific
if (type[i]== "parametric model" & dist[i]=="wei") {
cumHaz[[i]]<-exp(-exp(logshape[i])*logscale[i]+lp[[i]])*tt^exp(logshape[i])
}
if (type[i]== "parametric model" & dist[i]=="exp") {
cumHaz[[i]]<-exp(-logscale[i]+lp[[i]])*tt
}
if (type[i]== "parametric model" & dist[i]=="gom") {
cumHaz[[i]]<-exp(loglevel[i]+lp[[i]])*(1/rate[i])*(exp(rate[i]*tt) -1)
}
if (type[i]== "parametric model" & dist[i]=="logl") {
cumHaz[[i]]<- -log(1/(1+(exp(lp[[i]]-logscale[i])*tt)^
(1/exp(-logshape[i]))))
}
if (type[i]== "parametric model" & dist[i]=="logn") {
cumHaz[[i]]<- -log(1-pnorm((log(tt)-(logscale[i]-lp[[i]]))/
(exp(-logshape[i]))))
}
if (type[i]== "parametric model" & dist[i]=="gam") {
gamma[[i]]<-(abs(Q[i]))^(-2)
mu2[[i]]<- mu[i] +lp[[i]]
z[[i]]<-rep(0,length(tt))
z[[i]]<- sign(Q[i])*((log(tt)-mu2[[i]])/sigma[i])
u[[i]]<-gamma[[i]]*exp((abs(Q[i]))*z[[i]])
if(Q[i]>0){
cumHaz[[i]]<--log(1-pgamma(u[[i]],gamma[[i]]))
}
if(Q[i]==0){
cumHaz[[i]]<--log(1-pnorm(z[[i]]))
}
if(Q[i]<0){
cumHaz[[i]]<--log(pgamma(u[[i]],gamma[[i]]))
}
}
}
Haz<-unlist(cumHaz)
newtrans<-rep(1:ntrans,each=length(tt))
time<-rep(tt,ntrans)
Haz<-cbind(time=as.vector(time),Haz=as.vector(Haz),trans=as.vector(newtrans))
Haz<-as.data.frame(Haz)
##### state occupancy probabilities
msf<-list(Haz,trans)
msf$trans<-trans
msf$Haz<-Haz
class(msf)<-"msfit"
stateprobs <- probtrans(msf, predt = 0,variance=FALSE)
return(stateprobs)
}
##########################################################################################
##========================================================================================
## END OF MARKOV_NOIPD FUNCTION
##========================================================================================
##########################################################################################
HealthEconomicsHackathon/ceatutorial documentation built on Nov. 7, 2019, 12:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Markov_noipd
##########################################################################################
##========================================================================================
## START OF MARKOV_NOIPD FUNCTION
##========================================================================================
##########################################################################################
##========================================================================================
## THIS FUNCTION RETURNS THE STATE OCCUPANCY PROBABILITIES (PROBABILITIES OVER TIME OF
## BEING IN EACH STATE) BASED ON THE MODEL COEFFICIENTS OR MORTALITY RATES GIVEN AS INPUTS
## MODIFIBLE ARGUMENTS:
## ntrans NUMBER OF TRANSITIONS IN THE MODELLING
## type SOURCE OF INPUT. EITHER "parametric model" OR "mortality hazard".
## REQUIRED FOR EACH TRANSITION
## morhaz APPLICABLE IF TYPE="mortality hazard".
## ncovs NUMBER OF COVARIATE PARAMETERS IN THE MODEL (EXCLUDING INTERCEPT) FOR EACH
## TRANSITION I.E. ONE FOR EACH BINARY AND CONTINUOUS VARIABLE AND K-1 FOR
## EACH CATEGORICAL VARIABLE, WHERE K= NO OF CATEGORIES
## covs VARIABLE NAMES FOR THE COVARIATES
## covcoeff COEFFICIENTS RELATING TO EACH OF THE COVARIATES FOR EACH OF THE TRANSITIONS
## coveval VALUE AT WHICH TO EVALUATE EACH COVARIATE
## dist DISTRIBUTION TO USE FOR EACH TRANSITION. OPTIONS ARE:
## wei FOR WEIBULL, exp FOR EXPONENTIAL, gom FOR GOMPERTZ,
## logl FOR LOGLOGISTIC,logn FOR LOGNORMAL AND gam FOR GENERALISED GAMMA
## logscale COEFFICIENT FOR LOG(SCALE) FOR EACH OF THE TRANSITIONS - APPLICABLE TO
## WEIBULL, EXPONENTIAL, LOGLOGISTIC, LOGNORMAL
## logshape COEFFICIENT FOR LOG(SHAPE) FOR EACH OF THE TRANSITIONS - APPLICABLE TO
## WEIBULL, LOGLOGISTIC,LOGNORMAL
## rate COEFFICIENT FOR RATE - APPLICABLE TO GOMPERTZ
## loglevel COEFFICIENT FOR LOG(LEVEL) - APPLICABLE TO GOMPERTZ
## mu COEFFICIENT FOR MU - APPLICABLE TO GENERALISED GAMMA
## sigma COEFFICIENT FOR SIGMA - APPLICABLE TO GENERALISED GAMMA
## Q COEFFICIENT FOR Q - APPLICABLE TO GENERALISED GAMMA
## alt1end THIS SHOULD BE SET TO TRUE IF AN ADJUSTMENT IS REQUIRED TO TIMEPOINTS AT
## THE END OF SEQUENCE (SEE README FILE). OTHERWISE IT SHOULD BE FALSE
## alt2beg THIS SHOULD BE SET TO TRUE IF AN ADJUSTMENT IS REQUIRED TO TIMEPOINTS AT
## THE BEGINNING OF THE SEQUENCE (SEE README FILE).OTHERWISE IT SHOULD BE
## FALSE
## ttend APPLICABLE IF alt1end and alt2beg are FALSE. TIMEPOINT AT WHICH TO END
## EXTRAPOLATION
## ttinc APPLICABLE IF alt1end and alt2beg are FALSE. GAP BETWEEN TIMEPOINTS
## tt1end APPLICABLE IF alt1end = TRUE. THE ENDPOINT OF TIME BEFORE ADJUSTMENT
## TO THE LAST PART OF THE TIME SEQUENCE IS REQUIRED
## tt1inc APPLICABLE IF alt1end = TRUE. THE GAP BETWEEN TIMEPOINTS TO USE
## BEFORE ADJUSTMENT TO THE LAST PART OF THE TIME SEQUENCE
## tt1end2 APPLICABLE IF alt1end = TRUE. THE TIME HORIZON I.E. THE ENDPOINT OF
## THE ADJUSTMENT TO THE LAST PART OF THE TIME SEQUENCE
## tt1inc2 APPLICABLE IF alt1end = TRUE. THE GAP BETWEEN TIMEPOINTS TO USE FOR
## THE ADJUSTMENT TO THE LAST PART OF THE TIME SEQUENCE (WHICH SHOULD BE
## SMALLER THAN tt1inc)
## tt2beg APPLICABLE IF alt2beg = TRUE. THE TIMEPOINTS TO BE ADDED TO THE
## BEGINNING OF THE TIME SEQUENCE
## tt2start APPLICABLE IF alt2beg = TRUE. POINT AT WHICH TO START THE NEXT SECTION OF
## THE TIME SEQUENCE (AFTER THE ADJUSTMENT OF ADDITIONAL TIMEPOINTS TO THE
## BEGINNING OF THE SEQUENCE)
## tt2inc APPLICABLE IF alt2beg = TRUE and alt1end = FALSE. THE GAP BETWEEN
## TIMEPOINTS TO USE FOR THE LAST PART OF THE TIME SEQUENCE
## tt2end APPLICABLE IF alt2beg = TRUE and alt1end = FALSE. THE TIME HORIZON
## I.E. THE FINAL POINT OF THE TIME SEQUENCE
## data DATASET TO USE FOR MODELLING
## trans TRANSITION MATRIX
##========================================================================================
## NB. THIS FUNCTION USES THE SAME PARAMETERISATIONS AS phreg AND aftreg IN THE eha PACKAGE
## AND flexsurvreg IN THE flexsurv PACKAGE
Markov_noipd<-function(ntrans=3,
type=c("parametric model","parametric model", "mortality hazard"),
morhaz=c(0,0, 0.05),
ncovs=c(1,1,2),
covs=rbind("covariate1", "covariate1",c("covariate1", "covariate2")),
covcoeff=rbind(-0.537, -0.402, c(0.487, -0.687)),
coveval=rbind(0,0,c(0,1)),
dist=cbind("wei", "wei", "wei"),
logscale=c(1.289, 3.997, 0.174 ),
logshape=c(0.441, -0.236, -0.524 ),
rate=c(0.474, -0.487),
loglevel=c( -2.187, -2.825),
mu=c(1.300, 3.160, 1.370),
sigma=c(0.623, 5.780, 0.806),
Q=c(1.050, -3.280, 1.160 ),
alt1end=FALSE, alt2beg=FALSE,
ttend=15, ttinc=1/12,
tt1end=118/12, tt1inc=1/12,tt1end2=15, tt1inc2=0.007,
tt2beg=c(0,0.00001,0.005), tt2start=1/12, tt2inc=1/12,tt2end=15,
data=msmcancer, trans=tmat){
#### create the timepoints
if(alt1end==FALSE & alt2beg==FALSE){
tt<-seq(0,ttend,ttinc)
}
if(alt1end==TRUE & alt2beg==FALSE){
tt<-c(seq(0,tt1end,tt1inc), seq(tt1end+tt1inc2, tt1end2,tt1inc2))
}
if(alt2beg==TRUE & alt1end==FALSE){
tt<-c(tt2beg, seq(tt2start, tt2end,tt2inc))
}
if(alt2beg==TRUE & alt1end==TRUE){
tt<-c(tt2beg, seq(tt2start, tt1end,tt1inc), seq(tt1end+tt1inc2, tt1end2,tt1inc2))
}
#### set up required lists
covcoeffs<-vector("list", ntrans)
lp<-vector("list", ntrans)
x<-vector("list", ntrans)
cumHaz<-vector("list", ntrans)
gamma<-vector("list", ntrans)
mu2<-vector("list", ntrans)
z<-vector("list", ntrans)
u<-vector("list", ntrans)
for (i in 1:ntrans) {
#### calculations based on mortality hazards
if (type[i]=="mortality hazard"){
cumHaz[[i]] <-morhaz[i]*tt
}
#### calculations based on coefficients from each transition
x[[i]]<-coveval[i,1:ncovs[i]]
covcoeffs[[i]]<-covcoeff[i,1:ncovs[i]]
lp[[i]]<-sum(covcoeffs[[i]]* x[[i]] )
#### cumulative hazards for each transition - distribution specific
if (type[i]== "parametric model" & dist[i]=="wei") {
cumHaz[[i]]<-exp(-exp(logshape[i])*logscale[i]+lp[[i]])*tt^exp(logshape[i])
}
if (type[i]== "parametric model" & dist[i]=="exp") {
cumHaz[[i]]<-exp(-logscale[i]+lp[[i]])*tt
}
if (type[i]== "parametric model" & dist[i]=="gom") {
cumHaz[[i]]<-exp(loglevel[i]+lp[[i]])*(1/rate[i])*(exp(rate[i]*tt) -1)
}
if (type[i]== "parametric model" & dist[i]=="logl") {
cumHaz[[i]]<- -log(1/(1+(exp(lp[[i]]-logscale[i])*tt)^
(1/exp(-logshape[i]))))
}
if (type[i]== "parametric model" & dist[i]=="logn") {
cumHaz[[i]]<- -log(1-pnorm((log(tt)-(logscale[i]-lp[[i]]))/
(exp(-logshape[i]))))
}
if (type[i]== "parametric model" & dist[i]=="gam") {
gamma[[i]]<-(abs(Q[i]))^(-2)
mu2[[i]]<- mu[i] +lp[[i]]
z[[i]]<-rep(0,length(tt))
z[[i]]<- sign(Q[i])*((log(tt)-mu2[[i]])/sigma[i])
u[[i]]<-gamma[[i]]*exp((abs(Q[i]))*z[[i]])
if(Q[i]>0){
cumHaz[[i]]<--log(1-pgamma(u[[i]],gamma[[i]]))
}
if(Q[i]==0){
cumHaz[[i]]<--log(1-pnorm(z[[i]]))
}
if(Q[i]<0){
cumHaz[[i]]<--log(pgamma(u[[i]],gamma[[i]]))
}
}
}
Haz<-unlist(cumHaz)
newtrans<-rep(1:ntrans,each=length(tt))
time<-rep(tt,ntrans)
Haz<-cbind(time=as.vector(time),Haz=as.vector(Haz),trans=as.vector(newtrans))
Haz<-as.data.frame(Haz)
##### state occupancy probabilities
msf<-list(Haz,trans)
msf$trans<-trans
msf$Haz<-Haz
class(msf)<-"msfit"
stateprobs <- probtrans(msf, predt = 0,variance=FALSE)
return(stateprobs)
}
##########################################################################################
##========================================================================================
## END OF MARKOV_NOIPD FUNCTION
##========================================================================================
##########################################################################################
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.