Simulate_contact_model: Simulation of the epidemic process
In holaanna/contactsimulator: Contact model simulator

Description Usage Arguments Details Value References See Also Examples

Epidemic simulation using the contact type model.

1
2
3

Simulate_contact_model(param, grid_lines, pop_grid, grid_size = 5000,
  age_level = c(1, 1), age_dist = c(1, 0), m_start = 1,
  t_max = 118, t_intervention = 365, EI_model = 1, kern_model = 4)

`param`	Indicating a data frame containing a vector of parameters including: epsion The primary infection rate. See `func_time_beta` beta_0 Baseline or average transmission rate. See `func_time_beta` beta_1 Amplitude of the seasonality. See `func_time_beta` alpha1,alpha2 The dispersal kernel parameters. mu_lat,var_lat mean and variance of the latent period. See `E_to_I` for details. t0 Time at which the primary source became active . omega Period of the forcing. See `func_time_beta` gama The mean proportion of short range dispersal events. .
`pop_grid`	Population density of the grid a case resides. This is filled from bottom to top, then left to right.
`grid_size`	Grid resolution //@inheritParams circle_line_intersections
`age_level, age_dist`	Vectors of age level and the propportion of each age group respectively. See details.
`m_start`	The size of initial cases. Default is 1.
`t_max`	Final observation time.
`t_intervention`	Start of the intervention if any.
`EI_model`	Take integer values to specify the type of model used for the latent period. See `E_to_I`
`kern_model`	Take integer values to specify the type of dispersal kernel used. See `Samp_dis`

Simulate_contact_model provide the simulation of the epidemic process and the maximum distance the wave can travel at each particular obaservation date.

The contact model developped here only consider a short range dispersal kernal as BBTV (the pathogen considered) only spead at short range from the source (ref). Consequentely, we use an expoential dispersal parameter α of the form:

f(r,α)=α exp(-α r)

A data frame with components:

k: Index of the case
coor_x,coor_y: Location of the infection premisse
t_e,ti,tr: Exposure, infection and removal times respectively of the new premisses.
age: Age group category where the new premisse belongs to.
infected_source: Source of infection. 9999 for primary infection.
value_indx: Grid index in which lies the infection premisse in the raster.

\insertRef

KR08contactsimulator \insertRefMee11contactsimulator

E_to_I, func_time_beta.

data(bbtv)
attach(bbtv)
Dat<- bbtv[,c("longitude","latitude","BBTV","inspectiondate","leavesinfected","treatmentdate","location")]
Dat1<-subset(Dat,Dat$latitude> -27.4698 & Dat$BBTV%in%c("P&I","P", "NI") & difftime(as.Date(Dat$inspectiondate), as.Date("2010/01/01"), unit="days")>=0)  # data up in queensland
Dat1$treatmentdate[is.na(Dat1$treatmentdate)]<- Dat1$inspectiondate[is.na(Dat1$treatmentdate)]
Dat1$detection<-as.numeric(difftime(as.Date(Dat1$inspectiondate), as.Date("2010/01/01"), unit="days"))
Dat1$removal<-as.numeric(difftime(as.Date(Dat1$treatmentdate), as.Date("2010/01/01"), unit="days"))
Dat1$removal[which(Dat1$removal<0)]<- Dat1$detection[which(Dat1$removal<0)]
Datt<-Dat1[,c("longitude","latitude","BBTV","leavesinfected","detection","removal")]


Datt[which(Datt$leavesinfected=="LOTS"),"leavesinfected"]<- 45
Datt[which(Datt$leavesinfected=="1,2,4"),"leavesinfected"]<- 2.3
Datt[which(Datt$leavesinfected=="'3"),"leavesinfected"]<- 3
Datt[which(Datt$leavesinfected=="2 +bunch"),"leavesinfected"]<- 2
Datt[which(Datt$leavesinfected=="3 +bunch"),"leavesinfected"]<- 3
Datt[which(Datt$leavesinfected=="4+BUNCH"),"leavesinfected"]<- 4
Datt[which(Datt$leavesinfected=="avg 3.2"),"leavesinfected"]<- 3.2
Datt[which(Datt$leavesinfected=="1-6, avg 3.5"),"leavesinfected"]<- 3.5
Datt[which(Datt$leavesinfected=="all"),"leavesinfected"]<- 45


