R/DengueParmEst.R

Defines functions changetimeFSTP allmodelpredict plot.dengue summary.dengue allmodels CIturningpoint CIFinalsize estimTpointMA estimFSizeMA weight logistic5pm weibullm gompertzm SigmEmaxm logistic3pm Richardsm

Documented in allmodelpredict allmodels changetimeFSTP plot.dengue summary.dengue

###packages
require(nlme)
require(ggplot2)
require(NRAIA)
require(drc)



####Richards model funtions
Richardsm<-function(inc , time, start=NULL){
if(is.null(start)){
tt<-data.frame(inc,time)
gi<-stats::getInitial(cumsum(inc) ~ SSRichards(time, Asym, xmid, scal, lpow), data = tt)
c<-as.data.frame(gi)
start=list(alpha=c$gi[1],k=c$gi[3],gamma=c$gi[4],eta=c$gi[2])
start=start}
gnls.richards<-gnls(cumsum(inc) ~ alpha*((1+k*exp(-gamma*k*(time-eta)))^(-1/k)),start=start)
p<-summary(gnls.richards);
return(p)}



####3P logistic model 
logistic3pm <-function(inc, time, start=NULL){
if(is.null(start)){
tt<-data.frame(inc,time)
gi<-stats::getInitial(cumsum(inc) ~ SSlogis(time, Asym, xmid, scal), data = tt)
c<-as.data.frame(gi)
start=list(alpha=c$gi[1],gamma=1/c$gi[3],eta=c$gi[2])
start=start}
gnls.3plogistic<-gnls(cumsum(inc)~alpha/(1+exp(-gamma*(time-eta))),start=start)
p<-summary (gnls.3plogistic);
return(p)}


####Sigmoidal Emax Model
SigmEmaxm <-function(inc, time, start=NULL){
if(is.null(start)){
tt<-data.frame(inc,time)
gi<-drm(cumsum(inc) ~time, data=tt,fct=LL.4())
start=list(alpha=gi$fit$par[3],eta=gi$fit$par[4],beta=gi$fit$par[2],n=-gi$fit$par[1])
start=start}
gnls.Emax<-gnls(cumsum(inc)~beta + (time^n)*(alpha-beta)/(time^n + eta^n),start=start)
p<-summary (gnls.Emax);
return(p)}



####Gompertz model 4 parametrs
gompertzm <-function(inc, time, start=NULL){
if(is.null(start)){
tt<-data.frame(inc,time)
gi<-stats::getInitial(cumsum(inc) ~ SSgompertz(time, Asym, b2, b3), data = tt)
c<-as.data.frame(gi)
start=list(alpha=c$gi[1],eta=-log(c$gi[2])/log(c$gi[3]),beta=0, gamma=-log(c$gi[3]))
start=start}
gnls.gompertz<-gnls(cumsum(inc)~beta + (alpha-beta)*exp(-exp(-gamma*(time-eta))),start=start)
p<-summary (gnls.gompertz);
return(p)}



####Modelo de Weibull
weibullm <-function(inc, time, start=NULL){
if(is.null(start)){
tt<-data.frame(inc,time)
gi<-stats::getInitial(cumsum(inc) ~ SSweibull(time, Asym, Drop, lrc, pwr), data = tt)
c<-as.data.frame(gi)
start=list(alpha=c$gi[1],eta=exp(-c$gi[3]/c$gi[4]),beta=c$gi[1]-c$gi[2], k=c$gi[4])
start=start}
gnls.weibull<-gnls(cumsum(inc)~beta + (alpha-beta)*exp(-(time/eta)^-k),start=start)
p<-summary(gnls.weibull);
return(p)}


####5p logistic model
logistic5pm <-function(inc, time, start=NULL){
if(is.null(start)){
tt<-data.frame(inc,time)
gi<-drm(cumsum(inc) ~time, data=tt,fct=LL.5())
gi2<-drm(cumsum(inc) ~time, data=tt,fct=LL.4())
start=list(alpha=gi$fit$par[2],beta=gi$fit$par[3],g=gi$fit$par[5],eta=gi2$fit$par[4], k=-gi$fit$par[1])
start=start}
gnls.5plogistic<-gnls(cumsum(inc)~beta+ (alpha-beta)/(1+(2^(1/g)-1)*(time/eta)^-k)^g,start=start)
p<-summary (gnls.5plogistic);
return(p)}


#####model averaging weight function
weight<-function(inc , time, start=NULL){
aic<-AIC(Richardsm(inc , time, start[[1]]),logistic3pm(inc , time, start[[2]]),SigmEmaxm (inc , time, start[[3]]),
     gompertzm(inc , time, start[[4]]),weibullm(inc , time, start[[5]]),logistic5pm (inc , time, start[[6]]))
minAIC<-min(aic[,2])
deltaAIC<-aic[,2]-minAIC
w<-signif(exp(-deltaAIC/2)/sum(exp(-deltaAIC/2)),6)
names(w)<-c("Richards","3P Logistic","Sigmoidal Emax","Gompertz","Weibull","5P Logistic");
return(w)}



