• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# ebola-GravNet-border.R

a guest Sep 5th, 2019 87 in 139 days
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. library(Matrix)
2.
3. #Spatial decay function
4. decay<-function(x,mass,cross,beta0,beta1,beta2,beta3){
5.   return(c(1,beta3)[cross+1]/(1+exp(beta0 + (beta1*x)/(mass**beta2))))
6. }
7.
8. makecross <- function(regions){
9.   return(outer(regions, regions, get("!=")))
10. }
11.
12. makecross.group <- function(regions, group=c(5, 8, 13)){
13.   v <- regions %in% group
14.   return(outer(v, v, get("!=")))
15. }
16.
17. #returns probabilities that each node will be infected in the next time step
18. getprobs<-function(start,dist,mass,cross,b){
19.
20.   #Check Trivial Cases
21.   if(all(start)){
22.     return(rep(1,length(start)))
23.   }
24.   if(!any(start)){
25.     return(rep(0,length(start)))
26.   }
27.
28.   #calculate probability of infection from each infected node node to each uninfected node
29.   mass=mass%*%t(mass)
30.   tmp=decay(dist[which(start),which(!start)],mass[which(start),which(!start)],cross[which(start),which(!start)],b[1],b[2],b[3],b[4])
31.   tmp=tmp*(dist[which(start),which(!start)] != 0)
32.
33.   #combine probabilities of infected nodes for each uninfected node
34.   if(!is.null(dim(tmp))){
35.     tmp=log(1-tmp)
36.     tmp=rep(1,dim(tmp)[1]) %*% tmp
37.     tmp=1-exp(tmp[1,])
38.   }else if(sum(!start)==1){
39.     tmp=1-exp(sum(log(1-tmp)))
40.   }
41.
42.   #construct return vector
43.   ret=vector(length=length(start))
44.   ret[which(start)]=1
45.   ret[which(!start)]=tmp
46.
47.   return(ret)
48. }
49.
50. #calculate log likelihood for a single timestep
51. lglike<-function(p,d){
52.   for(i in 1:length(d)){
53.     if(!d[i]){
54.       p[i]=1-p[i]
55.     }
56.   }
57.   return(sum(log(p)))
58. }
59.
60. #simulates infection spread over multiple time steps from a given starting point
61. simforeward<-function(dist, mass, cross, beta, start, first, last){
62.   ret=start
63.   for(i in first:last){
64.     #if all nodes are infected stop the simulation
65.     if(all(ret!=Inf)){
66.       return(ret)
67.     }
68.
69.     #calculate probabilities of infection
70.     p=getprobs((ret<i),dist,mass,cross,beta)
71.     for(j in 1:length(ret)){
72.       #infect nodes with calculated probabilities
73.       if(start[j] == Inf & p[j] > runif(1,0,1)){
74.         ret[j] = i
75.       }
76.     }
77.     start=ret
78.   }
79.   return(ret)
80. }
81.
82. #creates a sparse distance matrix
83. makedist<-function(coords, cutoff=Inf){
84.   dmat=sqrt(((coords[1,] %*% t(rep(1,length(coords[1,])))) - (rep(1,length(coords[1,])) %*% t(coords[1,])))**2 + ((coords[2,] %*% t(rep(1,length(coords[2,])))) - (rep(1,length(coords[2,])) %*% t(coords[2,])))**2)
85.
86.   dmat=dmat*(dmat<cutoff)
87.
88.   return(Matrix(dmat))
89. }
90.
91. #calculates the negative log likelihood over all time steps (objective function for optimization)
92. objf<-function(beta, dist, mass, cross, data, days){
93.   ret=0
94.
95.   for(i in 1:(days)){
96.     p=getprobs((data<i),dist,mass,cross,beta)
97.     ret=ret+lglike(p,(data<(i+1)))
98.   }
99.
100.   return(-ret)
101. }
102.
103. #find the best fit for beta by minimizing negative log likelihood
104. getbeta<-function(dist,mass,cross,data,start=c(5,100,.1,.1),days=max(data[which(data!=Inf)])){
105.   ret=optim(par=start,fn=objf,dist=dist, mass=mass, cross=cross, data=data, days=days)
106.   return(ret)
107. }
108.
109. ##
110. getdelta<-function(beta, dist, mass,cross,  data, days=max(data[which(data!=Inf)])){
111.   require(mvtnorm)
112.   require(numDeriv)
113.
114.   #Estimate Hessian at mle
115.   hess=hessian(objf, beta, dist=dist, mass=mass, cross=cross, data=data, days=days)
116.   #Invert Hessian to get Variance-Covariance Matrix
117.   A=solve(hess,diag(rep(1,length(beta))))
118.   #Calculate the 95th equicoordinate quantile
119.   q=qmvnorm(.95,mean=beta,sigma=A,tail='both')\$quantile
120.   #Get bounds
121.   ret=q*sqrt(diag(A))/length(beta)
123.   names=vector()
124.   for(i in 0:(length(beta)-1)){
125.     names=c(names,paste('delta',i,sep=''))
126.   }
127.   names(ret)<-names
128.   return(ret)
129. }
130.
131. #Simulate multiple runs
132. makesims<-function(data, PID, pop, cross, dist, beta, first=0, last=157, nsims=1000, seed=802){
133.   set.seed(seed)
134.   ret <- matrix(NA, nrow=length(PID), ncol=nsims+1)
135.   ret[,1]=PID
136.   names=rep(NA,(nsims+1))
137.   names[1]='PID'
138.   start=data
139.   start[which(start>first)]<-Inf
140.   for(i in 1:nsims){