Fitting Conditional and Simultaneous Autoregressive Spatial Models in hglm

We present a new version (> 2.0) of the hglm package for fitting hierarchical generalized linear models (HGLMs) with spatially correlated random effects. CAR() and SAR() families for conditional and simultaneous autoregressive random effects were implemented. Eigen decomposition of the matrix describing the spatial structure (e.g., the neighborhood matrix) was used to transform the CAR/SAR random effects into an independent, but heteroscedastic, Gaussian random effect. A linear predictor is fitted for the random effect variance to estimate the parameters in the CAR and SAR models. This gives a computationally efficient algorithm for moderately sized problems.


Introduction
We present an algorithm for fitting spatial generalized linear models with conditional and simultaneous autoregressive (CAR & SAR; Besag, 1974;Cressie, 1993) random effects.The algorithm completely avoids the need to differentiate the spatial correlation matrix by transforming the model using an eigen decomposition of the precision matrix.This enables the use of already existing methods for hierarchical generalized linear models (HGLMs; Lee and Nelder, 1996) and this algorithm is implemented in the latest version ( 2.0) of the hglm package (Rönnegård et al., 2010).
The hglm package, up to version 1.2-8, provides the functionality for fitting HGLMs with uncorrelated random effects using an extended quasi likelihood method (EQL; Lee et al., 2006).The new version includes a first order correction of the fixed effects based on the current EQL fitting algorithm, which is more precise than EQL for models having non-normal outcomes (Lee and Lee, 2012).Similar to the h-likelihood correction of Noh and Lee (2007), it corrects the estimates of the fixed effects and thereby also reduces potential bias in the estimates of the dispersion parameters for the random effects (variance components).The improvement in terms of reduced bias is substantial for models with a large number of levels in the random effect (Noh and Lee, 2007), which is often the case for spatial generalized linear mixed models (GLMMs).
Earlier versions of the package allow modeling of the dispersion parameter of the conditional mean model, with fixed effects, but not the dispersion parameter(s) of the random effects.The current implementation, however, enables the user to specify a linear predictor for each dispersion parameter of the random effects.Adding this option was a natural extension to the package because it allowed implementation of our proposed algorithm that fits CAR and SAR random effects.
Though the CAR/SAR models are widely used for spatial data analysis there are not many software packages which can be used for their model fitting.GLMMs with CAR/SAR random effects are often fitted in a Bayesian way using BUGS software (e.g., WinBUGS;Lunn et al., 2000, or alike), which leads to extremely slow computation due to its dependence on Markov chain Monte Carlo simulations.The R package INLA (Martins et al., 2013) provides a relatively fast Bayesian computation of spatial HGLMs via Laplace approximation of the posterior distribution where the model development has focused on continuous domain spatial modeling (Lindgren and Rue, 2015) whereas less focus has been on discrete models including CAR and SAR.Recently, package spaMM (Rousset and Ferdy, 2014) was developed to fit spatial HGLMs but is rather slow even for moderately sized data.Here, we extend package hglm to also provide fast computation of HGLMs with CAR and SAR random effects and at the same time an attempt has been made to improve the accuracy of the estimates by including corrections for the fixed effects.
The rest of the paper is organized as follows.First we give an introduction to CAR and SAR structures and present the h-likelihood theory together with the eigen decomposition of the covariance matrix of the CAR random effects and show how it simplifies the computation of the model.Then we present the R code implementation, show the use of the implementation for two real data examples and evaluate the package using simulations.The last section concludes the article.

CAR and SAR structures in HGLMs
HGLMs with spatially correlated random effects are commonly used in spatial data analysis (Cressie, 1993;Wall, 2004).A spatial HGLM with Gaussian CAR random effects is given by (1) with z s |u s ⊥ z t |u t , ∀s = t and where s represents a location identified by the coordinates (x (s) , y (s)), β is a vector of fixed effects, X s is the fixed-effect design matrix for location s, u s being the location specific random effects and Z s is a design matrix associated with u s .While gives a Gaussian SAR random effects structure.The D matrix in Equations ( 3) and ( 4) is some known function of the location coordinates (see, e.g., Clayton and Kaldor, 1987) and ρ is often referred to as the spatial dependence parameter (Hodges, 2013).In application, D is often a neighborhood matrix whose diagonal elements are all 0, and off-diagonal elements (s, t) are 1 if locations s and t are neighbors.If Σ = τI, i.e., there is no spatial correlation, then this model may be estimated by using usual software packages for GLMMs.With Σ = τD, i.e., defining the spatial correlation directly via the D matrix, the model can be fitted using earlier versions of the hglm package.However, for the general structure in Equations ( 3) and (4), estimation using explicit maximization of the marginal or profile likelihood involves quite advanced derivations as a consequence of partial differentiation of the likelihood including Σ with ρ and τ as parameters (see, e.g., Lee and Lee, 2012).In this paper, we show that by using eigen decomposition of Σ, we can modify an already existing R package, e.g., hglm, with minor programming effort, to fit a HGLM with CAR or SAR random effects.