####Model averaging estimator Final size
estimFSizeMA<-function(inc , time, start=NULL){
alpha.MA<-weight(inc , time, start)[1]*coef(Richardsm(inc , time, start[[1]]))[1] + weight(inc , time, start)[2]*coef(logistic3pm(inc , time, start[[2]]))[1] + weight(inc , time, start)[3]*coef(SigmEmaxm (inc , time, start[[3]]))[1] +weight(inc , time, start)[4]*coef(gompertzm(inc , time, start[[4]]))[1] +weight(inc , time, start)[5]*coef(weibullm(inc , time, start[[5]]))[1] + weight(inc , time, start)[6]*coef(logistic5pm (inc , time, start[[6]]))[1]
stderror.alpha.MA<-weight(inc , time, start)[1]*sqrt(Richardsm(inc , time, start[[1]])$tTable[1, 2]^2 + (coef(Richardsm(inc , time, start[[1]]))[1]-alpha.MA)^2) + 
weight(inc , time, start)[2]*sqrt(logistic3pm(inc , time, start[[2]])$tTable[1, 2]^2 + (coef(logistic3pm(inc , time, start[[2]]))[1]-alpha.MA)^2) + 
weight(inc , time, start)[3]*sqrt(SigmEmaxm (inc , time, start[[3]])$tTable[1, 2]^2 + (coef(SigmEmaxm (inc , time, start[[3]]))[1]-alpha.MA)^2) + 
weight(inc , time, start)[4]*sqrt(gompertzm(inc , time, start[[4]])$tTable[1, 2]^2 + (coef(gompertzm(inc , time, start[[4]]))[1]-alpha.MA)^2) + 
weight(inc , time, start)[5]*sqrt(weibullm(inc , time, start[[5]])$tTable[1, 2]^2 + (coef(weibullm(inc , time, start[[5]]))[1]-alpha.MA)^2) + 
weight(inc , time, start)[6]*sqrt(logistic5pm (inc , time, start[[6]])$tTable[1, 2]^2 + (coef(logistic5pm (inc , time, start[[6]]))[1]-alpha.MA)^2)
low<-alpha.MA-qt(0.975,5)*stderror.alpha.MA
up<-alpha.MA+qt(0.975,5)*stderror.alpha.MA
names(alpha.MA)<-c("est.")
names(low)<-c("lower")
names(up)<-c("upper")
p<-c(low,alpha.MA,up);
return(p)}


####Model averaging estimator Turning point 

estimTpointMA<-function(inc , time, start=NULL){
eta.MA<-weight(inc , time, start)[1]*coef(Richardsm(inc , time, start[[1]]))[4] + weight(inc , time, start)[2]*coef(logistic3pm(inc , time, start[[2]]))[3] + weight(inc , time, start)[3]*coef(SigmEmaxm (inc , time, start[[3]]))[2] +weight(inc , time, start)[4]*coef(gompertzm(inc , time, start[[4]]))[2] +weight(inc , time, start)[5]*coef(weibullm(inc , time, start[[5]]))[2] + weight(inc , time, start)[6]*coef(logistic5pm (inc , time, start[[6]]))[4]
stderror.eta.MA<-weight(inc , time, start)[1]*sqrt(Richardsm(inc , time, start[[1]])$tTable[4, 2]^2 + (coef(Richardsm(inc , time, start[[1]]))[4]-eta.MA)^2) + 
weight(inc , time, start)[2]*sqrt(logistic3pm(inc , time, start[[2]])$tTable[3, 2]^2 + (coef(logistic3pm(inc , time, start[[2]]))[3]-eta.MA)^2) + 
weight(inc , time, start)[3]*sqrt(SigmEmaxm (inc , time, start[[3]])$tTable[2, 2]^2 + (coef(SigmEmaxm (inc , time, start[[3]]))[2]-eta.MA)^2) + 
weight(inc , time, start)[4]*sqrt(gompertzm(inc , time, start[[4]])$tTable[2, 2]^2 + (coef(gompertzm(inc , time, start[[4]]))[2]-eta.MA)^2) + 
weight(inc , time, start)[5]*sqrt(weibullm(inc , time, start[[5]])$tTable[2, 2]^2 + (coef(weibullm(inc , time, start[[5]]))[2]-eta.MA)^2) + 
weight(inc , time, start)[6]*sqrt(logistic5pm (inc , time, start[[6]])$tTable[4, 2]^2 + (coef(logistic5pm (inc , time, start[[6]]))[4]-eta.MA)^2)
low<-eta.MA-qt(0.975,5)*stderror.eta.MA
up<-eta.MA+qt(0.975,5)*stderror.eta.MA
names(eta.MA)<-c("est.")
names(low)<-c("lower")
names(up)<-c("upper")
p<-c(low,eta.MA,up);
return(p)
}


#### final size Confidence interval function
CIFinalsize<-function(inc , time, start=NULL, model){
switch(model,
       Richards = {intervals(Richardsm(inc , time, start))$coef[1,]},
       logistic3P = {intervals(logistic3pm(inc , time, start))$coef[1,]},
       SigmEmax = {intervals(SigmEmaxm (inc , time, start))$coef[1,]},
       Gompertz = {intervals(gompertzm(inc , time, start))$coef[1,]},
       Weibull = {intervals(weibullm(inc , time, start))$coef[1,]},
       logistic5P = {intervals(logistic5pm (inc , time, start))$coef[1,]},
       all = {p<-rbind(intervals(Richardsm(inc , time, start[[1]]))$coef[1,],
                  intervals(logistic3pm(inc , time, start[[2]]))$coef[1,],
                  intervals(SigmEmaxm (inc , time, start[[3]]))$coef[1,],
                  intervals(gompertzm(inc , time, start[[4]]))$coef[1,],
                  intervals(weibullm(inc , time, start[[5]]))$coef[1,],
                  intervals(logistic5pm (inc , time, start[[6]]))$coef[1,],
                  estimFSizeMA(inc , time, start))
                  rownames(p)<-c("Richards","3P logistic","SigmEmax","Gompertz",
                                  "Weibull","5P logistic","Model averaged")
                  p},
       {cat("This model name is not correct\n")})
}


#### turning point Confidence interval  function
CIturningpoint<-function(inc , time, start=NULL, model){
switch(model,
       Richards = {intervals(Richardsm(inc , time, start))$coef[4,]},
       logistic3P = {intervals(logistic3pm(inc , time, start))$coef[3,]},
       SigmEmax = {intervals(SigmEmaxm (inc , time, start))$coef[2,]},
       Gompertz = {intervals(gompertzm(inc , time, start))$coef[2,]},
       Weibull = {intervals(weibullm(inc , time, start))$coef[2,]},
       logistic5P = {intervals(logistic5pm (inc , time, start))$coef[4,]},
       all = {p<-rbind(intervals(Richardsm(inc , time, start[[1]]))$coef[4,],
                  intervals(logistic3pm(inc , time, start[[2]]))$coef[3,],
                  intervals(SigmEmaxm (inc , time, start[[3]]))$coef[2,],
                  intervals(gompertzm(inc , time, start[[4]]))$coef[2,],
                  intervals(weibullm(inc , time, start[[5]]))$coef[2,],
                  intervals(logistic5pm (inc , time, start[[6]]))$coef[4,],
                  estimTpointMA(inc , time, start))
                  rownames(p)<-c("Richards","3P logistic","SigmEmax","Gompertz",
                                  "Weibull","5P logistic","Model averaged")
                  p},
       {cat("This model name is not correct\n")})
}


