## INLA Code for model used in Parks et al. Nature Medicine 2020


## year.month, year.month2, year.month3 - month and year over time, values 1-456



## month1, month2, month3, month4, month5, month6, month7 - indexes for month, values 1-12


## ID1, ID2 - indexes for state  

## USA.adj - adjacency matrix for state units of contiguous USA


## anomaly - temperature anomaly 

# hyperparameter value
hyper_value = 0.001



#INLA formula



        fml  <- deaths ~

        # global terms
        	1 +                                                                     			# global intercept
        	year.month +                                                           				# global slope

        # month specific terms
        	f(month, model='rw1',cyclic = TRUE, hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))) +          		# month specific intercept
        	f(month2, year.month2, model='rw1', cyclic= TRUE, hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))) +          # month specific slope

        # state specific terms
        	f(ID, model="bym",graph=USA.adj, hyper = list(prec.unstruct = list(prior = "loggamma", param = c(1, hyper_value),fixed=FALSE),
        	  prec.spatial = list(prior = "loggamma", param = c(1, hyper_value),fixed=FALSE))) +                     # state specific intercept (BYM)
        	f(ID2, year.month2, model="bym",graph=USA.adj, hyper = list(prec.unstruct = list(prior = "loggamma", param = c(1, hyper_value),fixed=FALSE),
        	  prec.spatial = list(prior = "loggamma", param = c(1, hyper_value),fixed=FALSE))) +                  	 # state specific slope (BYM)

        # state-month specific terms
        	f(month3, model="rw1",cyclic = TRUE,group=ID,control.group=list(model='besag',graph=USA.adj, hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))),
        	  hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value),fixed=FALSE))) +            	# state-month specific intercept (spatially-correlated)
        	f(month4, year.month2, model="rw1",cyclic = TRUE,group=ID, control.group=list(model='besag',graph=USA.adj, hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))),
        	  hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value),fixed=FALSE))) +            	# state-month specific slope (spatially-correlated)
        	f(month6, model="rw1",cyclic = TRUE,group=ID,control.group=list(model='iid', hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))),
        	  hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))) +                           # state-month specific intercept (spatially-correlated)
        	f(month7, year.month2, model="rw1",cyclic = TRUE,group=ID, control.group=list(model='iid', hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))),
        	  hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))) +                    	# state-month specific slope (spatially-correlated)

        # temperature specific terms
        	f(month5, anomaly, model="rw1", cyclic=TRUE, hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))) +            	# month specific temperature slope

        # random walk across time
        	f(year.month3, model="rw1", hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value)))) +                           	# rw1
        
	# overdispersion term
        	f(e, model = "iid", hyper = list(prec = list(prior = "loggamma", param = c(1, hyper_value))))  	# overdispersion term



###############################################