h-likelihood estimation
In order to explain the specific model fitting algorithm implemented in hglm, we present a brief overview of the EQL algorithm for HGLMs (Lee and Nelder, 1996).First, we start with the standard HGLM containing only uncorrelated random effects.Then, we extend the discussion for CAR and SAR random effects.Though it is possible to have more than one random-effect term, in a HGLM, coming from different distributions among the conjugate distributions to GLM families, for presentational simplicity we consider only one random-effect in this section.The (log-)h-likelihood for a HGLM with independent random effects can be presented as where, i = 1, 2, . . ., n, j = 1, 2, . . .k, θ is the canonical parameter of the mean model, φ is the dispersion parameter of th mean model, b (.) a function which satisfies E y i,j |u i = µ i,j = b θ i,j and Var y i,j |u i = φb θ i,j = φV µ i,j , where V(.) is the GLM variance function.Furthermore, φ is the dispersion parameter of the mean model, θ R (u i ) = v i with θ R (.) being a so-called weak canonical link for u i (leading to conjugate HGLMs) or an identity link for Gaussian random effects leading to GLMMs (Lee et al., 2006, pp. 3-4), ψ = E (u i ), τ is the dispersion parameter, and b R and c R are some known functions depending on the distribution of u i .Equation ( 5) also implies that h can be defined uniquely by using the means and the mean-variance relations of y i,j |u i and the quasi-response ψ, allowing us representing it as a sum of two extended quasi-likelihoods (double EQL; Lee and Nelder, 2003) as The R Journal Vol.7/2, December 2015 ISSN 2073-4859 where d 0,i,j and d 1,i are the deviance components of y i,j |u i and ψ respectively and Var (ψ) = τV 1 (u i ).
It is worth noting that −2D is an approximation to h if y|u or u (or both) belongs/belong to Binomial, Poisson, double-Poisson (or quasi-Poisson), Beta, or Gamma families.However, no such equivalent relation (even approximately) is known for the quasi-Binomial family.
In order to estimate the model parameters, Lee and Nelder (1996) suggested a two-step procedure in which the first, for fixed dispersion parameters φ and τ, h (or equivalently −2D) is maximized w.r.t.β and v = {v i }.This maximization leads to an iteratively weighted least-square (IWLS) algorithm which solves where and y 1,a = are the vectors of GLM working responses for y i,j |u i and v i respectively, and = IW −1 with I being a diagonal matrix whose first n diagonal elements are all equal to φ and the remaining k diagonal elements are all equal to τ and W = diag (w 0 , w 1 ) where w 0 = and w 1 = In the second step, the following profile likelihood is maximized to estimate the dispersion parameters.
where β and v are obtained from the fist step.One needs to iterate between the two steps until convergence.This procedure is often referred to as "HL(0,1)" (Lee and Lee, 2012) and is available in the R packages spaMM and HGLMMM (Molas and Lesaffre, 2011).Because there is no unified algorithm to maximize p β,v computer packages often carry out the maximization by using general purpose optimization routines, e.g., package spaMM uses the optim function.
A unified algorithm for estimating the variance components can be derived by using profile likelihood adjustment in Equation ( 6) instead of ( 5).This leads to maximizing for estimating τ and φ.The corresponding score equation of the dispersion parameters, after ignoring the fact that β and v are functions of φ and τ, can be shown (see, e.g., Lee and Nelder, 2001) to have the form of Gamma family GLMs with d 0,i,j / 1 − h i,j and d 1,i / (1 − h i ) as the responses for φ and τ respectively, where h i,j and h i are the hat values corresponding to y i,j and v i in the first step, and 1 − h i,j /2 and (1 − h i ) /2 as the respective prior weights.This procedure is often referred to as EQL (Lee et al., 2006) or DEQL (Lee and Nelder, 2003) and was the only available procedure in hglm ( 1.2-8).An advantage of this algorithm is that fixed effects can also be fitted in the dispersion parameters (Lee et al., 2006) without requiring any major change in the algorithm.
HL(0,1) and EQL were found to be biased, especially for binary responses when the cluster size is small (Lee and Nelder, 2001;Noh and Lee, 2007) and when τ is large in both Binomial and Poisson GLMMs.Several adjustments are suggested and explained in the literature (see, e.g., Lee andNelder, 2001, 2003;Lee and Lee, 2012) to improve the performance of h-likelihood estimation.Among these alternative suggestions HL(1,1) is the most easily implementable and found to be computationally faster than the other alternatives (Lee and Lee, 2012).The HL(1,1) estimates the fixed effects, β, by maximizing p v (h) instead of h and can be implemented by adjusting the working response y a in the IWLS step, according to Lee and Lee (2012).In hglm ( 2.0), this correction for the β estimates has been implemented, through the method = "EQL1" option, but still p β,v (Q) is maximized for the estimation of the dispersion parameters, in order to make use of the unified algorithm for dispersion parameter estimates via Gamma GLMs.Table 1 shows available options for h-likelihood estimation (subject to the corrections mentioned above) in different R packages.
The above EQL method cannot be directly applied to HGLMs with CAR/SAR random effects because the resulting p β,v (Q) does not allow us to use the Gamma GLM to estimate the dispersion parameters, τ and ρ.In the following subsection we present a simplification which allows us to use the Gamma GLM for estimating τ and ρ.