allmodels<-function(inc , time, start=NULL, model){
switch(model,
       Richards = {w<-Richardsm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-Richardsm(inc , time, start)$tTable
                   r2<-intervals(Richardsm(inc , time, start))$coef[1,]
                   r3<-intervals(Richardsm(inc , time, start))$coef[4,]
                   r4<-predict(Richardsm(inc , time, start))
                   pp1<-D(expression(alpha*((1+k*exp(-gamma*k*(t-eta)))^(-1/k))),"t")
                   alpha<-coef(Richardsm(inc , time, start))[1]
                   k<-coef(Richardsm(inc , time, start))[2]
                   gamma<-coef(Richardsm(inc , time, start))[3]
                   eta<-coef(Richardsm(inc , time, start))[4]
                   t<-time
                   r5<-eval(pp1)
			output = list(Incidence=inc, Time=time,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodels",model.type="Richards")
			class(output) = "dengue"	
			return(output)
                  },
      logistic3P = {w<-logistic3pm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-logistic3pm(inc , time, start)$tTable
                   r2<-intervals(logistic3pm(inc , time, start))$coef[1,]
                   r3<-intervals(logistic3pm(inc , time, start))$coef[3,]
                   r4<-predict(logistic3pm(inc , time, start))
                   pp2<-D(expression(alpha/(1+exp(-gamma*(t-eta)))),"t")
                   alpha<-coef(logistic3pm(inc , time, start))[1]
                   gamma<-coef(logistic3pm(inc , time, start))[2]
                   eta<-coef(logistic3pm(inc , time, start))[3]
                   t<-time
                   r5<-eval(pp2)
			output = list(Incidence=inc, Time=time,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodels",model.type="logistic3P")
			class(output) = "dengue"	
			return(output)
                  },
       SigmEmax = {w<-SigmEmaxm (inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-SigmEmaxm (inc , time, start)$tTable
                   r2<-intervals(SigmEmaxm (inc , time, start))$coef[1,]
                   r3<-intervals(SigmEmaxm (inc , time, start))$coef[2,]
                   r4<-predict(SigmEmaxm (inc , time, start))
                   pp3<-D(expression(e0 + (t^n)*(alpha-e0)/(t^n + eta^n)),"t")
                   alpha<-coef(SigmEmaxm (inc , time, start))[1]
                   eta<-coef(SigmEmaxm (inc , time, start))[2]
                   e0<-coef(SigmEmaxm (inc , time, start))[3]
                   n<-coef(SigmEmaxm (inc , time, start))[4]
                   t<-time
                   r5<-eval(pp3)
			output = list(Incidence=inc, Time=time,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodels",model.type="SigmEmax")
			class(output) = "dengue"	
			return(output)
                  },
       Gompertz = {w<-gompertzm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-gompertzm(inc , time, start)$tTable
                   r2<-intervals(gompertzm(inc , time, start))$coef[1,]
                   r3<-intervals(gompertzm(inc , time, start))$coef[2,]
                   r4<-predict(gompertzm(inc , time, start))
                   pp4<-D(expression(e0 + (alpha-e0)*exp(-exp(-k*(t-eta)))),"t")
                   alpha<-coef(gompertzm(inc , time, start))[1]
                   eta<-coef(gompertzm(inc , time, start))[2]
                   e0<-coef(gompertzm(inc , time, start))[3]
                   k<-coef(gompertzm(inc , time, start))[4]
                   t<-time
                   r5<-eval(pp4)
			output = list(Incidence=inc, Time=time,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodels",model.type="Gompertz")
			class(output) = "dengue"	
			return(output)
                  },
       Weibull =  {w<-weibullm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-weibullm(inc , time, start)$tTable
                   r2<-intervals(weibullm(inc , time, start))$coef[1,]
                   r3<-intervals(weibullm(inc , time, start))$coef[2,]
                   r4<-predict(weibullm(inc , time, start))
                   pp5<-D(expression(e0 + (alpha-e0)*exp(-(t/eta)^-gamma)),"t")
                   alpha<-coef(weibullm(inc , time, start))[1]
                   eta<-coef(weibullm(inc , time, start))[2]
                   e0<-coef(weibullm(inc , time, start))[3]
                   gamma<-coef(weibullm(inc , time, start))[4]
                   t<-time
                   r5<-eval(pp5)
			output = list(Incidence=inc, Time=time,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodels",model.type="Weibull")
			class(output) = "dengue"	
			return(output)
                  },
     logistic5P = {w<-logistic5pm (inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-logistic5pm (inc , time, start)$tTable
                   r2<-intervals(logistic5pm (inc , time, start))$coef[1,]
                   r3<-intervals(logistic5pm (inc , time, start))$coef[4,]
                   r4<-predict(logistic5pm (inc , time, start))
                   pp6<-D(expression(d+ (alpha-d)/(1+(2**(1/g)-1)*(t/eta)**-b)**g),"t")
                   alpha<-coef(logistic5pm (inc , time, start))[1]
                   d<-coef(logistic5pm (inc , time, start))[2]
                   g<-coef(logistic5pm (inc , time, start))[3]
                   eta<-coef(logistic5pm (inc , time, start))[4]
                   b<-coef(logistic5pm (inc , time, start))[4]
                   t<-time
                   r5<-eval(pp6)
			output = list(Incidence=inc, Time=time,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodels",model.type="logistic5P")
			class(output) = "dengue"	
			return(output)
                  },
       all = {w<-AIC(Richardsm(inc , time, start[[1]]),logistic3pm(inc , time, start[[2]]),SigmEmaxm (inc , time, start[[3]]),
                     gompertzm(inc , time, start[[4]]),weibullm(inc , time, start[[5]]),logistic5pm (inc , time, start[[6]]))[,2]
                 names(w)<-c("Richards","3P Logistic","Sigmoidal Emax","Gompertz","Weibull","5P Logistic")
               p0<-w
               p1<-weight(inc , time, start)
               p2<-CIFinalsize(inc , time, start,model)
               p3<-CIturningpoint(inc , time, start, model)
               p41<-predict(Richardsm(inc , time, start[[1]]))
               p42<-predict(logistic3pm(inc , time, start[[2]]))
               p43<-predict(SigmEmaxm (inc , time, start[[3]]))
               p44<-predict(gompertzm(inc , time, start[[4]]))
               p45<-predict(weibullm(inc , time, start[[5]]))
               p46<-predict(logistic5pm (inc , time, start[[6]]))
               p4<-list(p41,p42,p43,p44,p45,p46)
               p5<-weight(inc , time, start)[1]*p41 + weight(inc , time, start)[2]*p42 + weight(inc , time, start)[3]*p43 +weight(inc , time, start)[4]*p44 +weight(inc , time, start)[5]*p45 + weight(inc , time, start)[6]*p46
               pp1<-D(expression(alpha*((1+k*exp(-gamma*k*(t-eta)))^(-1/k))),"t")
               alpha<-coef(Richardsm(inc , time, start[1]))[1]
               k<-coef(Richardsm(inc , time, start[1]))[2]
               gamma<-coef(Richardsm(inc , time, start[1]))[3]
               eta<-coef(Richardsm(inc , time, start[1]))[4]
               t<-time
               p61<-eval(pp1)
               pp2<-D(expression(alpha/(1+exp(-gamma*(t-eta)))),"t")
               alpha<-coef(logistic3pm(inc , time, start[2]))[1]
               gamma<-coef(logistic3pm(inc , time, start[2]))[2]
               eta<-coef(logistic3pm(inc , time, start[2]))[3]
               p62<-eval(pp2)
               pp3<-D(expression(e0 + (t^n)*(alpha-e0)/(t^n + eta^n)),"t")
               alpha<-coef(SigmEmaxm (inc , time, start[3]))[1]
               eta<-coef(SigmEmaxm (inc , time, start[3]))[2]
               e0<-coef(SigmEmaxm (inc , time, start[3]))[3]
                n<-coef(SigmEmaxm (inc , time, start[3]))[4]
               p63<-eval(pp3)
               pp4<-D(expression(e0 + (alpha-e0)*exp(-exp(-k*(t-eta)))),"t")
               alpha<-coef(gompertzm(inc , time, start[4]))[1]
               eta<-coef(gompertzm(inc , time, start[4]))[2]
               e0<-coef(gompertzm(inc , time, start[4]))[3]
               k<-coef(gompertzm(inc , time, start[4]))[4]
               p64<-eval(pp4)
               pp5<-D(expression(e0 + (alpha-e0)*exp(-(t/eta)^-gamma)),"t")
               alpha<-coef(weibullm(inc , time, start[5]))[1]
               eta<-coef(weibullm(inc , time, start[5]))[2]
               e0<-coef(weibullm(inc , time, start[5]))[3]
               gamma<-coef(weibullm(inc , time, start[5]))[4]
               p65<-eval(pp5)
               pp6<-D(expression(d+ (alpha-d)/(1+(2**(1/g)-1)*(t/eta)**-b)**g),"t")
               alpha<-coef(logistic5pm (inc , time, start[6]))[1]
               d<-coef(logistic5pm (inc , time, start[6]))[2]
               g<-coef(logistic5pm (inc , time, start[6]))[3]
               eta<-coef(logistic5pm (inc , time, start[6]))[4]
               b<-coef(logistic5pm (inc , time, start[6]))[5]
               p66<-eval(pp6)
               p6<-list(p61,p62,p63,p64,p65,p66)
               p7<-weight(inc , time, start)[1]*p61 + weight(inc , time, start)[2]*p62 + weight(inc , time, start)[3]*p63 +weight(inc , time, start)[4]*p64 +weight(inc , time, start)[5]*p65 + weight(inc , time, start)[6]*p66
               output = list(Incidence=inc, Time=time,AIC=p0,Weights=p1,FinalSize=p2,TurningPoint=p3,Predict=p4,PredictMA=p5,PredInc=p6, PredMAinc=p7,function.type="allmodels",model.type="all")
               class(output) = "dengue"
		   return(output)},
     {stop("This model name is not correct")})
}


### SUMMARY FUNCTION
summary.dengue <- function(object,...){
   if(object$function.type=="allmodels"){
           if(object$model.type!="all"){
                if(object$model.type=="Richards")
                {cat("\nRichards model\n")}
                if(object$model.type=="logistic3P")
                {cat("\n3P logistic model\n")}
                if(object$model.type=="SigmEmax")
                {cat("\nSigmoid Emax model\n")}
                if(object$model.type=="Gompertz")
                {cat("\nGompertz model\n")}
                if(object$model.type=="Weibull")
                {cat("\nWeibull model\n")}
                if(object$model.type=="logistic5P")
                {cat("\n5P logistic model\n")}
                cat("\nAIC")
                print(object$AIC)
                cat("\nParameter estimate\n")
                print(object$tTable)
                cat("\nFinal size estimate:\n")
                print(object$FinalSize)
                cat("\nTurning point estimate:\n")
                print(object$TurningPoint)}
          else {cat("\nAIC\n")
                print(object$AIC)
                cat("\nModel weights\n")
                print(object$Weights)
                cat("\nModel averaged Final size estimate\n")
                print(object$FinalSize)
                cat("\nModel averaged Turning point estimate\n")
                print(object$TurningPoint)
                }
          }
       if(object$function.type=="allmodelpredict"){
            if(object$model.type!="all"){
                if(object$model.type=="Richards")
                {cat("\nRichards model\n")}
                if(object$model.type=="logistic3P")
                {cat("\n3P logistic model\n")}
                if(object$model.type=="SigmEmax")
                {cat("\nSigmoid Emax model\n")}
                if(object$model.type=="Gompertz")
                {cat("\nGompertz model\n")}
                if(object$model.type=="Weibull")
                {cat("\nWeibull model\n")}
                if(object$model.type=="logistic5P")
                {cat("\n5P logistic model\n")}
                cat("\nAIC")
                print(object$AIC)
                cat("\nPrediction of the Final size at the end of epidemic:\n")
                print(object$FinalSize)
                cat("\nPrediction of Turning point at the end of epidemic:\n")
                print(object$TurningPoint)
                cat("\nPrediction of the incidence at the time point ",length(object$PredTime)," \n")
                print(object$PredInc[[length(object$PredTime)]])
                cat("\nPrediction of the cummulative number of cases at the time point ",length(object$PredTime)," \n")
                print(object$Predict[[length(object$PredTime)]])}
          else {cat("\nAIC\n")
                print(object$AIC)
                cat("\nModel weights\n")
                print(object$Weights)
                cat("\nModel averaged prediction of the Final size at the end of epidemic\n")
                print(object$FinalSize)
                cat("\nModel averaged prediction of the Turning point at the end of epidemic\n")
                print(object$TurningPoint)
                cat("\nModel averaged prediction of the incidence at the time point ",length(object$PredTime)," \n")
                print(object$PredMAinc[[length(object$PredTime)]])
                cat("\nModel averaged prediction of the cummulative number of cases at the time point ",length(object$PredTime)," \n")
                print(object$PredictMA[[length(object$PredTime)]])}
          }
       if(object$function.type=="changetimeFSTP"){
             cat("\nChanges over time of the final size\n")
              print(object$FSchangetime)
              cat("\nChanges over time of the turning point\n")
              print(object$TPchangetime)
           }
}


#### PLOT FUNCTION
plot.dengue= function(x,which=c(1,2),xlab="",...){
             L<-NULL
             U<-NULL
             t<-as.character(c("Richards","3P logistic","Sigmoidal Emax","Gompertz","Weibull","5p logistic","Model averaged"))
             t<-factor(t,levels=c("Richards","3P logistic","Sigmoidal Emax","Gompertz","Weibull","5p logistic","Model averaged"))
     if( 1 %in% which){
         if(x$function.type=="allmodels"){
             dev.new()
		if(x$model.type == "Richards"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict,col=2,lwd=2,lty=4)}
            if(x$model.type == "logistic3P"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict,col=3,lwd=2,lty=4)}
            if(x$model.type == "SigmEmax"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict,col=4,lwd=2,lty=4)}
            if(x$model.type == "Gompertz"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict,col=5,lwd=2,lty=4)}
            if(x$model.type == "Weibull"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict,col=6,lwd=2,lty=4)}
            if(x$model.type == "logistic5P"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict,col=7,lwd=2,lty=4)}
            if(x$model.type == "all"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(cumsum(x$Incidence))+10),...)
              lines(x$Time,x$Predict[[1]],col=2,lwd=2,lty=4)
              lines(x$Time,x$Predict[[2]],col=3,lwd=2,lty=4)
              lines(x$Time,x$Predict[[3]],col=4,lwd=2,lty=4)
              lines(x$Time,x$Predict[[4]],col=5,lwd=2,lty=4)
              lines(x$Time,x$Predict[[5]],col=6,lwd=2,lty=4)
              lines(x$Time,x$Predict[[6]],col=7,lwd=2,lty=4)
              lines(x$Time,x$PredictMA,col=8,lwd=3,lty=1)}
          }
          else {dev.new()
		if(x$model.type == "Richards"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$Predict[[1]],cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict,col=2,lwd=2,lty=4)}
            if(x$model.type == "logistic3P"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$Predict[[2]],cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict,col=3,lwd=2,lty=4)}
            if(x$model.type == "SigmEmax"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$Predict[[3]],cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict,col=4,lwd=2,lty=4)}
            if(x$model.type == "Gompertz"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$Predict[[4]],cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict,col=5,lwd=2,lty=4)}
            if(x$model.type == "Weibull"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$Predict[[5]],cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict,col=6,lwd=2,lty=4)}
            if(x$model.type == "logistic5P"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$Predict[[6]],cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict,col=7,lwd=2,lty=4)}
            if(x$model.type == "all"){
              plot(x$Time,cumsum(x$Incidence),pch=19,xlab=xlab,ylab="Cumulative number of cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredictMA,cumsum(x$Incidence))+100),...)
              lines(x$PredTime,x$Predict[[1]],col=2,lwd=2,lty=4)
              lines(x$PredTime,x$Predict[[2]],col=3,lwd=2,lty=4)
              lines(x$PredTime,x$Predict[[3]],col=4,lwd=2,lty=4)
              lines(x$PredTime,x$Predict[[4]],col=5,lwd=2,lty=4)
              lines(x$PredTime,x$Predict[[5]],col=6,lwd=2,lty=4)
              lines(x$PredTime,x$Predict[[6]],col=7,lwd=2,lty=4)
              lines(x$PredTime,x$PredictMA,col=8,lwd=3,lty=1)}
          }      
     }
       if( 2 %in% which){
           if(x$function.type=="allmodels"){
             dev.new()
		if(x$model.type == "Richards"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc,col=2,lwd=2,lty=4)}
            if(x$model.type == "logistic3P"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc,col=3,lwd=2,lty=4)}
            if(x$model.type == "SigmEmax"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc,col=4,lwd=2,lty=4)}
            if(x$model.type == "Gompertz"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc,col=5,lwd=2,lty=4)}
            if(x$model.type == "Weibull"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc,col=6,lwd=2,lty=4)}
            if(x$model.type == "logistic5P"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc,col=7,lwd=2,lty=4)}
            if(x$model.type == "all"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$Incidence)+3),ylim=c(0,max(x$Incidence)+5),...)
              lines(x$Time,x$PredInc[[1]],col=2,lwd=2,lty=4)
              lines(x$Time,x$PredInc[[2]],col=3,lwd=2,lty=4)
              lines(x$Time,x$PredInc[[3]],col=4,lwd=2,lty=4)
              lines(x$Time,x$PredInc[[4]],col=5,lwd=2,lty=4)
              lines(x$Time,x$PredInc[[5]],col=6,lwd=2,lty=4)
              lines(x$Time,x$PredInc[[6]],col=7,lwd=2,lty=4)
              lines(x$Time,x$PredMAinc,col=8,lwd=3,lty=1)}
          }
          else {dev.new()
		if(x$model.type == "Richards"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredInc[[1]],x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc,col=2,lwd=2,lty=4)}
            if(x$model.type == "logistic3P"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredInc[[2]],x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc,col=3,lwd=2,lty=4)}
            if(x$model.type == "SigmEmax"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredInc[[3]],x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc,col=4,lwd=2,lty=4)}
            if(x$model.type == "Gompertz"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredInc[[4]],x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc,col=5,lwd=2,lty=4)}
            if(x$model.type == "Weibull"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredInc[[5]],x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc,col=6,lwd=2,lty=4)}
            if(x$model.type == "logistic5P"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredInc[[6]],x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc,col=7,lwd=2,lty=4)}
            if(x$model.type == "all"){
              plot(x$Time,x$Incidence,type="h",lwd=6,xlab=xlab,ylab="Reported cases",xlim=c(0,length(x$PredTime)+3),ylim=c(0,max(x$PredMAinc,x$Incidence)+5),...)
              lines(x$PredTime,x$PredInc[[1]],col=2,lwd=2,lty=4)
              lines(x$PredTime,x$PredInc[[2]],col=3,lwd=2,lty=4)
              lines(x$PredTime,x$PredInc[[3]],col=4,lwd=2,lty=4)
              lines(x$PredTime,x$PredInc[[4]],col=5,lwd=2,lty=4)
              lines(x$PredTime,x$PredInc[[5]],col=6,lwd=2,lty=4)
              lines(x$PredTime,x$PredInc[[6]],col=7,lwd=2,lty=4)
              lines(x$PredTime,x$PredMAinc,col=8,lwd=3,lty=1)}
         }
       } 
     if( 3 %in% which){
              if(x$model.type!="all") {
                stop("Cannot produce this plot, because you should use all the models",call.=FALSE)
		  } 
             else {
             dev.new()
		 df <- data.frame(t = t,
                 F = x$FinalSize[,2],
                 L = x$FinalSize[,1],
                 U = x$FinalSize[,3])
              print(ggplot(df, aes(x = t, y = F)) +
               geom_point(size = 4) +  
               geom_errorbar(aes(ymax = U, ymin = L))+ 
               coord_flip(xlim = c(0,8), ylim = c(min(x$FinalSize[,1])-3,max(x$FinalSize[,3])+3)) +
               xlab("Models") + ylab("Final Size"))           
             }
         }
     if( 4 %in% which){
               if(x$model.type!="all") {
                 stop("Cannot produce this plot, because you should use all the models",call.=FALSE)
		     } 
               else {
		   dev.new()
               
              df<- data.frame(t =t,
                 F =x$TurningPoint[,2],
                 L =x$TurningPoint[,1],
                 U =x$TurningPoint[,3])
                print(ggplot(df, aes(x = t, y = F)) +
                 geom_point(size = 4) +  
                 geom_errorbar(aes(ymax = U, ymin = L))+ 
                  coord_flip(xlim = c(0,8), ylim = c(min(x$TurningPoint[,1])-0.1,max(x$TurningPoint[,3])+0.1)) +
                  xlab("Models") + ylab("Turning Point"))
                }
               } 
    if( 5 %in% which){
           if(x$function.type!="changetimeFSTP") {
                 stop("Cannot produce this plot, because you should use the function changetimeFSTP",call.=FALSE)
		     } 
               else {
		   dev.new()
            plot(x$Period,x$FSchangetime[,1],type="l",col=2,lwd=2,lty=4,xlab="Period",ylab="Final size of outbreak",xlim=c(min(x$Period)-1,max(x$Period)+1),ylim=c(min(x$FSchangetime)-10,max(x$FSchangetime)+10),...)
              lines(x$Period,x$FSchangetime[,2],col=3,lwd=2,lty=4)
              lines(x$Period,x$FSchangetime[,3],col=4,lwd=2,lty=4)
              lines(x$Period,x$FSchangetime[,4],col=5,lwd=2,lty=4)
              lines(x$Period,x$FSchangetime[,5],col=6,lwd=2,lty=4)
              lines(x$Period,x$FSchangetime[,6],col=7,lwd=2,lty=4)
              lines(x$Period,x$FSchangetime[,7],col=8,lwd=3,lty=1)}  
         } 
    if( 6 %in% which){
           if(x$function.type!="changetimeFSTP") {
                 stop("Cannot produce this plot, because you should use the function changetimeFSTP",call.=FALSE)
		     } 
               else {
		   dev.new()
            plot(x$Period,x$TPchangetime[,1],type="l",col=2,lwd=2,lty=4,xlab="Period",ylab="Turning point of outbreak",xlim=c(min(x$Period)-1,max(x$Period)+1),ylim=c(min(x$TPchangetime)-3,max(x$TPchangetime)+3),...)
              lines(x$Period,x$TPchangetime[,2],col=3,lwd=2,lty=4)
              lines(x$Period,x$TPchangetime[,3],col=4,lwd=2,lty=4)
              lines(x$Period,x$TPchangetime[,4],col=5,lwd=2,lty=4)
              lines(x$Period,x$TPchangetime[,5],col=6,lwd=2,lty=4)
              lines(x$Period,x$TPchangetime[,6],col=7,lwd=2,lty=4)
              lines(x$Period,x$TPchangetime[,7],col=8,lwd=3,lty=1)}  
         } 
                
}

 
     
