Thanks, sorry I should include the code I am trying to 'Shiny'! Please see below. The model outputs two different plots that are updated as the model runs which is simulating the passage of time. I am trying to convert this to a shiny app so the user can adjust a couple of the parameters and then see the impact they have on disease spread. I have got some simpler versions of the code to work but I am struggling to get the code to display the plots as they transition through time... it seems that Shiny runs the code and allows the model to finish and then just plots the final situation. I am experimenting with the reactive elements but not quite sure if this is something that can be done... can the plot continually update as the model runs?? I have tried to neaten up the code to a series of functions and then one which 'steps' the process forward one time step but I am struggling to get it to plot one stage at a time and continually update as such. The code I have added runs fine as stand alone... outside of shiny. But if anyone could point me in the right direction as to what functions/commands in Shiny I need to be using I would be very grateful. Many thanks.
########Create and plot associations of Simuland citizens who choose social contacts
########according to degree (number of citizens contacted per day) and loyalty (scales
########the size of the neighborhood explored by each citizen)
rm(list=ls())
set.seed<-1
library(igraph)
par(mfrow=c(2,1))
businessarea=0.1
degree<-5
loyalty<-0.9
#set colours for plotting
cols<-c("white","yellow","orange","red","purple","cyan","black")
pop_set_up<-function(family_size=5,households=100,xMin=-.5,xMax=.5,yMin=-.5,yMax=.5,businessarea){
pop<-family_size*households
#Extended simulation windows parameters
rExt=businessarea; #extension parameter -- use cluster radius
xMinExt=xMin-rExt;
xMaxExt=xMax+rExt;
yMinExt=yMin-rExt;
yMaxExt=yMax+rExt;
#rectangle dimensions
xDeltaExt=xMaxExt-xMinExt;
yDeltaExt=yMaxExt-yMinExt;
areaTotalExt=xDeltaExt*yDeltaExt; #area of extended rectangle
#homestead positions
#x and y (uniform) coordinates of Poisson points for the parent
xParent=xMinExt+xDeltaExt*runif(households);
yParent=yMinExt+yDeltaExt*runif(households);
#replicate for daughters
#replicate parent points (ie centres of disks/clusters)
xP=rep(xParent,each=family_size);
yP=rep(yParent,each=family_size);
###from epidemic model...setup
#age the citizens
age<-rbinom(length(xP),1,0.25)
output<-list(pop,xP,yP,age,households,xMin,xMax,yMin,yMax)
names(output)<-c("pop","xP","yP","age","households","xMin","xMax","yMin","yMax")
return(output)
}
disease_set_up<-function(pop){
#start with 2 infected individuals
expa<-sample(1:pop,2,replace=FALSE)
#a necessary disease parameter
E_I1<-5
#create dataframe to store disease state
S<-rep(1,pop)
E<-rep(0,pop)
I1<-rep(0,pop)
I2<-rep(0,pop)
I3<-rep(0,pop)
R<-rep(0,pop)
D<-rep(0,pop)
status<-data.frame(S,E,I1,I2,I3,R,D)
status$S[expa]<-0
status$E[expa]<-1
d_exp<-rep(NA,pop)
d_inf1<-rep(NA,pop)
d_inf2<-rep(NA,pop)
d_inf3<-rep(NA,pop)
d_exp[status$E==1]<-rpois(1,E_I1)
#create matrix to store epidemic progress
progression<-matrix(0,nr=1,nc=ncol(status))
progression[1,]<-colSums(status)
output<-list(status,progression,d_exp,d_inf1,d_inf2,d_inf3)
names(output)<-c("status","progression","d_exp","d_inf1","d_inf2","d_inf3")
return(output)
}
disease_model<-function(Spop=pop$pop,mat,age=pop$age,status=dis$status,d_exp=dis$d_exp,d_inf1=dis$d_inf1,d_inf2=dis$d_inf2,d_inf3=dis$d_inf3,progression=dis$progression){
#Probability of becoming exposed having contacted an infectious individual (daily)
#Need to work out R0 based on other parameters
S_E<-0.1
#lambda for a Poisson draw for the length of this period
E_I1<-5
#probability of transitioning to serious case for young (daily)
yI1_I2<-0.01
#and same for old
oI1_I2<-0.1
#diff between young and old
diffI1_I2<-oI1_I2-yI1_I2
#probability of transitioning to critical case for young (daily)
yI2_I3<-0.05
#and same for old
oI2_I3<-0.2
#diff between young and old
diffI2_I3<-oI2_I3-yI2_I3
#probability of death for young (daily)
yI3_D<-0.1
#and same for old
oI3_D<-0.4
#diff between young and old
diffI3_D<-oI3_D-yI3_D
#lambda for a Poisson draw for duration of a mild infection
yI1_R<-7
#and same for old
oI1_R<-7
#diff between young and old
diffI1_R<-oI1_R-yI1_R
#lambda for a Poisson draw for duration of a hospitalised infection
yI2_R<-5
#and same for old
oI2_R<-5
#diff between young and old
diffI2_R<-oI2_R-yI2_R
#lambda for a Poisson draw for duration of a critical infection
yI3_R<-5
#and same for old
oI3_R<-5
#diff between young and old
diffI3_R<-oI3_R-yI3_R
#this will block individuals changing twice in the same time-step
f2c<-rep(1,Spop)
inf<-sign(status$I1+status$I2+status$I3)
risk<-rep(0,Spop)
risk[as.numeric(names(table(which(inf*mat==1,arr.ind=TRUE)[,2])))]<-table(which(inf*mat==1,arr.ind=TRUE)[,2])
prob<-(1-S_E)^risk
infect<-rbinom(Spop,1,1-prob)*status$S
status$S[which(infect==1)]<-0
status$E[which(infect==1)]<-1
d_exp[!is.na(d_exp)&d_exp>0]<-d_exp[!is.na(d_exp)&d_exp>0]-1
d_exp[which(infect==1)]<-rpois(length(which(infect==1)),E_I1)
f2c[which(infect==1)]<-0
EI1<-which(status$E==1&!is.na(d_exp)&d_exp==0)
status$E[EI1]<-0
status$I1[EI1]<-1
f2c[EI1]<-0
I1I2<-rbinom(Spop,1,yI1_I2+age*diffI1_I2)*status$I1*f2c
status$I1[which(I1I2==1)]<-0
status$I2[which(I1I2==1)]<-1
f2c[which(I1I2==1)]<-0
I2I3<-rbinom(Spop,1,yI2_I3+age*diffI2_I3)*status$I2*f2c
status$I2[which(I2I3==1)]<-0
status$I3[which(I2I3==1)]<-1
f2c[which(I2I3==1)]<-0
I3D<-rbinom(Spop,1,yI3_D+age*diffI3_D)*status$I3*f2c
just.died<-which(I3D==1)
status$I3[just.died]<-0
status$D[just.died]<-1
f2c[just.died]<-0
EI1y<-EI1[EI1%in%which(age==0)]
EI1o<-EI1[EI1%in%which(age==1)]
d_inf1[!is.na(d_inf1)&d_inf1>0]<-d_inf1[!is.na(d_inf1)&d_inf1>0]-1
if(length(EI1y)>0){
d_inf1[EI1y]<-rpois(length(EI1y),yI1_R)
}
if(length(EI1o)>0){
d_inf1[EI1o]<-rpois(length(EI1o),oI1_R)
}
d_inf2[!is.na(d_inf2)&d_inf2>0]<-d_inf2[!is.na(d_inf2)&d_inf2>0]-1
d_inf2[which(I1I2==1&age==0)]<-rpois(length(which(I1I2==1&age==0)),yI2_R)
d_inf2[which(I1I2==1&age==1)]<-rpois(length(which(I1I2==1&age==1)),oI2_R)
d_inf3[!is.na(d_inf3)&d_inf3>0]<-d_inf3[!is.na(d_inf3)&d_inf3>0]-1
d_inf3[which(I2I3==1&age==0)]<-rpois(length(which(I2I3==1&age==0)),yI3_R)
d_inf3[which(I2I3==1&age==1)]<-rpois(length(which(I2I3==1&age==1)),oI3_R)
I1R<-which(status$I1==1&!is.na(d_inf1)&d_inf1==0)
status$I1[I1R]<-0
status$R[I1R]<-1
f2c[I1R]<-0
I2R<-which(status$I2==1&!is.na(d_inf2)&d_inf2==0)
status$I2[I2R]<-0
status$R[I2R]<-1
f2c[I2R]<-0
I3R<-which(status$I3==1&!is.na(d_inf3)&d_inf3==0)
status$I3[I3R]<-0
status$R[I3R]<-1
f2c[I3R]<-0
progression<-rbind(progression,colSums(status))
status.now<-as.matrix(status)%*%as.matrix(seq(1,7))
output<-list(status,status.now,progression,d_exp,d_inf1,d_inf2,d_inf3,just.died)
names(output)<-c("status","status.now","progression","d_exp","d_inf1","d_inf2","d_inf3","just.died")
return(output)
}
night_net<-function(pop=pop,households=households){
##########Create association matrix (Night)...will stay constant unless Sick or dead go to hospital
#set up association matrix
ass.night<-array(0,dim=c(pop,pop))
#a useful index vector
assseq<-seq(1,pop)
#make association matrix
#I feel like this could be speeded up, avoiding a for-loop
#just want bliocks of 1s for housemates
for(i in 1:households){
ass.night[(i*5-4):(i*5),(i*5-4):(i*5)]<-1}
return(ass.night)
}
forward<-function(Spop=pop$pop,businessarea,xP=pop$xP,yP=pop$yP,loyalty,degree,ass.night,status){
#####DAY
#Generate the (relative) locations in polar coordinates by
#simulating independent variables.
theta=2*pi*runif(Spop); #angular coordinates
rho=businessarea*sqrt(runif(Spop)); #radial coordinates
#Convert from polar to Cartesian coordinates
x0=rho*cos(theta);
y0=rho*sin(theta);
#translate points (ie parents points are the centres of cluster disks)
x=xP+x0;
y=yP+y0;
##########Bit of data handling
citizens<-data.frame(x,y)
##########Create association matrix (Business)
#calculate distance matrix
dm<-as.matrix(dist(citizens,method="euclidean"))
#set up (empty) association matrix
ass.day<-dm*0
###here we loop through each citizen, have them select "degree" partners, and associate them symmetrically.
###this is not the strict definition of "degree", but it does seem to describe human behaviour ("I seek n contacts per day")
#loop through citizens
for(assrow in 1:Spop){
dist_t<-dm[assrow,]
possible.partners<-seq(1,Spop)
possible.partners<-possible.partners[dist_t<=quantile(dist_t,1-loyalty)&dist_t>0]
partners<-sample(possible.partners,size=degree,replace=F)
ass.day[assrow,partners]<-ass.day[assrow,partners]+1
}
#ass.day<-sapply(1:nrow(ass.day),assign_contacts)
ass.day<-sign(ass.day+t(ass.day))
##MJS ADDITION##
#Disconnect individuals that have severe infection from day time network
ass.day[status$I2==1,]<-0
ass.day[,status$I2==1]<-0
ass.day[status$I3==1,]<-0
ass.day[,status$I3==1]<-0
#Dave query....doi we want to do something interesting with them like leave them at home so that folk could still interact with them?
#or have them iunteract with health professionals, hence affecting the health service carrying capacity?
#this bit makes an association network for plotting,removing any dead individs
business.ass<-ass.day
business.ass[status$D==1,]<-0
business.ass[,status$D==1]<-0
mat<-ass.night+ass.day
#simplify the total association matrix
mat<-sign(mat)
#remove any deaths from the total association matrix
mat[status$D==1,]<-0
mat[,status$D==1]<-0
dis<-disease_model(mat=mat)
output<-list(business.ass,mat,dis,citizens)
names(output)<-c("business.ass","mat","dis","citizens")
return(output)
}
####################################
####################################
####################################
pop<-pop_set_up(businessarea=businessarea)
dis<-disease_set_up(pop=pop$pop)
an<-night_net(pop=pop$pop,households=pop$households)
status<-dis$status
counter<-0
t<-0
while((counter<10)&(t<200)){
t<-t+1
step<-forward(businessarea=businessarea,degree=degree,loyalty=loyalty,ass.night=an,status=status)
dis<-step$dis
status<-dis$status
if(sum(status$R)+sum(status$D)==nrow(status)){counter<-counter+1}
print(dis$progression[nrow(dis$progression),])
#plot and add association lines
#cool code finds the array indices of the interacting dyads
dd<-which(step$business.ass>0,arr.ind=T)
#set up array of arrow vertices
arrow.x1<-arrow.x2<-arrow.y1<-arrow.y2<-numeric(dim(dd)[1])
#and for each dyad, populate the vertices. Might be able to avoid a for-loop here?
for(i in 1:dim(dd)[1]){
arrow.x1[i]<-step$citizens$x[dd[i,1]]
arrow.y1[i]<-step$citizens$y[dd[i,1]]
arrow.x2[i]<-step$citizens$x[dd[i,2]]
arrow.y2[i]<-step$citizens$y[dd[i,2]]}
#and for plotting purposes, checks whether alive
is.alive<-dis$status$D==0
is.ill<-(dis$status$I2==1|dis$status$I3==1)
step$x[is.ill]<-pop$xP[is.ill]
step$y[is.ill]<-pop$yP[is.ill]
dev.hold()
plot(NULL,xlim=c(pop$xMin,pop$xMax),ylim=c(pop$yMin,pop$yMax),xlab="latitude",ylab="longitude",pch=21,bg=cols[(dis$status.now[is.alive])],col="black",cex=1.5)
arrows(arrow.x1,arrow.y1,arrow.x2,arrow.y2,code=3,length=0,col="gray80")
points(step$citizens$x[is.alive],step$citizens$y[is.alive],pch=21,bg=cols[(dis$status.now[is.alive])],col="black",cex=1.5)
points(pop$xP[dis$just.died],pop$yP[dis$just.died],pch=21,bg="black",col="black",cex=1.5)
severe.cases<-(dis$progression[,4]+dis$progression[,5])/pop$pop
mort<-dis$progression[,7]/pop$pop
Time<-seq(1,201)
plot(NULL,xlab="Time",ylab="Cumulative mortality",ylim=c(0,0.1),xlim=c(1,Time[t+1]),type="n",las=1)
lines(mort[1:(t+1)]~Time[1:(t+1)],lwd=2)
lines(severe.cases[1:(t+1)]~Time[1:(t+1)],ylab="Severity",lwd=2,col="red")
dev.flush()
}
