# R/LLE_functions.R In tsiR: An Implementation of the TSIR Model

#### Documented in plotLLEpredicttsirTSIR_LETSIR_LLE

```#' @title predicttsir
#' @description function to predict incidence and susceptibles using the tsir model. This is different than
#' simulatetsir as you are inputting parameters as vectors. The output is a data frame I and S with mean and confidence intervals of predictions.
#' @param times The time vector to predict the model from. This assumes that the time step is equal to IP
#' @param births The birth vector (of length length(times) or a single element) where each element is the births in that given (52/IP) time step
#' @param beta The length(52/IP) beta vector of contact.
#' @param alpha A single numeric which acts as the homogeniety parameter.
#' @param S0 The starting initial condition for S. This should be greater than one, i.e. not a fraction.
#' @param I0 The starting initial condition for I. This should be greater than one, i.e. not a fraction.
#' @param nsim The number of simulations to perform.
#' @param stochastic A TRUE / FALSE argument where FALSE is the deterministic model, and TRUE is a negative binomial distribution.
#' @export
#' @examples
#' \dontrun{
#' require(kernlab)
#' require(ggplot2)
#' require(kernlab)
#' require(tsiR)
#' London <- twentymeas\$London
#'
#' London <- subset(London, time > 1950)
#'
#' IP <- 2
#' ## first estimate paramters from the London data
#' parms <- estpars(data=London, IP=2, regtype='gaussian')
#'
## look at beta and alpha estimate
#' plotbeta(parms)
#'
#' ## now lets predict forward 20 years using the mean birth rate,
#' ## starting from rough initial conditions
#' births <- min(London\$births)
#' times <- seq(1965,1985, by = 1/ (52/IP))
#' S0 <- parms\$sbar
#' I0 <- 1e-5*mean(London\$pop)
#'
#' pred <- predicttsir(times=times,births=births,
#'                     beta=parms\$contact\$beta,alpha=parms\$alpha,
#'                     S0=S0,I0=I0,
#'                     nsim=50,stochastic=T)
#'
#' ## plot this prediction
#' ggplot(pred\$I,aes(time,mean))+geom_line()+geom_ribbon(aes(ymin=low,ymax=high),alpha=0.3)
#'
#'
#'}

predicttsir <- function(times,births,beta,alpha,S0,I0,nsim,stochastic){

I.mat <- S.mat <-  matrix(NA,length(times),nsim)

alpha <- alpha[1]

if(length(beta) < length(times)){
beta <- rep(beta,length(times))[1:length(times)]
}

if(length(births) < length(times)){
births <- rep(births,length(times))
}

for(n in 1:nsim){

S <- I <- rep(NA,length(times))
S[1] <- round(S0)
I[1] <- round(I0)

for(t in 2:length(times)){

lambda <- min(S[t-1],unname(beta[t-1] * S[t-1] * (I[t-1])^alpha))

if(stochastic){
I[t] <- rnbinom(n=1,mu=lambda,size=I[t-1]+1e-10)
}else{
I[t] <- lambda
}
S[t] <- max(S[t-1] + births[t-1] - I[t],0)

}
I.mat[,n] <- I
S.mat[,n] <- S
}

I.mat <- data.frame(I.mat)
S.mat <- data.frame(S.mat)

I.error <- apply(as.matrix(I.mat), 1, function(x){mean(x)+c(-1.96,1.96)*sd(x)/sqrt(length(x))})

I.mat\$high <- I.error[2,]
I.mat\$low <- I.error[1,]

S.error <- apply(as.matrix(S.mat), 1, function(x){mean(x)+c(-1.96,1.96)*sd(x)/sqrt(length(x))})

S.mat\$high <- S.error[2,]
S.mat\$low <- S.error[1,]

I.mat\$mean <- rowMeans(I.mat[,1:nsim],na.rm=T)
S.mat\$mean <- rowMeans(S.mat[,1:nsim],na.rm=T)

I.mat\$time <- times
S.mat\$time <- times

return(list('I'=I.mat,'S'=S.mat))
}

#' @title TSIR_LE
#' @description A function to calculate the Lyapunov Exponennt (LE) from the TSIR model
#' @param I The I output from the simulated or predicted TSIR model
#' @param S The S output from the simulated or predicted TSIR model
#' @param alpha The homogeneity parameter from the simulated or predicted TSIR model
#' @param beta The inferred contact rate, use beta = contact\$beta where contact is an output from runtsir or simulatetsir
#' @param time The time vector from the data or simulated data
#' @param IP The generation interval of the pathogen (in weeks)
#'
#' @examples
#'
#' \dontrun{
#' require(kernlab)
#' require(ggplot2)
#' require(kernlab)
#' London <- twentymeas\$London
#' ## just analyze the biennial portion of the data
#' London <- subset(London, time > 1950)
#'
#' ## define the interval to be 2 weeks
#' IP <- 2
#'
#'## first estimate paramters from the London data
#' parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log')
#'
#' ## look at beta and alpha estimate
#' plotbeta(parms)
#'
#' ## simulate the fitted parameters
#' sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2)
#'
#' ## now lets predict forward 200 years using the mean birth rate,
#' ## starting from rough initial conditions
#' times <- seq(1965,2165, by = 1/ (52/IP))
#' births <- rep(mean(London\$births),length(times))
#' S0 <- parms\$sbar
#' I0 <- 1e-5*mean(London\$pop)
#'
#' pred <- predicttsir(times=times,births=births,
#'                    beta=parms\$contact\$beta,alpha=parms\$alpha,
#'                   S0=S0,I0=I0,
#'                   nsim=50,stochastic=T)
#'
#' ## take the last 10 years
#' pred <- lapply(pred, function(x)  tail(x, 52/IP * 20) )
#'
#' ## now compute the Lyapunov Exponent for the simulate and predicted model
#'
#' simLE <- TSIR_LE(
#' time=sim\$res\$time,
#' S=sim\$simS\$mean,
#' I=sim\$res\$mean,
#' alpha=sim\$alpha,
#'   beta=sim\$contact\$beta,
#' IP=IP
#' )
#'
#' predLE <- TSIR_LE(
#' time=pred\$I\$time,
#' S=pred\$S\$X3,
#' I=pred\$I\$X3,
#' alpha=parms\$alpha,
#' beta=parms\$contact\$beta,
#' IP=IP
#' )
#'
#' simLE\$LE
#' predLE\$LE
#'
#'
#' }

TSIR_LE = function(time, S, I, alpha, beta, IP) {
IT <- length(I)
j11 <- rep(NA, IT)
j12 <- rep(NA, IT)
j21 <- rep(NA, IT)
j22 <- rep(NA, IT)
# initial unit vector
J <- matrix(c(1, 0), ncol = 1)
# loop over the attractor
for (i in 1:IT) {
j11[i] <- 1 - beta[i %% (52/IP) + 1] * I[i]^alpha
j12[i] <- - beta[i %% (52/IP) + 1] * S[i] * (I[i]^(alpha - 1) * alpha)
j21[i] <- beta[i %% (52/IP) + 1] * I[i]^alpha
j22[i] <- beta[i %% (52/IP) + 1] * S[i] * I[i]^(alpha - 1) * alpha
J <- matrix(c(j11[i], j12[i], j21[i], j22[i]), ncol = 2, byrow = TRUE) %*% J
if(is.nan(J[1])) stop(paste0('Matrix product has become NaN at index ',i))
}
res <- list(LE = log(norm(J))/IT, j11 = j11, j12 = j12, j21 = j21, j22 = j22,time=time,
I = I, S = S, alpha = alpha, beta = beta)
return(res)
}

#' @title TSIR_LLE
#' @description A function to calculate the Local Lyapunov Exponennt (LLE) from the TSIR model
#' @param LE The output of TSIR_LE to pass the Jacobian elements
#' @param m The window to sweep the time-varying Jacobian elements. Defaults to one.
#' @examples
#'
#'
#' \dontrun{
#' require(kernlab)
#' require(ggplot2)
#' require(kernlab)
#' London <- twentymeas\$London
#' ## just analyze the biennial portion of the data
#' London <- subset(London, time > 1950)
#'
#' ## define the interval to be 2 weeks
#' IP <- 2
#'
#'## first estimate paramters from the London data
#' parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log')
#'
#' ## look at beta and alpha estimate
#' plotbeta(parms)
#'
#' ## simulate the fitted parameters
#' sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2)
#'
#' ## now lets predict forward 200 years using the mean birth rate,
#' ## starting from rough initial conditions
#' times <- seq(1965,2165, by = 1/ (52/IP))
#' births <- rep(mean(London\$births),length(times))
#' S0 <- parms\$sbar
#' I0 <- 1e-5*mean(London\$pop)
#'
#' pred <- predicttsir(times=times,births=births,
#'                    beta=parms\$contact\$beta,alpha=parms\$alpha,
#'                   S0=S0,I0=I0,
#'                   nsim=50,stochastic=T)
#'
#' ## take the last 10 years
#' pred <- lapply(pred, function(x)  tail(x, 52/IP * 20) )
#'
#' ## now compute the Lyapunov Exponent for the simulate and predicted model
#'
#' simLE <- TSIR_LE(
#' time=sim\$res\$time,
#' S=sim\$simS\$mean,
#' I=sim\$res\$mean,
#' alpha=sim\$alpha,
#'   beta=sim\$contact\$beta,
#' IP=IP
#' )
#'
#' predLE <- TSIR_LE(
#' time=pred\$I\$time,
#' S=pred\$S\$X3,
#' I=pred\$I\$X3,
#' alpha=parms\$alpha,
#' beta=parms\$contact\$beta,
#' IP=IP
#' )
#'
#' simLE\$LE
#' predLE\$LE
#'
#' simLLE <- TSIR_LLE(simLE)
#' predLLE <- TSIR_LLE(predLE)
#'
#' plotLLE(simLLE)
#' plotLLE(predLLE)
#' }
#'

TSIR_LLE = function(LE, m = 1) {
LLE <- rep(NA, length(LE\$I))
for (i in 1:(length(LE\$I) - m)) {
J <- matrix(c(1, 0, 0, 1), ncol = 2)
for (k in 0:(m - 1)) {
J <- matrix(c(LE\$j11[(i + k)], LE\$j12[(i + k)], LE\$j21[(i + k)], LE\$j22[(i + k)]), ncol = 2, byrow = TRUE) %*% J
}
LLE[i] = log(max(abs(eigen(J)\$values)))/m
}
res <- list(LLE = LLE, time=LE\$time, I = LE\$I, S = LE\$S)
return(res)
}

#' @title plotLLE
#' @description Function to plot the Local Lyapunov Exponents. The output is of class ggplot2 so you can add standard
#' ggplot2 options to it if desired.
#' @param LLE The output from TSIR_LLE
#' @export
#' @examples
#'
#'
#' \dontrun{
#' require(kernlab)
#' require(ggplot2)
#' require(kernlab)
#' London <- twentymeas\$London
#' ## just analyze the biennial portion of the data
#' London <- subset(London, time > 1950)
#'
#' ## define the interval to be 2 weeks
#' IP <- 2
#'
#'## first estimate paramters from the London data
#' parms <- estpars(data=London, IP=2, regtype='gaussian',family='poisson',link='log')
#'
#' ## look at beta and alpha estimate
#' plotbeta(parms)
#'
#' ## simulate the fitted parameters
#' sim <- simulatetsir(data=London,parms=parms,IP=2,method='deterministic',nsim=2)
#'
#' ## now lets predict forward 200 years using the mean birth rate,
#' ## starting from rough initial conditions
#' times <- seq(1965,2165, by = 1/ (52/IP))
#' births <- rep(mean(London\$births),length(times))
#' S0 <- parms\$sbar
#' I0 <- 1e-5*mean(London\$pop)
#'
#' pred <- predicttsir(times=times,births=births,
#'                    beta=parms\$contact\$beta,alpha=parms\$alpha,
#'                   S0=S0,I0=I0,
#'                   nsim=50,stochastic=T)
#'
#' ## take the last 10 years
#' pred <- lapply(pred, function(x)  tail(x, 52/IP * 20) )
#'
#' ## now compute the Lyapunov Exponent for the simulate and predicted model
#'
#' simLE <- TSIR_LE(
#' time=sim\$res\$time,
#' S=sim\$simS\$mean,
#' I=sim\$res\$mean,
#' alpha=sim\$alpha,
#'   beta=sim\$contact\$beta,
#' IP=IP
#' )
#'
#' predLE <- TSIR_LE(
#' time=pred\$I\$time,
#' S=pred\$S\$X3,
#' I=pred\$I\$X3,
#' alpha=parms\$alpha,
#' beta=parms\$contact\$beta,
#' IP=IP
#' )
#'
#' simLE\$LE
#' predLE\$LE
#'
#' simLLE <- TSIR_LLE(simLE)
#' predLLE <- TSIR_LLE(predLE)
#'
#' plotLLE(simLLE)
#' plotLLE(predLLE)
#' }
#'

plotLLE <- function(LLE){

local.ly <- LLE\$LLE

pm <- LLE\$LLE > 0

S <- LLE\$S

I <- LLE\$I

time <- LLE\$time

LLdf <- data.frame(time,S,I,'LLEsign'=pm,'LLE'=local.ly)

ggplot(LLdf,aes(time,I,colour=LLE,group=1))+geom_line()

p1 <- ggplot(LLdf,aes_string('S','I',size='abs(LLE)',shape='LLEsign',colour='LLEsign'))+geom_point(alpha=0.7)+ theme_bw()+
theme(text=element_text(size=12))+
scale_size_continuous(limits = c(0,1),range = c(2,10),name = 'abs(LLE)')+
scale_shape_manual(values=c(1,19),labels = c("LLE < 0", "LLE > 0"),name='sign(LLE)')+
# labs(size="abs(LLE)",
#      col='sign(LLE)', shape='sign(LLE)')+
#scale_color_manual(values=c('grey2','grey'))+
scale_color_manual(values=c('orangered3','dodgerblue'),labels = c("LLE < 0", "LLE > 0"),name='sign(LLE)')
#guides(shape=FALSE,size=FALSE,colour=FALSE)

return(p1)
}
```

## Try the tsiR package in your browser

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

tsiR documentation built on Jan. 21, 2021, 1:06 a.m.