install.packages("R0")

library(R0)

data(Germany.1918)

Germany.1918

mGT<-generation.time("gamma", c(3, 1.5)) #see also est.GT

#plot(mGT)

#data.frame(Germany.1918)

#plot(Germany.1918, col="red")

estR0<-estimate.R(Germany.1918, mGT, begin=1, end=27, methods=c("EG", "ML", "TD", "AR", "SB"),

pop.size=100000, nsim=1000)

estR0

#str(estR0)

par(mfrow=c(3,2), mar=c(2,2,2,2))

plot(estR0)

#individual plots

par(mfrow=c(2,3), mar=c(2,2,2,2))

plot(estR0$estimates$EG)

plot(estR0$estimates$ML)

plot(estR0$estimates$AR)

plot(estR0$estimates$TD)

plot(estR0$estimates$SB)

#par(mfrow=c(4,4), mar=c(2,2,2,2))

plotfit(estR0$estimates$SB) #works fine in R

##sensitivity analysis of the time window for exponential growth

sen=sensitivity.analysis(sa.type="time", incid=Germany.1918, GT=mGT, begin=1:15, end=16:30,

est.method="EG")

plot(sen, what=c("criterion","heatmap"))

##sensitivity analysis of the generation time

tmp<-sensitivity.analysis(sa.type="GT", incid=Germany.1918, GT.type="gamma", GT.mean=seq(1,5,1),

GT.sd.range=1, begin=1, end=27, est.method="EG")

par(mfrow=c(1,1), mar=c(4,4,4,4))

plot(x=tmp[,"GT.Mean"], xlab="mean GT (days)", y=tmp[,"R"], ylim=c(1.2, 2.1), ylab="R0 (95% CI)",

type="p", pch=19, col="black", main="Sensitivity of R0 to mean GT")

arrows(x0=as.numeric(tmp[,"GT.Mean"]), y0=as.numeric(tmp[,"CI.lower"]),

y1=as.numeric(tmp[,"CI.upper"]), angle=90, code=3, col="black", length=0.05)

^I am getting stuck on line 6, the first plot, but none of the other plots show either.