Simplification of model computation
The main difference between an ordinary GLMM with independent random effects and CAR/SAR random effects models is that the random effects in the later cases are not independent.In the following we show that by using the eigen decomposition we can reformulate a GLMM with CAR or SAR random effects to an equivalent GLMM with independent but heteroscedastic random effects.
Lemma 1 Let ω = {ω i } n i=1 be the eigenvalues of D, and V is the matrix whose columns are the corresponding orthonormal eigenvectors, then 1 where Λ is diagonal matrix whose ith diagonal element is given by It is worth noting that the relation between the eigenvalues of D and Σ was already known as early as in Ord (1975) and was used to simplify the likelihood function of the Gaussian CAR model.Similarly, for simplification of the SAR model, we have 1 An anonymous referee has pointed out that a somewhat extended version of Lemma 1, which deals with simultaneous diagonalization of two positive semi-definite matrices (Newcomb, 1969), has already been used to simplify the computation of Gaussian response models with intrinsic CAR random effects, especially in a Bayesian context (see, e.g., He et al., 2007, and the related discussions in Hodges 2013, Chap.5).Therefore, Lemma 1 is not an original contribution of this paper, however, to the authors' knowledge, no software package has yet utilized this convenient relationship to fit any HGLM by using interconnected GLMs, as discussed below.
Re-arranging Equation (1), we have where, η = {η s }, X = {X s }, Z = {Z s } V and u * = V T u.With the virtue of Lemma 1 and the properties of the multivariate normal distribution, we see that u * ∼ N 0, Λ −1 .This model can now easily be fitted using the hglm package in R. Further note that, for the CAR model, from Lemma 1, we have The R Journal Vol.7/2, December 2015 ISSN 2073-4859

Option Explanation
rand.family = CAR(D = nbr) The random effect has conditional or simulrand.family= SAR(D = nbr) taneous autoregressive covariance structure.
Here nbr is a matrix provided by the user.method = "EQL1" The first order correction of fixed effects (Lee and Lee, 2012) applied on the EQL estimates.rand.disp.X = X Linear predictor for the variance component of a random effect.The matrix X is provided by the user.

rand.family =
This option provides the possibility of having list(Gamma(),CAR()) different distributions for the random effects.where θ 0 = 1/τ and θ 1 = −ρ/τ.While for the SAR model, where θ 0 = 1/ √ τ and θ 1 = −ρ/ √ τ.Following Lee et al. (2006), we can use an inverse link for CAR or an inverse square-root link for SAR in a Gamma GLM with response u i * 2 / (1 − h i ), where h i is the corresponding hat value from the mean model (Equation ( 13)), (1 − h i ) /2 as the weight and the linear predictor given by Equations ( 14) and ( 15) to obtain an EQL estimate of θ 0 and θ 1 .

