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)
  122.   #Add names
  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){
  141.     add=simforeward(dist=dist,mass=pop,cross=cross,beta=unlist(beta[i,]),start=start,first=(first+1),last=last)
  142.     ret[,(i+1)] <- add
  143.     names[i+1] <- paste('Sim',i,sep='')
  144.     print(i)
  145.   }
  146.   colnames(ret)<-names
  147.   return(ret)
  148. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top