leav=sapply(Datt[,"leavesinfected"],function(x){
  gsub("all/","",x)
})

leav=sapply(leav,function(x){
  gsub("/all","",x)
})

leav[grepl("[+]",leav)]<- 45  # Assuming 45 leaves on a plant

Datt$leavesinfected<- leav

Datt=Datt[with(Datt,order(Datt$detection)),]

# Australian reference system
sp::coordinates(Datt) <- c("longitude", "latitude")
sp::proj4string(Datt) <- sp::CRS("+init=epsg:4326")
australianCRS <- sp::CRS("+init=epsg:3577")

pointsinaustraliangrid = sp::spTransform(Datt,australianCRS)

# Raster
rast <- raster::raster()
raster::extent(rast) <- raster::extent(pointsinaustraliangrid) # Set same extent

raster::res(rast)=5000 # Set resolution

 size<- raster::res(rast)[1]
dif=(raster::xmax(pointsinaustraliangrid)-raster::xmin(pointsinaustraliangrid))/size
cei= abs(ceiling(dif))

if(cei!=dif){
  if(raster::xmax(rast)!=raster::xmax(pointsinaustraliangrid)){
    raster::xmax(rast)<- raster::xmin(rast) + size*cei
  }
  if(raster::xmin(rast)!=raster::xmin(pointsinaustraliangrid)){
    raster::xmin(rast)<- raster::xmax(rast) - size*cei
  }

}

# Adding row at the top or bottom of the grid if raster leaves points out

dif1=(raster::ymax(pointsinaustraliangrid)-raster::ymin(pointsinaustraliangrid))/size
cei1= abs(ceiling(dif1))


if(cei1!=dif1){
  if(raster::ymax(rast)!=raster::ymax(pointsinaustraliangrid)){
    raster::ymax(rast)<- raster::ymin(rast) + size*cei1
  }
  if(raster::ymin(rast)!=raster::ymin(pointsinaustraliangrid)){
    raster::ymin(rast)<- raster::ymax(rast) - size*cei1
  }

}

raster::res(rast)=5000 # Set resolution

# And then ... rasterize it! This creates a grid version
# of your points using the cells of rast,

rast2 <- raster::rasterize(pointsinaustraliangrid, rast, 1, fun=sum)

# Extract infos on the grid


n_row_grid=nrow_grid=raster::nrow(rast)
n_col_grid=ncol_grid=raster::ncol(rast)
grid_size=raster::res(rast)[1]     # Resolution

n_line=(nrow_grid+1) + (ncol_grid +1)  # Number of grid  lines

x_min=raster::xmin(rast)  # min max of the bounding box
x_max=raster::xmax(rast)

y_min=raster::ymin(rast)
y_max=raster::ymax(rast)

da=as.data.frame(pointsinaustraliangrid)

pop_per_grid=raster::values(rast2)
pop_per_grid[is.na(pop_per_grid)]=0
mat=matrix(pop_per_grid,nrow = nrow_grid, byrow = TRUE)
pop_grid=apply(mat,2,rev)     # population per grid

# Structure of the grid
x=seq(x_min,x_max,grid_size)
y=seq(y_min,y_max,grid_size)

grid_lines=array(0,c(n_line,6))
for(i in 1:n_line){
  if(i<=(nrow_grid +1)){
    grid_lines[i,]=c(i,1,x[1],y[i],x[length(x)],y[i])
  }
  else{
    grid_lines[i,]=c(i,2,x[i-length(y)],y[1],x[i-length(y)],y[length(y)])
  }
}

grid_lines=as.data.frame(grid_lines)
colnames(grid_lines)<- c("indx","orient_line","coor_x_1","coor_y_1","coor_x_2","coor_y_2")
circle_x=2022230
circle_y=-3123109
r=10000
# Simulation with exponential kernel
alpha<- 30; beta<- 0.012; epsilon<- 0.02; omega<- 0.12; mu_lat<- 30; var_lat<- 20; t0<- 0; c<- 20; b1<-0; gama<- 0.5
param=data.frame(alpha1=alpha, alpha2=alpha, beta=beta, epsilon=epsilon, omega=omega, mu_lat=mu_lat, var_lat=var_lat, t0=t0, c=c, b1=b1, gama=gama)

Simulate_contact_model(param, grid_lines, pop_grid)

detach(bbtv)