Implementation
The spatial models above are implemented in the hglm package ( 2.0) by defining new families, CAR() and SAR(), for the random effects.The "CAR" and "SAR" families allow the user to define spatial structures by specifying the D matrix.Using Lemma 1, these families are created using D with default link function "identity" for the Gaussian random effects u, and "inverse" and "inverse.sqrt"for the dispersion models ( 14) and (15).When the "CAR" or "SAR" family is specified for the random effects, the parameter estimates ρ and τ are given in the output in addition to the other summary statistics of the dispersion models.The new input options in hglm are described in Table 2.
Throughout the examples below, we use Gaussian CAR random effects to introduce spatial correlation.However, one can also add additional independent non-Gaussian random effects along with the Gaussian CAR random effect.For example, an overdispersed count response can be fitted using a Poisson HGLM with an independent Gamma and Gaussian CAR random effects.

Examples and simulation study
In the hglm package vignette, we look at the improvement of EQL1 in comparison to EQL.Here, we focus on the precision of the parameter estimates with spatial HGLMs.

Poisson CAR & SAR model
We study the properties of the estimates produced by hglm using a simulation study built around the Scottish Lip Cancer example (see also the examples).We simulate data with the same X values, offset and neighborhood matrix as in the Scotthis Lip Cancer example data.We use the true values of the parameters, after Lee and Lee (2012), as intercept, β f pp , τ, ρ = (0.25, 0.35, 1.5, 0.1).The parameter estimates for 1000 Monte Carlo iterations are summarized in Table 3.
Both the EQL and the EQL1 correction are slightly biased for τ and ρ (Table 3) though the absolute amount of bias is small and may be negligible in practical applications.The EQL1 correction mainly improves the estimates of the intercept term.There were convergence problems for a small number of replicates, which was not surprising given the small number of observations (n = 56) and that the simulated value for the spatial autocorrelation parameter ρ connects the effects in the different Scottish districts rather weakly.Such convergence problems can be addressed by pre-specifying better starting values.For the converged estimates the bias was small and negligible.The simulation results also revealed that the distribution of ρ from CAR models is skewed (see Figure 1), which was also pointed out by Lee and Lee (2012).However, the distribution of − ρ/ τ turned out to be less skewed than the distribution of ρ.Similar observation was also found for SAR models.This suggests, we might draw any inference on spatial variance-covariance parameters in the transformed scale, θ 0 and θ 1 .

Computational efficiency
Fitting CAR and SAR models for large data sets could be computationally challenging, especially for the spatial variance-covariance parameters.Using our new algorithm in hglm, moderately sized problems can be fitted efficiently.
We re-sampled from the ohio data set (for more details on the data set see also the examples) for different number of locations, each with 10 replicates, executed on a single Intel© Xeon© E5520 2.27GHz CPU.In each replicate, the data was fitted using both hglm and spaMM, for the same model described above.The results regarding average computational time are summarized in Figure 2. hglm is clearly more usable for fitting larger sized data sets.In the package vignette, we also show the comparisons of the parameter estimates, where the "EQL" to the "EQL-" estimates from spaMM, with high correlation coefficients between the two methods: 1.0000, 0.9998, 0.9997 and 0.9999 for the residual variance, ρ, τ and the intercept, respectively.

Examples Scottish Lip Cancer data set
Here we analyze the cancer data set (source: Clayton and Kaldor, 1987)    For the CAR model, the hglm estimates are exactly the same as those labelled as "PQL" estimates in Lee and Lee (2012).To get the EQL1 correction the user has to add the option method = "EQL1" and the hglm function gives similar results to those reported in Lee and Lee (2012), e.g., their HL(1,1) estimates were intercept, β f pp , τ, ρ = (0.238, 0.0376, 0.155, 0.174) whereas our "EQL1" correction gives (0.234, 0.0377, 0.156, 0.174) (see also Section "Fitting a spatial Markov Random Field model using the CAR family" in the package vignette).A minor difference to the "EQL1" result appears because Lee and Lee (2012) used the HL(1,1) modification to their HL(0,1) whereas we apply such a correction directly to EQL which is slightly different from HL(0,1) (see Table 1).

