[R-es] sir

Jose Betancourt B. bet@n@ter @end|ng |rom gm@||@com
Jue Jun 11 22:09:28 CEST 2020


Estimados


La pregunta es como agregar una line de cas os reales a la de
prevalencia simulada?

saludos

Dr. Betancourt

2020-06-11 14:36 GMT-04:00, Jose Betancourt B. <betanster using gmail.com>:
> script completo y adjunto función
>
>
> rm(list = ls(all = TRUE))
> require("sfsmisc")
> source("sir_func.R")
> npop = 4400
> I_0 = 1
> R_0 = 0
> S_0 = npop-I_0-R_0
>
>
> tbegin = 0
> tend   = 300
> vt = seq(tbegin,tend,1)
> gamma = 1/5
> R0    = 1.20
> beta  = R0*gamma
>
> vparameters = c(gamma=gamma,beta=beta)
> inits = c(S=S_0,I=I_0,R=R_0)
>
> solved_model = as.data.frame(lsoda(inits, vt, derivative_calc_func,
> vparameters))
> cat("The item names in the solved_model object
> are:",names(solved_model),"\n")
>
> vS = solved_model$S
> vI = solved_model$I
> vR = solved_model$R
> vtime = solved_model$time
> vnpop = vS+vI+vR
>
> mult.fig(4,main="SIR model of pandemic coronavirus in Camaguey Bv with
> R0=1.2")
>
>
> ymax = 1.4*max(vI/npop)
> plot(vtime,vI/npop,type="l",xlab="time",ylab="fraction
> infected",ylim=c(0,ymax),lwd=3,col=4,main="Infected")
>
> n=length(vtime)
> lines(vtime[2:n],-diff(vS)/(diff(vtime)*vnpop[1:(n-1)]),type="l",lwd=3,col=2)
>
> legend("topright",legend=c("total infected (prevalence)","newly
> infected/day (incidence)"),bty="n",lwd=3,col=c(4,2))
>
>
> ymin = 0.9*min(vS/vnpop)
> plot(vtime,vS/vnpop,type="l",xlab="time",ylab="fraction
> susceptible",ylim=c(ymin,1),lwd=3,main="Susceptible")
>
> iind = which.min(abs(vS/vnpop-1/R0)) # find the index at which S/N is
> equal to 1/R0
> lines(c(vtime[iind],vtime[iind]),c(-1000,1000),col=3,lwd=3)
> legend("bottomleft",legend=c("time at which
> S=1/R0"),bty="n",lwd=3,col=c(3),cex=0.7)
>
> plot(vtime,log(vI/vnpop),type="l",xlab="time",ylab="log(fraction
> infected)",lwd=3,col=4,main="log(Infected)")
>
> text(40,-14,"Initial\n exponential\n rise",cex=0.7)
>
> lines(c(vtime[iind],vtime[iind]),c(-1000,1000),col=3,lwd=3)
> legend("topleft",legend=c("time at which
> S=1/R0","log(Infected)"),bty="n",lwd=3,col=c(3,4),cex=0.7)
>
> epsilon = 0.00001
> vsinf = seq(0,1-epsilon,epsilon)
> # LHS = -log(vsinf)+log(S_0/npop)
> # RHS = R0*(1-vsinf)
> LHS = -log(vsinf)+log(S_0/npop)
> RHS = R0*(1-vsinf)
> iind = which.min(abs(RHS-LHS))
> sinf_predicted = vsinf[iind]
>
> cat("The final fraction of susceptibles at the end of the epidemic
> from the model simulation is ",min(vS/vnpop),"\n")
> cat("The final fraction of susceptibles at the end of the epidemic
> predicted by the final size relation is ",sinf_predicted,"\n")
>
>
>
> --
> Dr. Jose A. Betancourt Bethencourt
> Universidad de Ciencias Medicas Carlos j. Finlay
>


-- 
Dr. Jose A. Betancourt Bethencourt
Universidad de Ciencias Medicas Carlos j. Finlay



Más información sobre la lista de distribución R-help-es