allmodelpredict<-function(inc , time, pred, start=NULL, model){
switch(model,
       Richards = {w<-Richardsm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-Richardsm(inc , time, start)$tTable
                   r2<-intervals(Richardsm(inc , time, start))$coef[1,]
                   r3<-intervals(Richardsm(inc , time, start))$coef[4,]
                   pp1<-D(expression(alpha*((1+k*exp(-gamma*k*(t-eta)))^(-1/k))),"t")
                   alpha<-coef(Richardsm(inc , time, start))[1]
                   k<-coef(Richardsm(inc , time, start))[2]
                   gamma<-coef(Richardsm(inc , time, start))[3]
                   eta<-coef(Richardsm(inc , time, start))[4]
                   t<-c(1:pred)
                   r5<-eval(pp1)
                   r4<-alpha*(1+k*exp(-gamma*k*(t-eta)))^(-1/k)
			output = list(Incidence=inc, Time=time,PredTime=t, AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodelpredict",model.type="Richards")
			class(output) = "dengue"	
			return(output)
                  },
      logistic3P = {w<-logistic3pm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-logistic3pm(inc , time, start)$tTable
                   r2<-intervals(logistic3pm(inc , time, start))$coef[1,]
                   r3<-intervals(logistic3pm(inc , time, start))$coef[3,]
                   pp2<-D(expression(alpha/(1+exp(-gamma*(t-eta)))),"t")
                   alpha<-coef(logistic3pm(inc , time, start))[1]
                   gamma<-coef(logistic3pm(inc , time, start))[2]
                   eta<-coef(logistic3pm(inc , time, start))[3]
                   t<-c(1:pred)
                   r4<-alpha/(1+exp(-gamma*(t-eta)))
                   r5<-eval(pp2)
			output = list(Incidence=inc, Time=time,PredTime=t,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodelpredict",model.type="logistic3P")
			class(output) = "dengue"	
			return(output)
                  },
       SigmEmax = {w<-SigmEmaxm (inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-SigmEmaxm (inc , time, start)$tTable
                   r2<-intervals(SigmEmaxm (inc , time, start))$coef[1,]
                   r3<-intervals(SigmEmaxm (inc , time, start))$coef[2,]
                   pp3<-D(expression(e0 + (t^n)*(alpha-e0)/(t^n + eta^n)),"t")
                   alpha<-coef(SigmEmaxm (inc , time, start))[1]
                   eta<-coef(SigmEmaxm (inc , time, start))[2]
                   e0<-coef(SigmEmaxm (inc , time, start))[3]
                   n<-coef(SigmEmaxm (inc , time, start))[4]
                   t<-c(1:pred)
                   r5<-eval(pp3)
                   r4<-e0 + (t^n)*(alpha-e0)/(t^n + eta^n)
			output = list(Incidence=inc, Time=time,PredTime=t,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodelpredict",model.type="SigmEmax")
			class(output) = "dengue"	
			return(output)
                  },
       Gompertz = {w<-gompertzm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-gompertzm(inc , time, start)$tTable
                   r2<-intervals(gompertzm(inc , time, start))$coef[1,]
                   r3<-intervals(gompertzm(inc , time, start))$coef[2,]
                   pp4<-D(expression(e0 + (alpha-e0)*exp(-exp(-k*(t-eta)))),"t")
                   alpha<-coef(gompertzm(inc , time, start))[1]
                   eta<-coef(gompertzm(inc , time, start))[2]
                   e0<-coef(gompertzm(inc , time, start))[3]
                   k<-coef(gompertzm(inc , time, start))[4]
                   t<-c(1:pred)
                   r5<-eval(pp4)
                   r4<-e0 + (alpha-e0)*exp(-exp(-k*(t-eta)))
			output = list(Incidence=inc, Time=time,PredTime=t,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodelpredict",model.type="Gompertz")
			class(output) = "dengue"	
			return(output)
                  },
       Weibull =  {w<-weibullm(inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-weibullm(inc , time, start)$tTable
                   r2<-intervals(weibullm(inc , time, start))$coef[1,]
                   r3<-intervals(weibullm(inc , time, start))$coef[2,]
                   pp5<-D(expression(e0 + (alpha-e0)*exp(-(t/eta)^-gamma)),"t")
                   alpha<-coef(weibullm(inc , time, start))[1]
                   eta<-coef(weibullm(inc , time, start))[2]
                   e0<-coef(weibullm(inc , time, start))[3]
                   gamma<-coef(weibullm(inc , time, start))[4]
                   t<-c(1:pred)
                   r5<-eval(pp5)
                   r4<-e0 + (alpha-e0)*exp(-(t/eta)^-gamma)
			output = list(Incidence=inc, Time=time,PredTime=t,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodelpredict",model.type="Weibull")
			class(output) = "dengue"	
			return(output)
                  },
     logistic5P = {w<-logistic5pm (inc , time, start)$AIC
                   names(w)<-c("")
                   r0<-w
                   r1<-logistic5pm (inc , time, start)$tTable
                   r2<-intervals(logistic5pm (inc , time, start))$coef[1,]
                   r3<-intervals(logistic5pm (inc , time, start))$coef[4,]
                   pp6<-D(expression(d+ (alpha-d)/(1+(2**(1/g)-1)*(t/eta)**-b)**g),"t")
                   alpha<-coef(logistic5pm (inc , time, start))[1]
                   d<-coef(logistic5pm (inc , time, start))[2]
                   g<-coef(logistic5pm (inc , time, start))[3]
                   eta<-coef(logistic5pm (inc , time, start))[4]
                   b<-coef(logistic5pm (inc , time, start))[4]
                   t<-c(1:pred)
                   r5<-eval(pp6)
                   r4<-d+ (alpha-d)/(1+(2**(1/g)-1)*(t/eta)**-b)**g
			output = list(Incidence=inc, Time=time,PredTime=t,AIC=r0,tTable=r1,FinalSize=r2,TurningPoint=r3,Predict=r4,PredInc=r5,function.type="allmodelpredict",model.type="logistic5P")
			class(output) = "dengue"	
			return(output)
                  },
       all = {w<-AIC(Richardsm(inc , time, start[[1]]),logistic3pm(inc , time, start[[2]]),SigmEmaxm (inc , time, start[[3]]),
                     gompertzm(inc , time, start[[4]]),weibullm(inc , time, start[[5]]),logistic5pm (inc , time, start[[6]]))[,2]
                 names(w)<-c("Richards","3P Logistic","Sigmoidal Emax","Gompertz","Weibull","5P Logistic")
               p0<-w
               p1<-weight(inc , time, start)
               p2<-CIFinalsize(inc , time, start,model)
               p3<-CIturningpoint(inc , time, start, model)
               pp1<-D(expression(alpha*((1+k*exp(-gamma*k*(t-eta)))^(-1/k))),"t")
               alpha<-coef(Richardsm(inc , time, start[1]))[1]
               k<-coef(Richardsm(inc , time, start[1]))[2]
               gamma<-coef(Richardsm(inc , time, start[1]))[3]
               eta<-coef(Richardsm(inc , time, start[1]))[4]
               t<-c(1:pred)
               p61<-eval(pp1)
               p41<-alpha*((1+k*exp(-gamma*k*(t-eta)))^(-1/k))
               pp2<-D(expression(alpha/(1+exp(-gamma*(t-eta)))),"t")
               alpha<-coef(logistic3pm(inc , time, start[2]))[1]
               gamma<-coef(logistic3pm(inc , time, start[2]))[2]
               eta<-coef(logistic3pm(inc , time, start[2]))[3]
               p62<-eval(pp2)
               p42<-alpha/(1+exp(-gamma*(t-eta)))
               pp3<-D(expression(e0 + (t^n)*(alpha-e0)/(t^n + eta^n)),"t")
               alpha<-coef(SigmEmaxm (inc , time, start[3]))[1]
               eta<-coef(SigmEmaxm (inc , time, start[3]))[2]
               e0<-coef(SigmEmaxm (inc , time, start[3]))[3]
                n<-coef(SigmEmaxm (inc , time, start[3]))[4]
               p63<-eval(pp3)
               p43<-e0 + (t^n)*(alpha-e0)/(t^n + eta^n)
               pp4<-D(expression(e0 + (alpha-e0)*exp(-exp(-k*(t-eta)))),"t")
               alpha<-coef(gompertzm(inc , time, start[4]))[1]
               eta<-coef(gompertzm(inc , time, start[4]))[2]
               e0<-coef(gompertzm(inc , time, start[4]))[3]
               k<-coef(gompertzm(inc , time, start[4]))[4]
               p64<-eval(pp4)
               p44<-e0 + (alpha-e0)*exp(-exp(-k*(t-eta)))
               pp5<-D(expression(e0 + (alpha-e0)*exp(-(t/eta)^-gamma)),"t")
               alpha<-coef(weibullm(inc , time, start[5]))[1]
               eta<-coef(weibullm(inc , time, start[5]))[2]
               e0<-coef(weibullm(inc , time, start[5]))[3]
               gamma<-coef(weibullm(inc , time, start[5]))[4]
               p65<-eval(pp5)
               p45<-e0 + (alpha-e0)*exp(-(t/eta)^-gamma)
               pp6<-D(expression(d+ (alpha-d)/(1+(2**(1/g)-1)*(t/eta)**-b)**g),"t")
               alpha<-coef(logistic5pm (inc , time, start[6]))[1]
               d<-coef(logistic5pm (inc , time, start[6]))[2]
               g<-coef(logistic5pm (inc , time, start[6]))[3]
               eta<-coef(logistic5pm (inc , time, start[6]))[4]
               b<-coef(logistic5pm (inc , time, start[6]))[5]
               p66<-eval(pp6)
               p46<-d+ (alpha-d)/(1+(2**(1/g)-1)*(t/eta)**-b)**g
               p6<-list(p61,p62,p63,p64,p65,p66)
               p7<-weight(inc , time, start)[1]*p61 + weight(inc , time, start)[2]*p62 + weight(inc , time, start)[3]*p63 +weight(inc , time, start)[4]*p64 +weight(inc , time, start)[5]*p65 + weight(inc , time, start)[6]*p66
               p4<-list(p41,p42,p43,p44,p45,p46)
               p5<-weight(inc , time, start)[1]*p41 + weight(inc , time, start)[2]*p42 + weight(inc , time, start)[3]*p43 +weight(inc , time, start)[4]*p44 +weight(inc , time, start)[5]*p45 + weight(inc , time, start)[6]*p46
               output = list(Incidence=inc, Time=time,PredTime=t,AIC=p0,Weights=p1,FinalSize=p2,TurningPoint=p3,Predict=p4,PredictMA=p5,PredInc=p6, PredMAinc=p7,function.type="allmodelpredict",model.type="all")
               class(output) = "dengue"
		   return(output)},
     {stop("This model name is not correct")})
}


changetimeFSTP<-function(inc, time,ini,start=NULL){
           y1<-NULL  
           y2<-NULL
           y3<-NULL
           y4<-NULL
           y5<-NULL
           y6<-NULL
           y7<-NULL 
           t1<-NULL  
           t2<-NULL
           t3<-NULL
           t4<-NULL
           t5<-NULL
           t6<-NULL
           t7<-NULL  
           p<-NULL
           for (i in ini:length(time))
           {p[[i]]<-allmodels(inc[1:i],time[1:i], start=start, model="all")
           y1[[i-ini+1]]<-p[[i]]$FinalSize[1,2]
           y2[[i-ini+1]]<-p[[i]]$FinalSize[2,2]
           y3[[i-ini+1]]<-p[[i]]$FinalSize[3,2]
           y4[[i-ini+1]]<-p[[i]]$FinalSize[4,2]
           y5[[i-ini+1]]<-p[[i]]$FinalSize[5,2]
           y6[[i-ini+1]]<-p[[i]]$FinalSize[6,2]
           y7[[i-ini+1]]<-p[[i]]$FinalSize[7,2]
           t1[[i-ini+1]]<-p[[i]]$TurningPoint[1,2]
           t2[[i-ini+1]]<-p[[i]]$TurningPoint[2,2]
           t3[[i-ini+1]]<-p[[i]]$TurningPoint[3,2]
           t4[[i-ini+1]]<-p[[i]]$TurningPoint[4,2]
           t5[[i-ini+1]]<-p[[i]]$TurningPoint[5,2]
           t6[[i-ini+1]]<-p[[i]]$TurningPoint[6,2]
           t7[[i-ini+1]]<-p[[i]]$TurningPoint[7,2]}
           y<-c(y1,y2,y3,y4,y5,y6,y7)
           t<-c(t1,t2,t3,t4,t5,t6,t7)
           d<-matrix(y,nrow=length(time)-ini+1, ncol=7, byrow=FALSE, dimnames=list(paste("1",1:length(time),sep="-")[ini:length(time)],c("Richards","logistic3P","SigmEmax","Gompertz","Weibull","logistic5P","Model averaged")))
           tt<-matrix(t,nrow=length(time)-ini+1, ncol=7, byrow=FALSE, dimnames=list(paste("1",1:length(time),sep="-")[ini:length(time)],c("Richards","logistic3P","SigmEmax","Gompertz","Weibull","logistic5P","Model averaged")))
           g<-c(ini:length(time))
           output = list(Incidence=inc, Time=time,Period=g,FSchangetime=d,TPchangetime=tt,function.type="changetimeFSTP",model.type="all")
           class(output) = "dengue"
	     return(output)
}

Try the DengueParmEst package in your browser

Any scripts or data that you put into this service are public.

DengueParmEst documentation built on May 31, 2017, 3:03 a.m.