Ohio elementary school grades data set
We analyze a data set consisting of the student grades of 1,967 Ohio Elementary Schools during the year 2001-2002.The data set is freely available on the internet (URL http://www.spatial-econometrics.com/) as a web supplement to LeSage and Pace (2009) but was not analyzed therein.The shape files were downloaded from http://www.census.gov/cgi-bin/geo/shapefiles2013/mainand the districts of 1,860 schools in these two files could be connected unambiguously.The data set contains information on, for instance, school building ID, Zip code of the location of the school, proportion of passing on five subjects, number of teachers, number of students, etc.We regress the median of 4th grade proficiency scores, y, on an intercept, based on school districts.The statistical model is given as where i = 1, 2, . . ., 1860 (observations), j = 1, 2, . . ., 616 (districts), i,j ∼ N 0, σ 2 e , {v j } = v ∼ N 0, τ (I − ρW) −1 and W = {w p,q } 616 p,q=1 is a spatial weight matrix (i.e., the neighborhood matrix).We construct w p,q = 1 if the two districts p and q are adjacent, and w p,q = 0 otherwise.
The above choice of constructing the weight matrix is rather simple.Because the aim of this paper is to demonstrate the use of hglm for fitting spatial models rather than drawing conclusions from real data analysis, we skip any further discussion on the construction of the weight matrix.Interested, readers are referred to LeSage and Pace (2009) for more discussion on the construction of the spatial weight matrices.With the spatial weight matrix defined as above, we can estimate model ( 16) by using our hglm package in the following way.> districtShape <-as.numeric(substr(as.character(ohioShape@data$UNSDIDFP),3, 7)) > > ## calculate fitted values from the CAR model > CARfit <-matrix(ohioCAR$ranef + ohioCAR$fixef, + dimnames = list(rownames(ohioDistrictDistMat), NULL)) > ohioShape@data$CAR <-CARfit[as.character(districtShape),]> is.na(ohioShape@data$CAR[353]) <-TRUE # remove estimate of Lake Erie > > ## visualize the results > spplot(ohioShape, zcol = "CAR", main = "Fitted values from CAR", + col.regions = heat.colors(1000)[1000:1],cuts = 1000) A predict() method is not available because predicting spatially correlated random effects for autoregressive models requires re-fitting the whole model.Thus standard kriging cannot be used because the covariance structure changes if the neighborhood matrix is altered, while keeping τ and ρ The R Journal Vol.7/2, December 2015 ISSN 2073-4859 unchanged.Instead, the fitted model needs to include the entire neighborhood matrix with districts having missing data as well.Consequently, the incidence matrix Z has more columns than rows.This method of predicting random effects is frequently used in animal breeding applications (Henderson, 1984) but, to our knowledge, has not been applied to spatial autoregressive models previously.
In the example above, ohioMedian$district has 616 levels and 54 of the districts have no records (Figure 3(A)).The incidence matrix Z, created using the model.matrixfunction, therefore has 616 columns and 54 of these are columns of zeros.Hence, there are 616 levels in the fitted spatial random effect giving predictions for the districts without records.

Conclusion
The hglm package is one of few non-Bayesian packages on CRAN to fit spatial HGLMs, where the fixed and random effects are estimated simultaneously.We have shown how the HGLM framework, allowing linear predictors to model variance components, can be exploited to fit CAR and SAR models.This gives a computationally efficient algorithm for moderately sized problems (number of locations < approx.5000).
0102 * * Significantly different from 0 at the 5% level.† The estimates are the means from 1000 replicates.Table3: Average bias in parameter estimates in the simulation using the Scottish Lip Cancer example.

Figure 1 :
Figure 1: Density plot of the spatial variance-covariance parameter estimates from CAR models via "EQL1" over 1000 Monte Carlo simulations for the Scottish Lip Cancer example.

Figure 2 :
Figure 2: Comparison of the execution time for fitting CAR models using hglm and spaMM.

Figure 3 :
Figure 3: Observed (A) and predicted (B) median 4th grade proficiency scores of the school districts in Ohio.Districts without any observations are displayed in gray.

Table 1 :
Implemented h-likelihood methods in hglm compared to spaMM and HGLMMM packages.

Table 2 :
New options in the hglm function.
from the hglm package.Calling data(cancer) loads a numeric vector E that represents the expected number of skin cancer patients in different districts in Scotland, a numeric vector O giving the corresponding observed counts, a numeric vector Paff giving proportion of people involved agriculture, farming, and fisheries, and matrix D giving the neighborhood structure among Scottish districts.Here we demonstrate how the data is fitted as a CAR model or a SAR model, using the hglm package with the EQL method.