Fit spatial econometrics models with INLA

Share:

Description

These functions fit some spatial econometrics models for a given value of rho (the spatial autocorrelation parameter). sem.inla fits a spatial error model, slm fits a spatial lag model and sdm.inla fits a spatial Durbin model.

Usage

1
2
3
4
5
6
sem.inla(formula, d, W, rho, improve = TRUE, impacts = FALSE, fhyper = NULL, 
   probit = FALSE, ...)
slm.inla(formula, d, W, rho, mmatrix = NULL, improve = TRUE, impacts = FALSE, 
   fhyper = NULL, probit = FALSE, ...)
sdm.inla(formula, d, W, rho, mmatrix = NULL, intercept = TRUE, impacts = FALSE,
   improve = TRUE, fhyper = NULL, probit = FALSE, ...)

Arguments

formula

Formula with the response variable, the fixed effects and, possibly, other non-linear effects.

d

Data.frame with the data.

W

Adjacency matrix.

rho

Value of the spatial autocorrelation paramter.

mmatrix

Design matrix of fixed effects.

intercept

Logical. Whether an intercept has been included in the model.

improve

Logical. Whether improve model fitting (this may require more computing time).

impacts

Logical. Whether impacts are computed.

fhyper

Options to be passed to the definition of the hyper-parameters in the spatial effects.

probit

Logical. Whether a probit model is used. Note this is only used when computing the impacts and that argument family must be set accordingly.

...

Other arguments passed to function inla.

Details

These functions fit a spatial econometrics model with a fixed value of the spatial autocorrelation parameter rho.

In addition, the marginal -log-likelihood is corrected to account for the variance-covariance matrix of the error term or random effects.

Value

An inla object.

Author(s)

Virgilio G<f3>mez-Rubio <virgilio.gomez@uclm.es>

References

Roger S. Bivand, Virgilio G<f3>mez-Rubio, H<e5>vard Rue (2014). Approximate Bayesian inference for spatial econometrics models. Spatial Statistics, Volume 9, 146-165.

Roger S. Bivand, Virgilio G<f3>mez-Rubio, H<e5>vard Rue (2015). Spatial Data Analysis with R-INLA with Some Extensions. Journal of Statistical Software, 63(20), 1-31. URL http://www.jstatsoft.org/v63/i20/.

See Also

leroux.inla

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
## Not run: 

if(requireNamespace("INLA", quietly = TRUE)) {
require(INLA)
require(spdep)

data(columbus)

lw <- nb2listw(col.gal.nb, style="W")

#Maximum Likelihood (ML) estimation
colsemml <- errorsarlm(CRIME ~ INC + HOVAL, data=columbus, lw, method="eigen", 
	quiet=FALSE)
colslmml <- lagsarlm(CRIME ~ INC + HOVAL, data=columbus, lw, method="eigen", 
	type="lag", quiet=FALSE)
colsdmml <- lagsarlm(CRIME ~ INC + HOVAL, data=columbus, lw, method="eigen", 
	type="mixed", quiet=FALSE)

#Define grid on rho
rrho<-seq(-1, .95, length.out=40)

#Adjacency matrix
W <- as(as_dgRMatrix_listw(nb2listw(col.gal.nb)), "CsparseMatrix")
#Index for spatial random effects
columbus$idx<-1:nrow(columbus)

#Formula
form<- CRIME ~ INC + HOVAL

zero.variance = list(prec=list(initial = 25, fixed=TRUE))



seminla<-mclapply(rrho, function(rho){

                sem.inla(form, d=columbus, W=W, rho=rho,
                        family = "gaussian", impacts=FALSE,
                        control.family = list(hyper = zero.variance),
                        control.predictor=list(compute=TRUE),
                        control.compute=list(dic=TRUE, cpo=TRUE),
                        control.inla=list(print.joint.hyper=TRUE), 
				#tolerance=1e-20, h=1e-6),
			verbose=FALSE
                )

})



slminla<-mclapply(rrho, function(rho){

                slm.inla(form, d=columbus, W=W, rho=rho,
                        family = "gaussian", impacts=FALSE,
                        control.family = list(hyper = zero.variance),
                        control.predictor=list(compute=TRUE),
                        control.compute=list(dic=TRUE, cpo=TRUE),
                        control.inla=list(print.joint.hyper=TRUE), 
				#tolerance=1e-20, h=1e-6),
			verbose=FALSE
                )
})


sdminla<-mclapply(rrho, function(rho){

                sdm.inla(form, d=columbus, W=W, rho=rho,
                        family = "gaussian", impacts=FALSE,
                        control.family = list(hyper = zero.variance),
                        control.predictor=list(compute=TRUE),
                        control.compute=list(dic=TRUE, cpo=TRUE),
                        control.inla=list(print.joint.hyper=TRUE), 
				#tolerance=1e-20, h=1e-6),
			verbose=FALSE
                )
})

#BMA using a uniform prior (in the log-scale) and using a Gaussian 
#approximation to the marginal
sembma<-INLABMA(seminla, rrho, 0, usenormal=TRUE)
slmbma<-INLABMA(slminla, rrho, 0, usenormal=TRUE)
sdmbma<-INLABMA(sdminla, rrho, 0, usenormal=TRUE)

#Display results
plot(sembma$rho$marginal, type="l", ylim=c(0,5))
lines(slmbma$rho$marginal, lty=2)
lines(sdmbma$rho$marginal, lty=3)
#Add ML estimates
abline(v=colsemml$lambda, col="red")
abline(v=colslmml$rho, col="red", lty=2)
abline(v=colsdmml$rho, col="red", lty=3)
#Legend
legend(-1,5, c("SEM", "SLM", "SDM"), lty=1:3)
}


## End(Not run)