diff --git a/docs/404.html b/docs/404.html index 366c8a6..7e41807 100644 --- a/docs/404.html +++ b/docs/404.html @@ -32,7 +32,7 @@
@@ -106,7 +106,7 @@vignettes/SpatialFiltering.Rmd
@@ -96,13 +97,36 @@ The Moran eigenvector approach (Dray, Legendre, and Peres-Neto 2006; Griffith and Peres-Neto 2006) involved the spatial patterns represented by maps of eigenvectors; by choosing suitable orthogonal patterns and adding them to a linear or generalised linear model, the spatial dependence present in the residuals can be moved into the model.
-It uses brute force to search the set of eigenvectors of the matrix \(\mathbf{M W M}\), where
-\[\mathbf{M} = \mathbf{I} - \mathbf{X}(\mathbf{X}^{\rm T}
-\mathbf{X})^{-1}\mathbf{X}^{\rm T}\] is a symmetric and idempotent projection matrix and \(\mathbf{W}\) are the spatial weights. In the spatial lag form of SpatialFiltering
and in the GLM ME
form below, \(\mathbf{X}\) is an \(n\)-vector of ones, that is the intercept only.
In its general form, SpatialFiltering
chooses the subset of the \(n\) eigenvectors that reduce the residual spatial autocorrelation in the error of the model with covariates. The lag form adds the covariates in assessment of which eigenvectors to choose, but does not use them in constructing the eigenvectors. SpatialFiltering
was implemented and contributed by Yongwan Chun and Michael Tiefelsdorf, and is presented in Tiefelsdorf and Griffith (2007); ME
is based on Matlab code by Pedro Peres-Neto and is discussed in Dray, Legendre, and Peres-Neto (2006) and Griffith and Peres-Neto (2006).
The Moran eigenvector approach (Dray, Legendre, and Peres-Neto +2006; Griffith and Peres-Neto +2006) involved the spatial patterns represented by maps of +eigenvectors; by choosing suitable orthogonal patterns and adding them +to a linear or generalised linear model, the spatial dependence present +in the residuals can be moved into the model.
+It uses brute force to search the set of eigenvectors of the matrix +\(\mathbf{M W M}\), where
+\[\mathbf{M} = \mathbf{I} -
+\mathbf{X}(\mathbf{X}^{\rm T}
+\mathbf{X})^{-1}\mathbf{X}^{\rm T}\] is a symmetric and
+idempotent projection matrix and \(\mathbf{W}\) are the spatial weights. In
+the spatial lag form of SpatialFiltering
and
+in the GLM ME
form below, \(\mathbf{X}\) is an \(n\)-vector of ones, that is the intercept
+only.
In its general form, SpatialFiltering
+chooses the subset of the \(n\)
+eigenvectors that reduce the residual spatial autocorrelation in the
+error of the model with covariates. The lag form adds the covariates in
+assessment of which eigenvectors to choose, but does not use them in
+constructing the eigenvectors.
+SpatialFiltering
was implemented and
+contributed by Yongwan Chun and Michael Tiefelsdorf, and is presented in
+Tiefelsdorf and Griffith (2007);
+ME
is based on Matlab code by Pedro Peres-Neto
+and is discussed in Dray, Legendre, and
+Peres-Neto (2006)
+and Griffith and Peres-Neto (2006).
library(spdep)
require("sf", quietly=TRUE)
@@ -144,7 +168,7 @@ Introduction## fitted(nySFE)vec14 -1.05105 0.60534 -1.736 0.083662 .
## fitted(nySFE)vec75 1.90600 0.60534 3.149 0.001826 **
## fitted(nySFE)vec21 -1.06331 0.60534 -1.757 0.080138 .
-## fitted(nySFE)vec36 -1.17861 0.60534 -1.947 0.052578 .
+## fitted(nySFE)vec36 1.17861 0.60534 1.947 0.052578 .
## fitted(nySFE)vec61 -1.08582 0.60534 -1.794 0.073986 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
@@ -164,8 +188,23 @@ Introduction## 2 267 97.837 10 21.782 5.9444 3.988e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
-Since the SpatialFiltering
approach does not allow weights to be used, we see that the residual autocorrelation of the original linear model is absorbed, or ‘whitened’ by the inclusion of selected eigenvectors in the model, but that the covariate coefficients change little. The addition of these eigenvectors – each representing an independent spatial pattern – relieves the residual autocorrelation, but otherwise makes few changes in the substantive coefficient values.
The ME
function also searches for eigenvectors from the spatial lag variant of the underlying model, but in a GLM framework. The criterion is a permutation bootstrap test on Moran’s \(I\) for regression residuals, and in this case, because of the very limited remaining spatial autocorrelation, is set at \(\alpha = 0.5\). Even with this very generous stopping rule, only few eigenvectors are chosen; their combined contribution only just improves the fit of the GLM model.
Since the SpatialFiltering
approach does
+not allow weights to be used, we see that the residual autocorrelation
+of the original linear model is absorbed, or ‘whitened’ by the inclusion
+of selected eigenvectors in the model, but that the covariate
+coefficients change little. The addition of these eigenvectors – each
+representing an independent spatial pattern – relieves the residual
+autocorrelation, but otherwise makes few changes in the substantive
+coefficient values.
The ME
function also searches for
+eigenvectors from the spatial lag variant of the underlying model, but
+in a GLM framework. The criterion is a permutation bootstrap test on
+Moran’s \(I\) for regression residuals,
+and in this case, because of the very limited remaining spatial
+autocorrelation, is set at \(\alpha =
+0.5\). Even with this very generous stopping rule, only few
+eigenvectors are chosen; their combined contribution only just improves
+the fit of the GLM model.
Figure \[fig:geigen2\] shows the spatial patterns chosen to match the very small amount of spatial autocorrelation remaining in the model. As with the other Poisson regressions, the closeness to TCE sites is highly significant. Since, however, many TCE sites are also in or close to more densely populated urban areas with the possible presence of both point-source and non-point-source pollution, it would be premature to take such results simply at their face value. There is, however, a potentially useful contrast between the cities of Binghamton in the south of the study area with several sites in its vicinity, and Syracuse in the north without TCE sites in this data set.
+Figure \[fig:geigen2\] shows the +spatial patterns chosen to match the very small amount of spatial +autocorrelation remaining in the model. As with the other Poisson +regressions, the closeness to TCE sites is highly significant. Since, +however, many TCE sites are also in or close to more densely populated +urban areas with the possible presence of both point-source and +non-point-source pollution, it would be premature to take such results +simply at their face value. There is, however, a potentially useful +contrast between the cities of Binghamton in the south of the study area +with several sites in its vicinity, and Syracuse in the north without +TCE sites in this data set.
This vignette formed pp. 302–305 of the first edition of Bivand, R. S., Pebesma, E. and Gómez-Rubio V. (2008) Applied Spatial Data Analysis with R, Springer-Verlag, New York. It was retired from the second edition (2013) to accommodate material on other topics, and is made available in this form with the understanding of the publishers.↩︎
This vignette formed pp. 302–305 of the first edition of +Bivand, R. S., Pebesma, E. and Gómez-Rubio V. (2008) Applied Spatial +Data Analysis with R, Springer-Verlag, New York. It was retired from the +second edition (2013) to accommodate material on other topics, and is +made available in this form with the understanding of the publishers.↩︎
Site built with pkgdown 2.0.6.
+Site built with pkgdown 2.0.7.
vignettes/nb_igraph.Rmd
@@ -98,27 +100,77 @@ Since the spdep package was created, spatial weights objects have been constructed as lists with three components and a few attributes, in old-style class listw
objects. The first component of a listw
object is an nb
object, a list of n
integer vectors, with at least a character vector region.id
attribute with n
unique values (like the row.names
of a data.frame
object); n
is the number of spatial entities. Component i
of this list contains the integer identifiers of the neighbours of i
as a sorted vector with no duplication and values in 1:n
; if i
has no neighbours, the component is a vector of length 1
with value 0L
. The nb
object may contain an attribute indicating whether it is symmetric or not, that is whether i
is a neighbour of j
implies that j
is a neighbour of i
. Some neighbour definitions are symmetric by construction, such as contiguities or distance thresholds, others are asymmetric, such as k
-nearest neighbours. The nb
object redundantly stores both i
-j
and j
-i
links.
The second component of a listw
object is a list of n
numeric vectors, each of the same length as the corresponding non-zero vectors in the nb
object. These give the values of the spatial weights for each i
-j
neighbour pair. It is often the case that while the neighbours are symmetric by construction, the weights are not, as for example when weights are row-standardised by dividing each row of input weights by the count of neighbours or cardinality of the neighbour set of i
. In the nb2listw
function, it is also possible to pass through general weights, such as inverse distances, shares of boundary lengths and so on.
The third component of a listw
object records the style
of the weights as a character code, with "B"
for binary weights taking values zero or one (only one is recorded), "W"
for row-standardised weights, and so on. In order to subset listw
objects, knowledge of the style
may be necessary
It is obvious that this is similar to the way in which sparse matrices are stored, either by row - like the listw
object, or by column. The key insight is that storing zero values is unnecessary, as we only need to store the row and column locations of non-zero values. Early on, a Netlib library was used to provide limited support in spdep for sparse matrices, followed by functionality in SparseM, spam, and Matrix.
From spdep and spatialreg versions 1.2, this vignette and accompanying functionality has been moved to spatialreg.
+Since the spdep package was created, spatial
+weights objects have been constructed as lists with three
+components and a few attributes, in old-style class listw
+objects. The first component of a listw
object is an
+nb
object, a list of n
integer vectors, with
+at least a character vector region.id
attribute with
+n
unique values (like the row.names
of a
+data.frame
object); n
is the number of spatial
+entities. Component i
of this list contains the integer
+identifiers of the neighbours of i
as a sorted vector with
+no duplication and values in 1:n
; if i
has no
+neighbours, the component is a vector of length 1
with
+value 0L
. The nb
object may contain an
+attribute indicating whether it is symmetric or not, that is whether
+i
is a neighbour of j
implies that
+j
is a neighbour of i
. Some neighbour
+definitions are symmetric by construction, such as contiguities or
+distance thresholds, others are asymmetric, such as
+k
-nearest neighbours. The nb
object
+redundantly stores both i
-j
and
+j
-i
links.
The second component of a listw
object is a list of
+n
numeric vectors, each of the same length as the
+corresponding non-zero vectors in the nb
object. These give
+the values of the spatial weights for each i
-j
+neighbour pair. It is often the case that while the neighbours are
+symmetric by construction, the weights are not, as for example when
+weights are row-standardised by dividing each row of input
+weights by the count of neighbours or cardinality of the neighbour set
+of i
. In the nb2listw
function, it is also
+possible to pass through general weights, such as inverse distances,
+shares of boundary lengths and so on.
The third component of a listw
object records the
+style
of the weights as a character code, with
+"B"
for binary weights taking values zero or one (only one
+is recorded), "W"
for row-standardised weights, and so on.
+In order to subset listw
objects, knowledge of the
+style
may be necessary
It is obvious that this is similar to the way in which sparse
+matrices are stored, either by row - like the listw
object,
+or by column. The key insight is that storing zero values is
+unnecessary, as we only need to store the row and column locations of
+non-zero values. Early on, a Netlib library was used to provide limited
+support in spdep for sparse matrices, followed by
+functionality in SparseM, spam, and
+Matrix.
From spdep and spatialreg versions +1.2, this vignette and accompanying functionality has been moved to +spatialreg.
Since Matrix is a recommended package, its functionality has increasingly been used over time, and it has become one of two packages on which spatialreg depends. This is reported on loading:
+Since Matrix is a recommended package, its +functionality has increasingly been used over time, and it has become +one of two packages on which spatialreg depends. This +is reported on loading:
## Loading required package: spData
## Loading required package: Matrix
## Loading required package: sf
-## Linking to GEOS 3.11.0, GDAL 3.5.2, PROJ 9.0.1; sf_use_s2() is TRUE
+## Linking to GEOS 3.12.1, GDAL 3.8.0, PROJ 9.3.0; sf_use_s2() is TRUE
The legacy Columbus OH data set has 49 spatial entities, polygons, defined as the boundaries of policing districts in the city. spatialreg imports from spdep which in turn depends on sf.
+The legacy Columbus OH data set has 49 spatial entities, polygons, +defined as the boundaries of policing districts in the city. +spatialreg imports from spdep which in +turn depends on sf.
dothis <- TRUE
if (!suppressPackageStartupMessages(require(sf, quietly=TRUE))) {
@@ -129,9 +181,9 @@ Getting some datasf_extSoftVersion()
}
## GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H
-## "3.11.0" "3.5.2" "9.0.1" "true" "true"
+## "3.12.1" "3.8.0" "9.3.0" "true" "true"
## PROJ
-## "9.0.1"
+## "9.3.0"
library(sf)
columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1])
Contiguous neighbours are often used for polygonal spatial entities, here with the poly2nb function defaulting to the queen criterion - entities are neighbours if they share a boundary point. We see that the entity IDs are copied across to the nb
object:
Contiguous neighbours are often used for polygonal spatial entities,
+here with the poly2nb function defaulting to the
+queen criterion - entities are neighbours if they share a
+boundary point. We see that the entity IDs are copied across to the
+nb
object:
nb_q <- spdep::poly2nb(columbus)
nb_q
spdep::is.symmetric.nb(nb_q)
## [1] TRUE
-In order to make the object more complicated, let us drop the neighbour links for the 21st entity (noting that the print method reports the ID of the entity with no neighbours, not its number in 1:n
), and plot the resulting map of neighbours:
In order to make the object more complicated, let us drop the
+neighbour links for the 21st entity (noting that the print method
+reports the ID of the entity with no neighbours, not its number in
+1:n
), and plot the resulting map of neighbours:
col2 <- spdep::droplinks(nb_q, 21)
nb_q[[21]]
spdep::is.symmetric.nb(col2)
## [1] TRUE
@@ -194,7 +254,15 @@ At present only listw
objects can be coerced to objects of classes defined in Matrix. Because the style
is lost on coercion, it may not be possible to reconstruct spatial weights as the sparse matrix representation does not preserve it. We will start with symmetric binary weights, first creating a spatial weights object, and signalling that one entity has no neighbours with the zero.policy
argument (default false). The matrix and graph representations of no-neighbour entities are not obvious.
At present only listw
objects can be coerced to objects
+of classes defined in Matrix. Because the
+style
is lost on coercion, it may not be possible to
+reconstruct spatial weights as the sparse matrix representation does not
+preserve it. We will start with symmetric binary weights, first creating
+a spatial weights object, and signalling that one entity has no
+neighbours with the zero.policy
argument (default false).
+The matrix and graph representations of no-neighbour entities are not
+obvious.
nb_B <- spdep::nb2listw(col2, style="B", zero.policy=TRUE)
nb_B$style
spdep provides coercion methods from listw
to the "symmetricMatrix"
, "RsparseMatrix"
and "CsparseMatrix"
classes defined in Matrix. The "RsparseMatrix"
is the representation that is most similar to listw
, as it is row-based, but it is used less frequently in operations on sparse matrices. The entity IDs are passed using sparse matrix row and column names at present. Here we believe that our listw
object can be represented as a symmetric matrix, storing only a triangle rather than both i
-j
and j
-i
weights. The coercion method does check whether symmetry is present before proceeding:
spdep provides coercion methods from
+listw
to the "symmetricMatrix"
,
+"RsparseMatrix"
and "CsparseMatrix"
classes
+defined in Matrix. The "RsparseMatrix"
is
+the representation that is most similar to listw
, as it is
+row-based, but it is used less frequently in operations on sparse
+matrices. The entity IDs are passed using sparse matrix row and column
+names at present. Here we believe that our listw
object can
+be represented as a symmetric matrix, storing only a triangle rather
+than both i
-j
and
+j
-i
weights. The coercion method does check
+whether symmetry is present before proceeding:
## [1] TRUE
@@ -223,53 +302,92 @@
rownames(B)[1:10]
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
-Let us now try to retreive the list of neighbours from the symmetric sparse matrix. At present, we have to coerce from one Matrix internal representation to another in order to get to the "dgCMatrix"
format used inside mat2listw
, so we coerce to "dgTMatrix"
from "dsTMatrix"
. The style is not retreived automatically, but is set to "M"
to indicate conversion from a matrix. The neighbour links are retreived correctly, as are the IDs:
Let us now try to retreive the list of neighbours from the symmetric
+sparse matrix. At present, we have to coerce from one
+Matrix internal representation to another in order to
+get to the "dgCMatrix"
format used inside
+mat2listw
, so we coerce to "dgTMatrix"
from
+"dsTMatrix"
. The style is not retreived automatically, but
+is set to "M"
to indicate conversion from a matrix. The
+neighbour links are retreived correctly, as are the IDs:
## Warning in sn2listw(df): 21 is not an origin
-++## Warning in spdep::mat2listw(as(as(B, "TsparseMatrix"), "CsparseMatrix")): style +## is M (missing); style should be set to a valid value
+## Warning in sn2listw(df, style = style, zero.policy = zero.policy, +## from_mat2listw = TRUE): style is M (missing); style should be set to a valid +## value
+## Warning in sn2listw(df, style = style, zero.policy = zero.policy, +## from_mat2listw = TRUE): 21 is not an origin
+## Warning in sn2listw(df, style = style, zero.policy = zero.policy, from_mat2listw = TRUE): no-neighbour observations found, set zero.policy to TRUE; +## this warning will soon become an error
+## Warning in spdep::mat2listw(as(as(B, "TsparseMatrix"), "CsparseMatrix")): no-neighbour observations found, set zero.policy to TRUE; +## this warning will soon become an error
nb_B1$style
-## [1] "M"
+-all.equal(nb_B1$neighbours, col2, check.attributes=FALSE)
+all.equal(nb_B1$neighbours, col2, check.attributes=FALSE)
- +## [1] TRUE
## [1] TRUE
An initial reason for implementing support for sparse weights matrices in spdep was to permit the calculation of the log determinant term in spatial regressions for larger data sets. Using the eigenvalue approach with for example spatialreg::eigenw
is limited by the need to operate on dense matrices in memory to solve the eigenproblem:
+An initial reason for implementing support for sparse weights +matrices in spdep was to permit the calculation of the +log determinant term in spatial regressions for larger data sets. Using +the eigenvalue approach with for example
+spatialreg::eigenw
+is limited by the need to operate on dense matrices in memory to solve +the eigenproblem:rho <- 0.1 do_spatialreg <- FALSE if (requireNamespace("spatialreg", quietly=TRUE)) do_spatialreg <- TRUE if (do_spatialreg) sum(log(1 - rho * spatialreg::eigenw(nb_B)))
-## [1] -1.44787
When
-n
is large, this may become impractical and/or time-consuming, but does permit the rapid calculation of values of the log determinant for differing values of the spatial coefficient \(\rho\). The Matrix package provides manydeterminant
methods, here for a"dsCMatrix"
resulting from subtracting a"dsCMatrix"
, the product of a scalar and a"dsTMatrix"
, from a"ddiMatrix"
. The value of the log determinant follows, calling a sparse Cholesky decomposition internally for suitable input matrices.+When
+n
is large, this may become impractical and/or +time-consuming, but does permit the rapid calculation of values of the +log determinant for differing values of the spatial coefficient \(\rho\). The Matrix package +provides manydeterminant
methods, here for a +"dsCMatrix"
resulting from subtracting a +"dsCMatrix"
, the product of a scalar and a +"dsTMatrix"
, from a"ddiMatrix"
. The value of +the log determinant follows, calling a sparse Cholesky decomposition +internally for suitable input matrices.-## [1] "dgCMatrix" ## attr(,"package") ## [1] "Matrix"
@@ -977,7 +418,7 @@+c(determinant(I - rho * B, logarithm=TRUE)$modulus)
-## [1] -1.44787
The computation of a sparse Cholesky decomposition for each value of the spatial coefficient \(\rho\) may be avoided by updating a pre-computed object; this approach provides fast and accurate log determinants for larger (but not very large) data sets:
-+The computation of a sparse Cholesky decomposition for each value of +the spatial coefficient \(\rho\) may be +avoided by updating a pre-computed object; this approach provides fast +and accurate log determinants for larger (but not very large) data +sets:
++nW <- -B nChol <- Cholesky(nW, Imult=8) n * log(rho) + (2 * c(determinant(update(nChol, nW, 1/rho))$modulus))
## Warning in .local(x, logarithm, ...): the default value of argument 'sqrt' of +## method 'determinant(<CHMfactor>, <logical>)' may change from TRUE to FALSE as +## soon as the next release of Matrix; set 'sqrt' when programming
## [1] 15.40218
Asymmetric sparse matrices
-The use of row-standardisation leads to asymmetry even if the underlying neighbours are symmetric, unless all entities have matching numbers of neighbours (for example a regular grid on a torus):
-@@ -290,9 +290,6 @@+-The use of row-standardisation leads to asymmetry even if the +underlying neighbours are symmetric, unless all entities have matching +numbers of neighbours (for example a regular grid on a torus):
+ @@ -282,42 +400,53 @@Asymmetric sparse matrices## .. ..$ : chr [1:49] "1" "2" "3" "4" ... ## ..@ x : num [1:230] 0.333 0.25 0.5 0.25 0.25 ... ## ..@ factors : list()
diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png index 5fe2fe3..881792b 100644 Binary files a/docs/reference/Rplot001.png and b/docs/reference/Rplot001.png differ diff --git a/docs/reference/SET_MCMC.html b/docs/reference/SET_MCMC.html index c6e8af8..c63be90 100644 --- a/docs/reference/SET_MCMC.html +++ b/docs/reference/SET_MCMC.html @@ -17,7 +17,7 @@+-## [1] FALSE
The
-lag
method forlistw
objects is often used to create spatially lagged values, and returns the same values as the vector given by the product of the sparse general matrix and an input numeric vector. Note that by settingzero.policy
toTRUE
, the spatial lag of entity 21, which has no neighbours, is zero by construction:+The
+lag
method forlistw
objects is often +used to create spatially lagged values, and returns the same values as +the vector given by the product of the sparse general matrix and an +input numeric vector. Note that by settingzero.policy
to +TRUE
, the spatial lag of entity 21, which has no +neighbours, is zero by construction:+all.equal(r1, r2, check.attributes=FALSE)set.seed(1) x <- runif(n) r1 <- as.numeric(W %*% x) r2 <- lag(nb_W, x, zero.policy=TRUE) -all.equal(r1, r2, check.attributes=FALSE)
-## [1] TRUE
+ -+c(x[21], r1[21])
## [1] 0.9347052 0.0000000
-@@ -1874,7 +1055,7 @@Log determinants (asymmetric weights) used in spatial regression +
Log determinants (asymmetric weights) used in spatial +regression
-Calculating the log determinant for asymmetric weights (here with symmetric neighbours and symmetry induced by non-constant numbers of neighbours) may be carried out using eigenvalues as before, but the result may be a complex vector (here it is not, as discussed below). The appropriate
-determinant
method for"dgCMatrix"
objects uses an LU decomposition internally:+Calculating the log determinant for asymmetric weights (here with +symmetric neighbours and symmetry induced by non-constant numbers of +neighbours) may be carried out using eigenvalues as before, but the +result may be a complex vector (here it is not, as discussed below). The +appropriate
+determinant
method for"dgCMatrix"
+objects uses an LU decomposition internally:-## [1] -1.594376
+class(I - rho * W)
-## [1] "dgCMatrix" ## attr(,"package") ## [1] "Matrix"
@@ -385,6 +385,7 @@+c(determinant(I - rho * W, logarithm=TRUE)$modulus)
## [1] -1.594376
We can show the internal workings of the method as:
-diff --git a/docs/reference/ML_models.html b/docs/reference/ML_models.html index 6074b24..b1b5c37 100644 --- a/docs/reference/ML_models.html +++ b/docs/reference/ML_models.html @@ -29,7 +29,7 @@+@@ -325,55 +454,76 @@## [1] -1.594376
Log dete
Log determinants: symmetric by similarity
-The
-nb2listw
function stores components that can be employed to transform the asymmetric weights matrix to symmetry by similarity, permitting the same log determinant to be computed using less costly numerical methods. The"W"
style used the cardinalities of neighbour sets (row sums) to introduce row standardisation, and they are stored as an attribute:+The
+ +all.equal(d, spdep::card(col2))nb2listw
function stores components that can be +employed to transform the asymmetric weights matrix to symmetry by +similarity, permitting the same log determinant to be computed using +less costly numerical methods. The"W"
style used the +cardinalities of neighbour sets (row sums) to introduce row +standardisation, and they are stored as an attribute:-## [1] TRUE
If we first restore the row-standarised matrix to its binary form (which must be symmetric), we can pre- and post-multiply by the square roots of the inverted neighbour counts, yielding a symmetric matrix with the appropriate properties:
-@@ -359,7 +250,7 @@+If we first restore the row-standarised matrix to its binary form +(which must be symmetric), we can pre- and post-multiply by the square +roots of the inverted neighbour counts, yielding a symmetric matrix with +the appropriate properties:
+-## [1] TRUE
+-## [1] Inf
+ --## 21 -## 0
@@ -165,17 +165,17 @@++## [1] 0
class(Ws)
-## [1] "dsCMatrix" ## attr(,"package") ## [1] "Matrix"
+c(determinant(I - rho * Ws, logarithm=TRUE)$modulus)
-## [1] -1.594376
As can be seen, the division by the square root of zero for entity 21 is not a problem as the row of
+dW
is zero. The transformation by similarity permits the use of numerical methods for sparse symmetric matrices (and equivalently for eigenvalues and dense matrices). Note that this transformation is not available for intrinsically asymmetric neighbours, or for intrinsically asymmetric general weights.As can be seen, the division by the square root of zero for entity 21 +is not a problem as the row of
dW
is zero. The +transformation by similarity permits the use of numerical methods for +sparse symmetric matrices (and equivalently for eigenvalues and dense +matrices). Note that this transformation is not available for +intrinsically asymmetric neighbours, or for intrinsically asymmetric +general weights.-diff --git a/docs/reference/ME.html b/docs/reference/ME.html index 7ee206b..8f190b5 100644 --- a/docs/reference/ME.html +++ b/docs/reference/ME.html @@ -17,7 +17,7 @@Using
eigs
in RSpectra for finding some eigenvalues +Using
-eigs
in RSpectra for finding +some eigenvaluesIn spatial regression, the domain of the spatial coefficient is given by the inverse of the maximum and minimum eigenvalues. When
-n
is moderate, we have the eigenvalues anyway, so the interval for line search is available without extra effort. Whenn
is somewhat larger, use may be made of theeigs
function in RSpectra:+In spatial regression, the domain of the spatial coefficient is given +by the inverse of the maximum and minimum eigenvalues. When +
+n
is moderate, we have the eigenvalues anyway, so the +interval for line search is available without extra effort. When +n
is somewhat larger, use may be made of the +eigs
function in RSpectra:-## [1] -0.3212551 0.1638329
@@ -883,7 +320,7 @@+-if (!require("RSpectra", quietly=TRUE)) dothis <- FALSE
+## [1] -0.3212551 0.1638329
In this case, the results are trivial with small
-k
.+-## [1] -1.544645 1.000000
+-## [1] -1.544645 1.000000
Using row-standardisation has the nice feature of setting the upper bound to unity, and there are graph methods for finding out whether the lower bound is
+-1
.Using row-standardisation has the nice feature of setting the upper +bound to unity, and there are graph methods for finding out whether the +lower bound is
-1
.@@ -382,17 +532,18 @@Using igraph
Converting from symmetric adjacency matrix to graph
-First we’ll see how to get from sparse matrices to graphs. The mode of a symmetric matrix is
-"undirected"
by definition:@@ -138,6 +138,7 @@+First we’ll see how to get from sparse matrices to graphs. The mode +of a symmetric matrix is
+"undirected"
by definition:class(B)
-## [1] "dgCMatrix" ## attr(,"package") ## [1] "Matrix"
@@ -70,7 +70,7 @@+object.size(B)
- +## 10824 bytes
## Loading required package: igraph
@@ -402,41 +553,48 @@## ## Attaching package: 'igraph'
Converting from sym
-## The following object is masked from 'package:base': ## ## union
diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 0103057..5976ad2 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,9 +1,9 @@ -pandoc: 2.14.0.3 -pkgdown: 2.0.6 +pandoc: 3.1.3 +pkgdown: 2.0.7 pkgdown_sha: ~ articles: SpatialFiltering: SpatialFiltering.html nb_igraph: nb_igraph.html sids_models: sids_models.html -last_built: 2022-10-07T18:46Z +last_built: 2023-11-23T10:35Z diff --git a/docs/reference/GMerrorsar-1.png b/docs/reference/GMerrorsar-1.png index 1c44f0d..b35453b 100644 Binary files a/docs/reference/GMerrorsar-1.png and b/docs/reference/GMerrorsar-1.png differ diff --git a/docs/reference/GMerrorsar.html b/docs/reference/GMerrorsar.html index 1deb508..9b81b7a 100644 --- a/docs/reference/GMerrorsar.html +++ b/docs/reference/GMerrorsar.html @@ -17,7 +17,7 @@-g1 <- graph.adjacency(B, mode="undirected") +
+g1 <- graph.adjacency(B, mode="undirected") class(g1)
-## [1] "igraph"
+-object.size(g1)
+## 9216 bytes
## 6544 bytes
Converting from graph to symmetric adjacency matrix
-We can also convert this graph pack to the same matrix, but note that
-get.adjacency
chooses a particular class of sparse matrix to be returned, so that the conversion process typically leads many matrices to fewer graph types, and back to fewer matrix types:+We can also convert this graph pack to the same matrix, but note that +
+get.adjacency
chooses a particular class of sparse matrix +to be returned, so that the conversion process typically leads many +matrices to fewer graph types, and back to fewer matrix types:# Matrix 1.4-2 vulnerability work-around ow <- options("warn")$warn options("warn"=2L) -B1 <- try(get.adjacency(g1), silent=TRUE) +B1 <- try(get.adjacency(g1), silent=TRUE) if (!inherits(B1, "try-error")) print(class(B1))
-## [1] "dgCMatrix" ## attr(,"package") ## [1] "Matrix"
- ++if (!inherits(B1, "try-error")) print(object.size(B1))
- +## 10824 bytes
-## [1] TRUE
+options("warn"=ow)
@@ -63,7 +63,24 @@Graph components in spdep
-A simple example of using igraph to do the same as an existing spdep function is Nicholas Lewin-Koh’s
-n.comp.nb
from the early days of the package. It is useful to know whether annb
object is divided up into separate subgraphs, and which entities are members of which such subgraph.@@ -149,7 +143,7 @@+A simple example of using igraph to do the same as +an existing spdep function is Nicholas Lewin-Koh’s +
+n.comp.nb
from the early days of the package. It is useful +to know whether annb
object is divided up into separate +subgraphs, and which entities are members of which such subgraph.## @@ -446,38 +604,47 @@
Graph components in spdep
Graph components in igraph
-The same result can be obtained using the
-clusters
function in igraph:-c1 <- clusters(g1) +
The same result can be obtained using the
+clusters
+function in igraph:+c1 <- clusters(g1) c1$no == res$nc
-## [1] TRUE
+-all.equal(c1$membership, res$comp.id)
+all.equal(c1$membership, res$comp.id)
- +## [1] "names for target but not for current"
## [1] TRUE
The same holds for the row-standardised variant:
-++g1W <- graph.adjacency(W, mode="directed", weighted="W") +c1W <- clusters(g1W) +all.equal(c1W$membership, res$comp.id)W <- as(spdep::nb2listw(col2, style="W", zero.policy=TRUE), "CsparseMatrix") -g1W <- graph.adjacency(W, mode="directed", weighted="W") -c1W <- clusters(g1W) -all.equal(c1W$membership, res$comp.id)
## [1] "names for target but not for current"
@@ -86,7 +86,7 @@Shortest paths in weights matrices: igraph
-Finding shortest paths between spatial entities across a given graph is a way to express closeness. If the graph is connected, that is that it is possible to traverse the graph edges from any node to any other, the longest shortest path is then a useful measure. In igraph, the
-is.connected
function tells us tells us that our graph is not connected, as we know from the figure above. The diameter measure is then the diameter of the largest component subgraph. Note that this generates ann
xn
matrix:+-is.connected(g1)
Finding shortest paths between spatial entities across a given graph +is a way to express closeness. If the graph is connected, that is that +it is possible to traverse the graph edges from any node to any other, +the longest shortest path is then a useful measure. In +igraph, the
+is.connected
function tells us +tells us that our graph is not connected, as we know from the figure +above. The diameter measure is then the diameter of the largest +component subgraph. Note that this generates ann
x +n
matrix:+is.connected(g1)
-## [1] FALSE
diff --git a/docs/index.html b/docs/index.html index 8803c65..8dfd9b2 100644 --- a/docs/index.html +++ b/docs/index.html @@ -33,7 +33,7 @@-dg1 <- diameter(g1) +
+dg1 <- diameter(g1) dg1
-## [1] 7
-sp_mat <- shortest.paths(g1) +
+sp_mat <- shortest.paths(g1) str(sp_mat)
## num [1:49, 1:49] 0 1 1 2 2 3 4 3 3 4 ... ## - attr(*, "dimnames")=List of 2 @@ -487,45 +654,79 @@
Shortest paths in weights mat
Shortest paths in weights matrices: spdep
-If we do the same in spdep, using
-nblag
to a maximum number of lag orders - the diameter, but which is unknown in advance (the largest lag order for which the number of links is greater than zero), we run into the problem of how to represent missing neighbour information.+If we do the same in spdep, using
+nblag
+to a maximum number of lag orders - the diameter, but which is unknown +in advance (the largest lag order for which the number of links is +greater than zero), we run into the problem of how to represent missing +neighbour information.nbl10 <- spdep::nblag(col2, maxlag=10) vals <- sapply(nbl10, function(x) sum(spdep::card(x))) zero <- which(vals == 0) zero[1]-1
-## [1] 7
If we insert zero into the weights matrix where there is no connection using
-zero.policy=TRUE
, we generate a zero shortest path. If we are to create a matrix that matches the one produced byshortest.paths
, we need to set all these non-structural zeros to infinity (the length of the path between unconnected nodes), and re-instate structural zeros on the diagonal:+If we insert zero into the weights matrix where there is no +connection using
+zero.policy=TRUE
, we generate a zero +shortest path. If we are to create a matrix that matches the one +produced byshortest.paths
, we need to set all these +non-structural zeros to infinity (the length of the path between +unconnected nodes), and re-instate structural zeros on the diagonal:+all.equal(mat, sp_mat, check.attributes=FALSE)lmat <- lapply(nbl10[1:(zero[1]-1)], spdep::nb2mat, style="B", zero.policy=TRUE) mat <- matrix(0, n, n) for (i in seq(along=lmat)) mat = mat + i*lmat[[i]] mat[mat==0] <- Inf diag(mat) <- 0 -all.equal(mat, sp_mat, check.attributes=FALSE)
## [1] TRUE
@@ -152,22 +152,6 @@Smirnov/Anselin (2009) cyclical matrices
-Another area in which a graph representation might prove useful is in trying to establish the domain of the spatial coefficient when spatial weights are row-standardised. In that case by construction we know that the maximum eigenvalue is 1. If there are multiple blocks, that is graph components, where the numbers of nodes per block are greater than 1, then each will have a maximum eigenvalue of 1. The remaining problems are the numbers of zero eigenvalues (at least the singleton graph components), and whether any non-singleton component fulfills the condition termed by Smirnov and Anselin (2009) a cyclical matrix, for which the minimum eigenvalue is -1. The term cyclical appears to be used in many different ways, and it is not clear that its use here after Smirnov and Anselin (2009, pp. 2984-2985) indicates which meaning should be used to find the relevant graph function. The definition used here is that a block matrix (subgraph) is cyclical if: “for every location, every pair of its neighbours are not connected.” That is, if w[i,j] and w[i,k] are greater than zero, w[j,k] must be zero to meet the condition.
-The internal function find_q1_q2 returns the number of non-singleton components, and the number of these that meet this condition. It does this for each block/subgraph by testing the condition until it meets w[j,k] > 0, at which point it breaks. Smirnov and Anselin (2009) state that rook neighbours on a regular grid meet the condition:
-@@ -85,7 +85,8 @@+Another area in which a graph representation might prove useful is in +trying to establish the domain of the spatial coefficient when spatial +weights are row-standardised. In that case by construction we know that +the maximum eigenvalue is 1. If there are multiple blocks, that is graph +components, where the numbers of nodes per block are greater than 1, +then each will have a maximum eigenvalue of 1. The remaining problems +are the numbers of zero eigenvalues (at least the singleton graph +components), and whether any non-singleton component fulfills the +condition termed by Smirnov and Anselin (2009) a cyclical matrix, for +which the minimum eigenvalue is -1. The term cyclical appears to be used +in many different ways, and it is not clear that its use here after +Smirnov and Anselin (2009, pp. 2984-2985) indicates which meaning should +be used to find the relevant graph function. The definition used here is +that a block matrix (subgraph) is cyclical if: “for every location, +every pair of its neighbours are not connected.” That is, if w[i,j] and +w[i,k] are greater than zero, w[j,k] must be zero to meet the +condition.
+The internal function find_q1_q2 returns the number of non-singleton +components, and the number of these that meet this condition. It does +this for each block/subgraph by testing the condition until it meets +w[j,k] > 0, at which point it breaks. Smirnov and Anselin (2009) +state that rook neighbours on a regular grid meet the condition:
+nb_r <- spdep::cell2nb(7, 7, type="rook") nb_rW <- spdep::nb2listw(nb_r, style="W") spdep:::find_q1_q2(nb_rW)
-## [1] 1 1
One block/graph component is found, and this one meets the cyclical matrix condition, as also shown by the domain:
-@@ -547,7 +748,7 @@+One block/graph component is found, and this one meets the cyclical +matrix condition, as also shown by the domain:
+1/range(Re(eigenw(similar.listw(nb_rW))))
-## [1] -1 1
This does not apply to the spatial weights we have been using above, with two non-singleton components, neither meeting the cyclical matrix condition:
-+This does not apply to the spatial weights we have been using above, +with two non-singleton components, neither meeting the cyclical matrix +condition:
+spdep:::find_q1_q2(nb_W)
-## [1] 2 0
+1/range(Re(eigenw(similar.listw(nb_W))))
-## [1] -1.544645 1.000000
By construction, all two-node connected graph components also meet the condition, as the eigenvalues sum to zero, and the maximum is unity, so the minimum must be -1.
+By construction, all two-node connected graph components also meet +the condition, as the eigenvalues sum to zero, and the maximum is unity, +so the minimum must be -1.
Smirnov/Anselin (2009) cyclical m diff --git a/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-13-1.png index 35ad15d..68226ef 100644 Binary files a/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-13-1.png differ diff --git a/docs/articles/sids_models.html b/docs/articles/sids_models.html index 3490b0f..0a8cecc 100644 --- a/docs/articles/sids_models.html +++ b/docs/articles/sids_models.html @@ -33,7 +33,7 @@
diff --git a/docs/articles/sids_models_files/figure-html/unnamed-chunk-13-1.png b/docs/articles/sids_models_files/figure-html/unnamed-chunk-13-1.png index c25e312..79e8f99 100644 Binary files a/docs/articles/sids_models_files/figure-html/unnamed-chunk-13-1.png and b/docs/articles/sids_models_files/figure-html/unnamed-chunk-13-1.png differ diff --git a/docs/articles/sids_models_files/figure-html/unnamed-chunk-14-1.png b/docs/articles/sids_models_files/figure-html/unnamed-chunk-14-1.png index f38d076..94021a1 100644 Binary files a/docs/articles/sids_models_files/figure-html/unnamed-chunk-14-1.png and b/docs/articles/sids_models_files/figure-html/unnamed-chunk-14-1.png differ diff --git a/docs/authors.html b/docs/authors.html index e423470..89dfe06 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -17,7 +17,7 @@-North Carolina SIDS data set (models)
-Roger Bivand
+Roger +Bivand
Source:vignettes/sids_models.Rmd
@@ -98,12 +99,37 @@Roger Bivand
Introduction
-This data set was presented first in Symons, Grimson, and Yuan (1983), analysed with reference to the spatial nature of the data in Cressie and Read (1985), expanded in Cressie and Chan (1989), and used in detail in Cressie (1991). It is for the 100 counties of North Carolina, and includes counts of numbers of live births (also non-white live births) and numbers of sudden infant deaths, for the July 1, 1974 to June 30, 1978 and July 1, 1979 to June 30, 1984 periods. In Cressie and Read (1985), a listing of county neighbours based on shared boundaries (contiguity) is given, and in Cressie and Chan (1989), and in Cressie (1991, 386–89), a different listing based on the criterion of distance between county seats, with a cutoff at 30 miles. The county seat location coordinates are given in miles in a local (unknown) coordinate reference system. The data are also used to exemplify a range of functions in the spatial statistics module user’s manual (Kaluzny et al. 1996).
+This data set was presented first in Symons, +Grimson, and Yuan (1983), +analysed with reference to the spatial nature of the data in Cressie and Read (1985), expanded in Cressie and Chan (1989), and used in detail in +Cressie (1991). It is for the 100 counties of +North Carolina, and includes counts of numbers of live births (also +non-white live births) and numbers of sudden infant deaths, for the July +1, 1974 to June 30, 1978 and July 1, 1979 to June 30, 1984 periods. In +Cressie and Read (1985), a listing of county +neighbours based on shared boundaries (contiguity) is given, and in +Cressie and Chan (1989), and in Cressie (1991, +386–89), a different listing based on the criterion of +distance between county seats, with a cutoff at 30 miles. The county +seat location coordinates are given in miles in a local (unknown) +coordinate reference system. The data are also used to exemplify a range +of functions in the spatial statistics module user’s manual (Kaluzny et al. +1996).
Getting the data into R
-We will be using the spdep and spatialreg packages, here version: spdep, version 1.2-7, 2022-10-01, the sf package and the tmap package. The data from the sources referred to above is documented in the help page for the
+nc.sids
data set in spData. The actual data, included in a shapefile of the county boundaries for North Carolina were made available in the maptools package.1 These data are known to be geographical coordinates (longitude-latitude in decimal degrees) and are assumed to use the NAD27 datum.We will be using the spdep and +spatialreg packages, here version: spdep, version +1.3-1, 2023-11-03, the sf package and the +tmap package. The data from the sources referred to +above is documented in the help +page for the
nc.sids
data set in +spData. The actual data, included in a shapefile of the +county boundaries for North Carolina were made available in the +maptools package 1. These data are known to be geographical +coordinates (longitude-latitude in decimal degrees) and are assumed to +use the NAD27 datum.library(spdep) nc <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE) @@ -113,15 +139,26 @@
Getting the data into Rnc$ft.SID74 <- sqrt(1000)*(sqrt(nc$SID74/nc$BIR74) + sqrt((nc$SID74+1)/nc$BIR74)) nc$both <- factor(paste(nc$L_id, nc$M_id, sep=":"))
-## [1] TRUE
We will now examine the data set reproduced from Cressie and collaborators, included in spData (formerly in spdep), and add the neighbour relationships used in Cressie and Chan (1989) to the background map as a graph shown in Figure :
+We will now examine the data set reproduced from Cressie and +collaborators, included in spData (formerly in +spdep), and add the neighbour relationships used in +Cressie and Chan (1989) to the background map as +a graph shown in Figure \(\ref{plot-CC89.nb}\):
gal_file <- system.file("weights/ncCC89.gal", package="spData")[1] ncCC89 <- read.gal(gal_file, region.id=nc$FIPSNO)
+CAR model fitting
-We will now try to replicate three of the four models fitted by (Cressie and Chan 1989) to the transformed rates variable. The first thing to do is to try to replicate their 30 mile distance between county seats neighbours, which almost works. From there we try to reconstruct three of the four models they fit, concluding that we can get quite close, but that a number of questions are raised along the way.
-Building the weights is much more complicated, because they use a combination of distance-metric and population-at-risk based weights, but we can get quite close (see also Kaluzny et al. 1996):
+We will now try to replicate three of the four models fitted by (Cressie and Chan +1989) to the transformed rates variable. The first thing to +do is to try to replicate their 30 mile distance between county seats +neighbours, which almost works. From there we try to reconstruct three +of the four models they fit, concluding that we can get quite close, but +that a number of questions are raised along the way.
+Building the weights is much more complicated, because they use a +combination of distance-metric and population-at-risk based weights, but +we can get quite close (see also Kaluzny et al. 1996):
-sids.nhbr30.dist <- nbdists(ncCC89, cbind(nc$east, nc$north)) sids.nhbr <- listw2sn(nb2listw(ncCC89, glist=sids.nhbr30.dist, style="B", zero.policy=TRUE)) @@ -130,11 +167,21 @@
CAR model fittingel1 <- min(dij)/dij el2 <- sqrt(n[sids.nhbr$to]/n[sids.nhbr$from]) sids.nhbr$weights <- el1*el2 -sids.nhbr.listw <- sn2listw(sids.nhbr)
The first model (I) is a null model with just an intercept, the second (II) includes all the 12 parcels of contiguous counties in 4 east-west and 4 north-south bands, while the fourth (IV) includes the transformed non-white birth-rate:
+if (packageVersion("spdep") >= "1.3.1") { + sids.nhbr.listw <- sn2listw(sids.nhbr, style="B", zero.policy=TRUE) +} else { + sids.nhbr.listw <- sn2listw(sids.nhbr) +}The first model (I) is a null model with just an intercept, the +second (II) includes all the 12 parcels of contiguous counties in 4 +east-west and 4 north-south bands, while the fourth (IV) includes the +transformed non-white birth-rate:
-Cressie identifies Anson county as an outlier, and drops it from further analysis. Because the weights are constructed in a complicated way, they will be subsetted by dropping the row and column of the weights matrix:
+Cressie identifies Anson county as an outlier, and drops it from +further analysis. Because the weights are constructed in a complicated +way, they will be subsetted by dropping the row and column of the +weights matrix:
-lm_nc <- lm(ft.SID74 ~ 1, data=nc) outl <- which.max(rstandard(lm_nc)) @@ -145,7 +192,12 @@
CAR model fittingW.4 <- W[-outl, -outl] sids.nhbr.listw.4 <- mat2listw(W.4) nc2 <- nc[!(1:length(nc$CNTY_ID) %in% outl),]
It appears that both numerical issues (convergence in particular) and uncertainties about the exact spatial weights matrix used make it difficult to reproduce the results of Cressie and Chan (1989), also given in Cressie (1991). We now try to replicate them for the null weighted CAR model (Cressie has intercept 2.838, \(\hat{\theta}\) 0.833, for k=1):
+It appears that both numerical issues (convergence in particular) and +uncertainties about the exact spatial weights matrix used make it +difficult to reproduce the results of Cressie and +Chan (1989), also given in +Cressie (1991). We now try to replicate them +for the null weighted CAR model (Cressie has intercept 2.838, \(\hat{\theta}\) 0.833, for k=1):
-library(spatialreg) ecarIaw <- spautolm(ft.SID74 ~ 1, data=nc2, listw=sids.nhbr.listw.4, weights=BIR74, family="CAR") @@ -169,8 +221,10 @@
CAR model fitting## ML residual variance (sigma squared): 1266.5, (sigma: 35.588) ## Number of observations: 99 ## Number of parameters estimated: 3 -## AIC: 243.69
The spatial parcels model seems not to work, with Cressie’s \(\hat{\theta}\) 0.710, and failure in the numerical Hessian used to calculate the standard error of the spatial coefficient:
+## AIC: NA (not available for weighted model) +The spatial parcels model seems not to work, with Cressie’s \(\hat{\theta}\) 0.710, and failure in the +numerical Hessian used to calculate the standard error of the spatial +coefficient:
@@ -205,8 +259,12 @@ecarIIaw <- spautolm(ft.SID74 ~ both - 1, data=nc2, listw=sids.nhbr.listw.4, weights=BIR74, family="CAR") summary(ecarIIaw)
CAR model fitting## ML residual variance (sigma squared): 891.48, (sigma: 29.858) ## Number of observations: 99 ## Number of parameters estimated: 14 -## AIC: 226.51 -
Finally, the non-white model repeats Cressie’s finding that much of the variance of the transformed SIDS rate for 1974–8 can be accounted for by the transformed non-white birth variable (Cressie intercept 1.644, \(\hat{b}\) 0.0346, \(\hat{\theta}\) 0.640 — not significant), but the estimate of the spatial coefficient is not close here:
+## AIC: NA (not available for weighted model) +Finally, the non-white model repeats Cressie’s finding that much of +the variance of the transformed SIDS rate for 1974–8 can be accounted +for by the transformed non-white birth variable (Cressie intercept +1.644, \(\hat{b}\) 0.0346, \(\hat{\theta}\) 0.640 — not significant), +but the estimate of the spatial coefficient is not close here:
@@ -221,23 +279,24 @@ecarIVaw <- spautolm(ft.SID74 ~ ft.NWBIR74, data=nc2, listw=sids.nhbr.listw.4, weights=BIR74, family="CAR") summary(ecarIVaw)
CAR model fitting## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) -## (Intercept) 1.4371519 0.2252729 6.3796 1.775e-10 +## (Intercept) 1.4371520 0.2252729 6.3796 1.775e-10 ## ft.NWBIR74 0.0408354 0.0062919 6.4902 8.572e-11 ## ## Lambda: 0.22391 LR test value: 1.1577 p-value: 0.28194 -## Numerical Hessian standard error of lambda: 0.55128 +## Numerical Hessian standard error of lambda: 0.55147 ## ## Log likelihood: -114.0376 ## ML residual variance (sigma squared): 1201.5, (sigma: 34.663) ## Number of observations: 99 ## Number of parameters estimated: 4 -## AIC: 236.08 +## AIC: NA (not available for weighted model)
-nc2$fitIV <- fitted.values(ecarIVaw) sf_use_s2(FALSE) tm_shape(nc2) + tm_fill("fitIV")
The final figure shows the value of the log likelihood function for the null model (I):
+The final figure shows the value of the log likelihood function for +the null model (I):
@@ -252,23 +311,32 @@ecarIawll <- spautolm(ft.SID74 ~ 1, data=nc2, listw=sids.nhbr.listw.4, weights=BIR74, family="CAR", llprof=seq(-0.1, 0.9020532358, length.out=100)) plot(ll ~ lambda, ecarIawll$llprof, type="l")
References -Cressie, N., and N. H. Chan. 1989. “Spatial Modelling of Regional Variables.” Journal of the American Statistical Association 84: 393–401. +Cressie, N., and N. H. Chan. 1989. “Spatial Modelling of Regional +Variables.” Journal of the American Statistical +Association 84: 393–401.
-Cressie, N., and T. R. C. Read. 1985. “Do Sudden Infant Deaths Come in Clusters?” Statistics and Decisions Supplement Issue 2: 333–49. +Cressie, N., and T. R. C. Read. 1985. “Do Sudden Infant Deaths +Come in Clusters?” Statistics and Decisions Supplement +Issue 2: 333–49.-Kaluzny, S. P., S. C. Vega, T. P. Cardoso, and A. A. Shelly. 1996. S-PLUS SPATIALSTATS User’s Manual Version 1.0. Seattle: MathSoft Inc. +Kaluzny, S. P., S. C. Vega, T. P. Cardoso, and A. A. Shelly. 1996. +S-PLUS SPATIALSTATS User’s Manual Version 1.0. +Seattle: MathSoft Inc.-Symons, M. J., R. C. Grimson, and Y. C. Yuan. 1983. “Clustering of Rare Events.” Biometrics 39: 193–205. +Symons, M. J., R. C. Grimson, and Y. C. Yuan. 1983. “Clustering of +Rare Events.” Biometrics 39: 193–205.+@@ -290,7 +358,7 @@
-
- +
These data were taken with permission from a now-offline link: [sal.agecon.uiuc.edu/datasets/sids.zip]; see also GeoDa Center for a contemporary source.↩︎
These data were taken with permission from a now-offline +link: [sal.agecon.uiuc.edu/datasets/sids.zip]; see also GeoDa Center for +a contemporary source.↩︎
References -
Site built with pkgdown 2.0.6.
+Site built with pkgdown 2.0.7.
Authors
Michael Tiefelsdorf. Contributor.
-- -Bendix Carstensen. Contributor. -
-- -Gregory Warnes. Contributor. -
-- -Soren Hojsgaard. Contributor. -
-- Randal Johnson. Contributor. -
-@@ -230,7 +242,7 @@@@ -177,7 +161,11 @@-Citation
Bivand, R., G. Millo, and G. Piras. 2021. A Review of Software for Spatial Econometrics in R. Mathematics 9 (11):1276. https://doi.org/10.3390/math9111276.
+Bivand R, Millo G, Piras G (2021). +“A Review of Software for Spatial Econometrics in R.” +Mathematics, 9(11). +doi:10.3390/math9111276, https://www.mdpi.com/2227-7390/9/11/1276. +
@Article{, title = {A Review of Software for Spatial Econometrics in {R}}, author = {Roger Bivand and Giovanni Millo and Gianfranco Piras}, @@ -188,7 +176,11 @@-Citation
url = {https://www.mdpi.com/2227-7390/9/11/1276}, doi = {10.3390/math9111276}, }Roger Bivand, Gianfranco Piras (2015). Comparing Implementations of Estimation Methods for Spatial Econometrics. Journal of Statistical Software, 63(18), 1-36. URL https://www.jstatsoft.org/v63/i18/.
+Bivand R, Piras G (2015). +“Comparing Implementations of Estimation Methods for Spatial Econometrics.” +Journal of Statistical Software, 63(18), 1–36. +doi:10.18637/jss.v063.i18. +
@Article{, title = {Comparing Implementations of Estimation Methods for Spatial Econometrics}, author = {Roger Bivand and Gianfranco Piras}, @@ -199,7 +191,11 @@-Citation
pages = {1--36}, doi = {10.18637/jss.v063.i18}, }Bivand, R. S., Hauke, J., and Kossowski, T. (2013). Computing the Jacobian in Gaussian spatial autoregressive models: An illustrated comparison of available methods. Geographical Analysis, 45(2), 150-179. URL https://doi.org/10.1111/gean.12008
+Bivand R, Hauke J, Kossowski T (2013). +“Computing the Jacobian in Gaussian spatial autoregressive models: An illustrated comparison of available methods.” +Geographical Analysis, 45(2), 150–179. +doi:10.1111/gean.12008. +
@Article{, title = {Computing the Jacobian in Gaussian spatial autoregressive models: An illustrated comparison of available methods}, author = {Roger Bivand and Jan Hauke and Tomasz Kossowski}, @@ -210,14 +206,30 @@-Citation
pages = {150--179}, doi = {10.1111/gean.12008}, }Roger S. Bivand, Edzer Pebesma, Virgilio Gomez-Rubio, 2013. Applied spatial data analysis with R, Second edition. Springer, NY. https://asdar-book.org/
+Bivand R, Pebesma E, Gómez-Rubio V (2013). +Applied spatial data analysis with R, Second edition. +Springer, NY. +https://asdar-book.org/. +
@Book{, - author = {Roger S. Bivand and Edzer Pebesma and Virgilio Gomez-Rubio}, + author = {Roger S. Bivand and Edzer Pebesma and Virgilio Gómez-Rubio}, title = {Applied spatial data analysis with {R}, Second edition}, year = {2013}, publisher = {Springer, NY}, url = {https://asdar-book.org/}, }+Pebesma E, Bivand R (2023). +Spatial Data Science With Applications in R. +Chapman & Hall. +https://r-spatial.org/book/. +
+@Book{, + author = {Edzer Pebesma and Roger S. Bivand}, + title = {Spatial Data Science With Applications in {R}}, + year = {2023}, + publisher = {Chapman & Hall}, + url = {https://r-spatial.org/book/}, +}Citation
Dev status
diff --git a/docs/news/index.html b/docs/news/index.html index 7849cce..870b0f3 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -17,7 +17,7 @@Changelog
Source:NEWS.md
++Version 1.3-1 (development)
+
- +
move
expm
from Imports to Suggests #42- +
added
zero.policy
pass-through tospdep::mat2listw
calls inpredict.Sarlm
and tospdep::sn2listw
insids_models.Rmd
; setspdep
requirement to1.3-1
- +
corrected #19 because the fitted model weights component is never NULL, but may have a single unique value
++Version 1.2-9 (2023-05-25)2023-05-25
+
- +
- +
address #37; #38 remains (no formula Durbin support for prediction using any Sarlm object)
- +
address #19 by not reporting
AIC
where case weights are used inspautolm
orerrorsarlm
- +
address bug in
predict()
for new data, SDEM. Others in #37, #38 need work.- +
Further added checking for SLX/SDEM impacts and edge/corner cases; starting transition to use multcomp in place og gmodels
+Version 1.2-6 (2022-10-07)2022-10-07
- @@ -101,7 +118,8 @@
make local copy of
gmodels::estimable()
for lm objects only, add authors to package contributorsVersi
- #2 Split spatialreg from spdep; spdep functions still present there as deprecated
+Spatial simultaneous autoregressive error model estimation by GMM
diff --git a/docs/reference/MCMCsamp.html b/docs/reference/MCMCsamp.html index d4c23fe..20013ef 100644 --- a/docs/reference/MCMCsamp.html +++ b/docs/reference/MCMCsamp.html @@ -17,7 +17,7 @@GMerrorsar(formula, data = list(), listw, na.action = na.fail, - zero.policy = NULL, method="nlminb", arnoldWied=FALSE, + zero.policy = attr(listw, "zero.policy"), method="nlminb", arnoldWied=FALSE, control = list(), pars, scaleU=FALSE, verbose=NULL, legacy=FALSE, se.lambda=TRUE, returnHcov=FALSE, pWOrder=250, tol.Hcov=1.0e-10) # S3 method for Gmsar @@ -297,9 +297,9 @@
Examples
#> coef(x) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 59.8932192 5.36616252 11.161276 0.0000000000 -#> INC -0.9413120 0.33056857 -2.847554 0.0044056569 -#> HOVAL -0.3022502 0.09047605 -3.340665 0.0008357788 +#> (Intercept) 59.8932190 5.36616256 11.161276 0.0000000000 +#> INC -0.9413120 0.33056857 -2.847554 0.0044056574 +#> HOVAL -0.3022502 0.09047605 -3.340665 0.0008357787 COL.errW.GM <- GMerrorsar(CRIME ~ INC + HOVAL, data=COL.OLD, spdep::nb2listw(COL.nb, style="W"), returnHcov=TRUE) (x <- summary(COL.errW.GM, Hausman=TRUE)) @@ -371,6 +371,7 @@Examples
identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10)) #> [1] TRUE nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY)) +#> Warning: style is M (missing); style should be set to a valid value listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") @@ -392,7 +393,7 @@Examples
#> PCTOWNHOME -0.419890 0.191329 -2.1946 0.0281930 #> #> Lambda: 0.040487 LR test value: 5.2438 p-value: 0.022026 -#> Numerical Hessian standard error of lambda: 0.017209 +#> Numerical Hessian standard error of lambda: 0.017201 #> #> Log likelihood: -276.1069 #> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333) @@ -464,7 +465,7 @@Examples
Examples
#> [1] TRUE #require("spdep", quietly=TRUE) nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY)) +#> Warning: style is M (missing); style should be set to a valid value listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") @@ -159,7 +160,7 @@Examples
#> PCTOWNHOME -0.419890 0.191329 -2.1946 0.0281930 #> #> Lambda: 0.040487 LR test value: 5.2438 p-value: 0.022026 -#> Numerical Hessian standard error of lambda: 0.017209 +#> Numerical Hessian standard error of lambda: 0.017201 #> #> Log likelihood: -276.1069 #> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333) @@ -178,244 +179,44 @@Examples
#> 1. Empirical mean and standard deviation for each variable, #> plus standard error of the mean: #> -#> Mean SD Naive SE Time-series SE -#> lambda 0.04001 0.01707 0.0005399 0.001952 -#> (Intercept) -0.64306 0.18227 0.0057638 0.021874 -#> PEXPOSURE 0.07714 0.04402 0.0013921 0.005096 -#> PCTAGE65P 3.77115 0.66208 0.0209367 0.078911 -#> PCTOWNHOME -0.41005 0.20560 0.0065016 0.025344 +#> Mean SD Naive SE Time-series SE +#> lambda 0.04501 0.01711 0.000541 0.002089 +#> (Intercept) -0.62524 0.18031 0.005702 0.022321 +#> PEXPOSURE 0.07085 0.04851 0.001534 0.006552 +#> PCTAGE65P 3.74204 0.60659 0.019182 0.068848 +#> PCTOWNHOME -0.39873 0.21332 0.006746 0.025439 #> #> 2. Quantiles for each variable: #> #> 2.5% 25% 50% 75% 97.5% -#> lambda 0.003697 0.02876 0.04096 0.05108 0.07377 -#> (Intercept) -1.010162 -0.75178 -0.62651 -0.50945 -0.33213 -#> PEXPOSURE -0.012870 0.04708 0.08023 0.10886 0.15518 -#> PCTAGE65P 2.377758 3.33159 3.82741 4.16499 5.08140 -#> PCTOWNHOME -0.763610 -0.56285 -0.41434 -0.28717 0.04188 +#> lambda 0.009947 0.03345 0.04621 0.05710 0.07604 +#> (Intercept) -0.975566 -0.74111 -0.62524 -0.49983 -0.28177 +#> PEXPOSURE -0.005411 0.03271 0.06415 0.09889 0.17297 +#> PCTAGE65P 2.728875 3.32066 3.64211 4.14748 4.94058 +#> PCTOWNHOME -0.816134 -0.52762 -0.40539 -0.25284 0.01614 #> -# \dontrun{ +if (FALSE) { esar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8, family="SAR", method="eigen") summary(esar1fw) -#> -#> Call: -#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY, weights = POP8, family = "SAR", method = "eigen") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.48488 -0.26823 0.09489 0.46552 4.28343 -#> -#> Coefficients: -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.797063 0.144054 -5.5331 3.146e-08 -#> PEXPOSURE 0.080545 0.028334 2.8428 0.004473 -#> PCTAGE65P 3.816731 0.576037 6.6258 3.453e-11 -#> PCTOWNHOME -0.380778 0.156507 -2.4330 0.014975 -#> -#> Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764 -#> Numerical Hessian standard error of lambda: 0.016466 -#> -#> Log likelihood: -251.6017 -#> ML residual variance (sigma squared): 1104.1, (sigma: 33.229) -#> Number of observations: 281 -#> Number of parameters estimated: 6 -#> AIC: 515.2 -#> res <- MCMCsamp(esar1fw, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> lambda 0.01432 0.01645 0.0002327 0.0009528 -#> (Intercept) -0.82781 0.15037 0.0021266 0.0091370 -#> PEXPOSURE 0.08532 0.03057 0.0004323 0.0018015 -#> PCTAGE65P 3.82149 0.58470 0.0082689 0.0340386 -#> PCTOWNHOME -0.35000 0.16243 0.0022972 0.0091124 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> lambda -0.01695 0.003279 0.01484 0.0253 0.04706 -#> (Intercept) -1.13038 -0.918410 -0.81946 -0.7302 -0.53934 -#> PEXPOSURE 0.02943 0.064686 0.08326 0.1050 0.14704 -#> PCTAGE65P 2.61095 3.438180 3.84949 4.1896 4.93237 -#> PCTOWNHOME -0.67443 -0.464258 -0.34844 -0.2349 -0.03659 -#> ecar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="CAR", method="eigen") summary(ecar1f) -#> -#> Call: -#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY, family = "CAR", method = "eigen") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.539732 -0.384311 -0.030646 0.335126 3.808848 -#> -#> Coefficients: -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.648362 0.181129 -3.5796 0.0003442 -#> PEXPOSURE 0.077899 0.043692 1.7829 0.0745986 -#> PCTAGE65P 3.703830 0.627185 5.9055 3.516e-09 -#> PCTOWNHOME -0.382789 0.195564 -1.9574 0.0503053 -#> -#> Lambda: 0.084123 LR test value: 5.8009 p-value: 0.016018 -#> Numerical Hessian standard error of lambda: 0.030868 -#> -#> Log likelihood: -275.8283 -#> ML residual variance (sigma squared): 0.40758, (sigma: 0.63842) -#> Number of observations: 281 -#> Number of parameters estimated: 6 -#> AIC: 563.66 -#> res <- MCMCsamp(ecar1f, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> lambda 0.08622 0.02920 0.0004129 0.001653 -#> (Intercept) -0.68322 0.19789 0.0027987 0.011727 -#> PEXPOSURE 0.08282 0.05247 0.0007421 0.003270 -#> PCTAGE65P 3.75704 0.63493 0.0089793 0.036382 -#> PCTOWNHOME -0.35288 0.22991 0.0032515 0.014529 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> lambda 0.02522 0.0658 0.08841 0.1086 0.13455 -#> (Intercept) -1.08604 -0.8154 -0.67904 -0.5458 -0.30916 -#> PEXPOSURE -0.02079 0.0487 0.07909 0.1156 0.19020 -#> PCTAGE65P 2.57120 3.2835 3.75527 4.2142 5.01250 -#> PCTOWNHOME -0.82748 -0.4938 -0.34169 -0.2025 0.09923 -#> esar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8, family="SAR", method="eigen") summary(esar1fw) -#> -#> Call: -#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY, weights = POP8, family = "SAR", method = "eigen") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.48488 -0.26823 0.09489 0.46552 4.28343 -#> -#> Coefficients: -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.797063 0.144054 -5.5331 3.146e-08 -#> PEXPOSURE 0.080545 0.028334 2.8428 0.004473 -#> PCTAGE65P 3.816731 0.576037 6.6258 3.453e-11 -#> PCTOWNHOME -0.380778 0.156507 -2.4330 0.014975 -#> -#> Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764 -#> Numerical Hessian standard error of lambda: 0.016466 -#> -#> Log likelihood: -251.6017 -#> ML residual variance (sigma squared): 1104.1, (sigma: 33.229) -#> Number of observations: 281 -#> Number of parameters estimated: 6 -#> AIC: 515.2 -#> res <- MCMCsamp(esar1fw, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> lambda 0.01510 0.01693 0.0002394 0.001083 -#> (Intercept) -0.81406 0.15443 0.0021840 0.009162 -#> PEXPOSURE 0.08481 0.03188 0.0004509 0.001959 -#> PCTAGE65P 3.81810 0.55635 0.0078680 0.030659 -#> PCTOWNHOME -0.37343 0.16429 0.0023234 0.010419 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> lambda -0.01840 0.004274 0.0153 0.02593 0.04866 -#> (Intercept) -1.13018 -0.923852 -0.8057 -0.70461 -0.53477 -#> PEXPOSURE 0.02828 0.061177 0.0845 0.10433 0.15129 -#> PCTAGE65P 2.71673 3.436993 3.8196 4.18692 4.87663 -#> PCTOWNHOME -0.67435 -0.492926 -0.3811 -0.26100 -0.03220 -#> ecar1fw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8, family="CAR", method="eigen") summary(ecar1fw) -#> -#> Call: -#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY, weights = POP8, family = "CAR", method = "eigen") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.491042 -0.270906 0.081435 0.451556 4.198134 -#> -#> Coefficients: -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.790154 0.144862 -5.4545 4.910e-08 -#> PEXPOSURE 0.081922 0.028593 2.8651 0.004169 -#> PCTAGE65P 3.825858 0.577720 6.6223 3.536e-11 -#> PCTOWNHOME -0.386820 0.157436 -2.4570 0.014010 -#> -#> Lambda: 0.022419 LR test value: 0.38785 p-value: 0.53343 -#> Numerical Hessian standard error of lambda: 0.038543 -#> -#> Log likelihood: -251.5711 -#> ML residual variance (sigma squared): 1102.9, (sigma: 33.21) -#> Number of observations: 281 -#> Number of parameters estimated: 6 -#> AIC: 515.14 -#> res <- MCMCsamp(ecar1fw, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> lambda 0.04093 0.04223 0.0005973 0.002928 -#> (Intercept) -0.84314 0.19274 0.0027257 0.014578 -#> PEXPOSURE 0.09392 0.04040 0.0005713 0.003673 -#> PCTAGE65P 3.73730 0.58584 0.0082850 0.034082 -#> PCTOWNHOME -0.35679 0.18917 0.0026752 0.012137 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> lambda -0.04111 0.01262 0.04172 0.06947 0.12566 -#> (Intercept) -1.28030 -0.94566 -0.82057 -0.71741 -0.52485 -#> PEXPOSURE 0.02852 0.06854 0.08920 0.11372 0.19814 -#> PCTAGE65P 2.60227 3.33711 3.73379 4.11492 4.95998 -#> PCTOWNHOME -0.69743 -0.48868 -0.36696 -0.23354 0.05357 -#> -# } +} esar0 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(esar0) @@ -458,417 +259,53 @@Examples
#> plus standard error of the mean: #> #> Mean SD Naive SE Time-series SE -#> lambda 0.04646 0.01589 0.0005025 0.002026 -#> (Intercept) -0.62387 0.19521 0.0061731 0.028032 -#> PEXPOSURE 0.07337 0.03977 0.0012575 0.004163 -#> PCTAGE65P 3.77748 0.67223 0.0212579 0.078551 -#> PCTOWNHOME -0.42259 0.19976 0.0063168 0.027140 +#> lambda 0.04276 0.01673 0.0005291 0.002137 +#> (Intercept) -0.64231 0.19940 0.0063056 0.027545 +#> PEXPOSURE 0.08298 0.04821 0.0015246 0.007180 +#> PCTAGE65P 3.76613 0.59918 0.0189478 0.070409 +#> PCTOWNHOME -0.40409 0.19845 0.0062755 0.026634 #> #> 2. Quantiles for each variable: #> -#> 2.5% 25% 50% 75% 97.5% -#> lambda 0.015933 0.03450 0.04683 0.05826 0.07461 -#> (Intercept) -1.036934 -0.74613 -0.62918 -0.50290 -0.22351 -#> PEXPOSURE 0.006434 0.04402 0.07379 0.09781 0.14820 -#> PCTAGE65P 2.287060 3.30526 3.82691 4.21965 4.95288 -#> PCTOWNHOME -0.794631 -0.55436 -0.42774 -0.29149 -0.02868 +#> 2.5% 25% 50% 75% 97.5% +#> lambda 0.0141671 0.03110 0.04116 0.05429 0.07622 +#> (Intercept) -1.0578824 -0.77036 -0.61555 -0.49229 -0.28781 +#> PEXPOSURE 0.0007409 0.04776 0.07779 0.12221 0.17562 +#> PCTAGE65P 2.6714807 3.33132 3.75638 4.27260 4.92908 +#> PCTOWNHOME -0.8107680 -0.53926 -0.39020 -0.25915 -0.01036 #> -# \dontrun{ +if (FALSE) { esar0w <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, weights=POP8) summary(esar0) -#> -#> Call:errorsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, -#> data = nydata, listw = listw_NY) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.56754 -0.38239 -0.02643 0.33109 4.01219 -#> -#> Type: error -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.618193 0.176784 -3.4969 0.0004707 -#> PEXPOSURE 0.071014 0.042051 1.6888 0.0912635 -#> PCTAGE65P 3.754200 0.624722 6.0094 1.862e-09 -#> PCTOWNHOME -0.419890 0.191329 -2.1946 0.0281930 -#> -#> Lambda: 0.040487, LR test value: 5.2438, p-value: 0.022026 -#> Asymptotic standard error: 0.016214 -#> z-value: 2.4971, p-value: 0.01252 -#> Wald statistic: 6.2356, p-value: 0.01252 -#> -#> Log likelihood: -276.1069 for error model -#> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333) -#> Number of observations: 281 -#> Number of parameters estimated: 6 -#> AIC: 564.21, (AIC for lm: 567.46) -#> res <- MCMCsamp(esar0w, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> lambda 0.01264 0.01693 0.0002394 0.0009194 -#> (Intercept) -0.79572 0.14919 0.0021098 0.0085953 -#> PEXPOSURE 0.08218 0.02967 0.0004196 0.0017624 -#> PCTAGE65P 3.77097 0.56713 0.0080204 0.0321372 -#> PCTOWNHOME -0.37425 0.15969 0.0022584 0.0095249 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> lambda -0.01898 0.00148 0.01220 0.02415 0.04619 -#> (Intercept) -1.07804 -0.89673 -0.80093 -0.68897 -0.50028 -#> PEXPOSURE 0.02375 0.06145 0.08308 0.10270 0.14032 -#> PCTAGE65P 2.61747 3.40744 3.76515 4.13820 4.86343 -#> PCTOWNHOME -0.68124 -0.48248 -0.36610 -0.26440 -0.06796 -#> esar1 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, etype="emixed") summary(esar1) -#> -#> Call:errorsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, -#> data = nydata, listw = listw_NY, etype = "emixed") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.81562 -0.37641 -0.02224 0.33638 4.00054 -#> -#> Type: error -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -1.118019 0.247425 -4.5186 6.225e-06 -#> PEXPOSURE 0.218279 0.079245 2.7545 0.005879 -#> PCTAGE65P 3.416477 0.645587 5.2920 1.210e-07 -#> PCTOWNHOME 0.036593 0.249835 0.1465 0.883551 -#> lag.(Intercept) 0.121515 0.057636 2.1083 0.035003 -#> lag.PEXPOSURE -0.035075 0.015943 -2.2000 0.027808 -#> lag.PCTAGE65P 0.263096 0.220118 1.1953 0.231989 -#> lag.PCTOWNHOME -0.155680 0.059213 -2.6291 0.008560 -#> -#> Lambda: 0.022723, LR test value: 1.6846, p-value: 0.19432 -#> Asymptotic standard error: 0.017169 -#> z-value: 1.3235, p-value: 0.18567 -#> Wald statistic: 1.7516, p-value: 0.18567 -#> -#> Log likelihood: -269.5398 for error model -#> ML residual variance (sigma squared): 0.39759, (sigma: 0.63055) -#> Number of observations: 281 -#> Number of parameters estimated: 10 -#> AIC: 559.08, (AIC for lm: 558.76) -#> res <- MCMCsamp(esar1, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> lambda 0.02664 0.01731 0.0002448 0.001445 -#> (Intercept) -1.11779 0.24615 0.0034810 0.018477 -#> PEXPOSURE 0.21666 0.07903 0.0011177 0.006156 -#> PCTAGE65P 3.44403 0.66613 0.0094205 0.054200 -#> PCTOWNHOME 0.02791 0.24636 0.0034841 0.017966 -#> lag.(Intercept) 0.11899 0.05413 0.0007655 0.003895 -#> lag.PEXPOSURE -0.03408 0.01569 0.0002219 0.001188 -#> lag.PCTAGE65P 0.26541 0.22248 0.0031464 0.017271 -#> lag.PCTOWNHOME -0.15371 0.05829 0.0008243 0.004495 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> lambda -0.007983 0.01561 0.02677 0.03815 0.060478 -#> (Intercept) -1.630741 -1.26755 -1.12615 -0.95492 -0.658527 -#> PEXPOSURE 0.060952 0.16303 0.21795 0.27110 0.369123 -#> PCTAGE65P 2.179985 2.98375 3.46602 3.94588 4.721013 -#> PCTOWNHOME -0.454307 -0.13862 0.02277 0.17105 0.556899 -#> lag.(Intercept) 0.009429 0.08303 0.11691 0.15236 0.225551 -#> lag.PEXPOSURE -0.064503 -0.04472 -0.03393 -0.02320 -0.006341 -#> lag.PCTAGE65P -0.207992 0.12362 0.25575 0.43532 0.690657 -#> lag.PCTOWNHOME -0.276539 -0.19200 -0.15699 -0.11334 -0.041910 -#> lsar0 <- lagsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(lsar0) -#> -#> Call:lagsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.586752 -0.391580 -0.022469 0.338017 4.029430 -#> -#> Type: lag -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.514495 0.156154 -3.2948 0.000985 -#> PEXPOSURE 0.047627 0.034509 1.3801 0.167542 -#> PCTAGE65P 3.648198 0.599046 6.0900 1.129e-09 -#> PCTOWNHOME -0.414601 0.169554 -2.4453 0.014475 -#> -#> Rho: 0.038893, LR test value: 6.9683, p-value: 0.0082967 -#> Asymptotic standard error: 0.015053 -#> z-value: 2.5837, p-value: 0.0097755 -#> Wald statistic: 6.6754, p-value: 0.0097755 -#> -#> Log likelihood: -275.2447 for lag model -#> ML residual variance (sigma squared): 0.41166, (sigma: 0.6416) -#> Number of observations: 281 -#> Number of parameters estimated: 6 -#> AIC: 562.49, (AIC for lm: 567.46) -#> LM test for residual autocorrelation -#> test value: 1.4633, p-value: 0.22641 -#> res <- MCMCsamp(lsar0, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> rho 0.03941 0.01438 0.0002033 0.0008153 -#> (Intercept) -0.50459 0.16645 0.0023540 0.0103933 -#> PEXPOSURE 0.04431 0.03807 0.0005383 0.0025311 -#> PCTAGE65P 3.58732 0.63463 0.0089750 0.0395655 -#> PCTOWNHOME -0.40438 0.17782 0.0025147 0.0104785 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> rho 0.01139 0.02985 0.03940 0.04977 0.06749 -#> (Intercept) -0.83475 -0.61489 -0.50587 -0.39908 -0.16192 -#> PEXPOSURE -0.03039 0.01908 0.04369 0.06867 0.12181 -#> PCTAGE65P 2.40599 3.17100 3.60271 4.00568 4.84911 -#> PCTOWNHOME -0.75254 -0.51872 -0.40299 -0.28250 -0.06687 -#> lsar1 <- lagsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, type="mixed") summary(lsar1) -#> -#> Call:lagsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY, type = "mixed") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.799308 -0.390125 -0.021371 0.346128 3.965251 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -1.131233 0.249631 -4.5316 5.853e-06 -#> PEXPOSURE 0.218364 0.079301 2.7536 0.005894 -#> PCTAGE65P 3.361158 0.654123 5.1384 2.771e-07 -#> PCTOWNHOME 0.071903 0.253967 0.2831 0.777085 -#> lag.(Intercept) 0.132544 0.056175 2.3595 0.018300 -#> lag.PEXPOSURE -0.035239 0.015536 -2.2681 0.023322 -#> lag.PCTAGE65P 0.161685 0.223690 0.7228 0.469798 -#> lag.PCTOWNHOME -0.140681 0.058529 -2.4036 0.016234 -#> -#> Rho: 0.026981, LR test value: 2.558, p-value: 0.10974 -#> Asymptotic standard error: 0.016766 -#> z-value: 1.6093, p-value: 0.10755 -#> Wald statistic: 2.5899, p-value: 0.10755 -#> -#> Log likelihood: -269.1031 for mixed model -#> ML residual variance (sigma squared): 0.39587, (sigma: 0.62918) -#> Number of observations: 281 -#> Number of parameters estimated: 10 -#> AIC: 558.21, (AIC for lm: 558.76) -#> LM test for residual autocorrelation -#> test value: 4.908, p-value: 0.026732 -#> res <- MCMCsamp(lsar1, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> rho 0.02639 0.01731 0.0002447 0.001376 -#> (Intercept) -1.09340 0.24932 0.0035259 0.018650 -#> PEXPOSURE 0.21187 0.08271 0.0011697 0.006697 -#> PCTAGE65P 3.33757 0.62894 0.0088945 0.046186 -#> PCTOWNHOME 0.02806 0.25799 0.0036485 0.019903 -#> lag.(Intercept) 0.12638 0.05644 0.0007982 0.004347 -#> lag.PEXPOSURE -0.03369 0.01644 0.0002324 0.001375 -#> lag.PCTAGE65P 0.16969 0.23805 0.0033665 0.018638 -#> lag.PCTOWNHOME -0.13408 0.06147 0.0008693 0.004898 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> rho -0.009348 0.014863 0.02659 0.03817 0.059014 -#> (Intercept) -1.566398 -1.247599 -1.11066 -0.93527 -0.592916 -#> PEXPOSURE 0.051143 0.154615 0.21446 0.26488 0.379335 -#> PCTAGE65P 2.070955 2.946428 3.33766 3.75663 4.666186 -#> PCTOWNHOME -0.452715 -0.146853 0.02290 0.20161 0.563357 -#> lag.(Intercept) 0.018492 0.088830 0.12357 0.16471 0.232746 -#> lag.PEXPOSURE -0.066219 -0.044733 -0.03338 -0.02204 -0.003191 -#> lag.PCTAGE65P -0.341037 0.001146 0.17672 0.35102 0.622026 -#> lag.PCTOWNHOME -0.254849 -0.177764 -0.13151 -0.08932 -0.020180 -#> ssar0 <- sacsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(ssar0) -#> -#> Call:sacsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.468382 -0.375687 -0.034996 0.314714 3.833950 -#> -#> Type: sac -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -0.386572 0.123188 -3.1381 0.001701 -#> PEXPOSURE 0.026684 0.024013 1.1112 0.266479 -#> PCTAGE65P 3.089824 0.562851 5.4896 4.029e-08 -#> PCTOWNHOME -0.323052 0.137449 -2.3503 0.018756 -#> -#> Rho: 0.089451 -#> Asymptotic standard error: 0.019427 -#> z-value: 4.6046, p-value: 4.1325e-06 -#> Lambda: -0.08192 -#> Asymptotic standard error: 0.033201 -#> z-value: -2.4674, p-value: 0.01361 -#> -#> LR test value: 10.114, p-value: 0.0063661 -#> -#> Log likelihood: -273.672 for sac model -#> ML residual variance (sigma squared): 0.3766, (sigma: 0.61368) -#> Number of observations: 281 -#> Number of parameters estimated: 7 -#> AIC: 561.34, (AIC for lm: 567.46) -#> res <- MCMCsamp(ssar0, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> rho -0.007526 0.08591 0.001215 0.02618 -#> lambda 0.028126 0.08464 0.001197 0.02421 -#> (Intercept) -0.632441 0.27560 0.003898 0.04505 -#> PEXPOSURE 0.089391 0.07677 0.001086 0.01375 -#> PCTAGE65P 3.296382 0.68594 0.009701 0.05142 -#> PCTOWNHOME -0.291066 0.21055 0.002978 0.02161 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> rho -0.15552 -0.09457 0.008668 0.07076 0.1114 -#> lambda -0.12817 -0.04749 0.035462 0.11419 0.1398 -#> (Intercept) -1.21573 -0.81705 -0.585620 -0.41254 -0.2203 -#> PEXPOSURE -0.01643 0.02928 0.070437 0.13724 0.2618 -#> PCTAGE65P 2.06841 2.82620 3.299895 3.72438 4.6615 -#> PCTOWNHOME -0.67524 -0.43644 -0.296051 -0.15542 0.1733 -#> ssar1 <- sacsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, type="sacmixed") summary(ssar1) -#> -#> Call:sacsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata, -#> listw = listw_NY, type = "sacmixed") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -1.633958 -0.363826 -0.019927 0.348238 3.655509 -#> -#> Type: sacmixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -1.133298 0.247495 -4.5791 4.670e-06 -#> PEXPOSURE 0.206963 0.074480 2.7788 0.005456 -#> PCTAGE65P 3.083983 0.671081 4.5955 4.316e-06 -#> PCTOWNHOME 0.174800 0.256280 0.6821 0.495196 -#> lag.(Intercept) 0.153427 0.050817 3.0192 0.002534 -#> lag.PEXPOSURE -0.033400 0.013817 -2.4173 0.015634 -#> lag.PCTAGE65P -0.079738 0.222144 -0.3589 0.719634 -#> lag.PCTOWNHOME -0.102502 0.056760 -1.8059 0.070940 -#> -#> Rho: 0.092495 -#> Asymptotic standard error: 0.023829 -#> z-value: 3.8817, p-value: 0.00010375 -#> Lambda: -0.091069 -#> Asymptotic standard error: 0.038431 -#> z-value: -2.3697, p-value: 0.017804 -#> -#> LR test value: 22.379, p-value: 0.0010335 -#> -#> Log likelihood: -267.5392 for sacmixed model -#> ML residual variance (sigma squared): 0.35617, (sigma: 0.5968) -#> Number of observations: 281 -#> Number of parameters estimated: 11 -#> AIC: 557.08, (AIC for lm: 567.46) -#> res <- MCMCsamp(ssar1, mcmc=5000, burnin=500, listw=listw_NY) summary(res) -#> -#> Iterations = 1:5000 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 5000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> rho -0.006243 0.07349 0.0010394 0.016978 -#> lambda 0.023036 0.07006 0.0009907 0.014064 -#> (Intercept) -1.074786 0.25700 0.0036346 0.020044 -#> PEXPOSURE 0.213653 0.07987 0.0011296 0.006275 -#> PCTAGE65P 3.369138 0.64267 0.0090887 0.048578 -#> PCTOWNHOME 0.035386 0.26783 0.0037877 0.021115 -#> lag.(Intercept) 0.106733 0.06761 0.0009562 0.006909 -#> lag.PEXPOSURE -0.033731 0.01826 0.0002582 0.001782 -#> lag.PCTAGE65P 0.239723 0.32512 0.0045978 0.037622 -#> lag.PCTOWNHOME -0.138544 0.07558 0.0010688 0.007509 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> rho -0.15259 -0.058974 -0.001855 0.05403 0.110005 -#> lambda -0.12122 -0.034805 0.030348 0.07995 0.127609 -#> (Intercept) -1.63521 -1.245895 -1.066876 -0.92132 -0.594111 -#> PEXPOSURE 0.05949 0.156813 0.216836 0.26607 0.373617 -#> PCTAGE65P 2.10893 2.932956 3.351803 3.78560 4.733524 -#> PCTOWNHOME -0.48089 -0.148442 0.034779 0.22580 0.583465 -#> lag.(Intercept) -0.03405 0.065889 0.107361 0.15007 0.231268 -#> lag.PEXPOSURE -0.07338 -0.045326 -0.034166 -0.02251 0.004342 -#> lag.PCTAGE65P -0.39085 -0.002034 0.212151 0.47009 0.871781 -#> lag.PCTOWNHOME -0.29853 -0.183541 -0.139221 -0.08629 0.008422 -#> -# } +}Examples
Examples
set.seed(123) system.time(MEbinom1 <- ME(c(hopkins_part) ~ 1, family="binomial", listw=lw, alpha=0.05, verbose=TRUE, nsim=49)) -#> eV[,2], I: 0.07978424 ZI: NA, pr(ZI): 0.04 -#> eV[,54], I: 0.05704227 ZI: NA, pr(ZI): 0.12 +#> eV[,1], I: 0.08290518 ZI: NA, pr(ZI): 0.04 +#> eV[,9], I: 0.06426565 ZI: NA, pr(ZI): 0.14 #> user system elapsed -#> 1.295 0.003 1.304 +#> 1.335 0.006 1.352 glmME <- glm(c(hopkins_part) ~ 1 + fitted(MEbinom1), family="binomial") #anova(glmME, test="Chisq") coef(summary(glmME)) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) -1.194892 0.1600429 -7.466073 8.262357e-14 -#> fitted(MEbinom1)vec2 9.030818 2.5381594 3.558019 3.736629e-04 -#> fitted(MEbinom1)vec54 -8.759610 2.5205844 -3.475230 5.104161e-04 +#> Estimate Std. Error z value Pr(>|z|) +#> (Intercept) -1.146132 0.1542253 -7.431543 1.073378e-13 +#> fitted(MEbinom1)vec1 8.293309 2.4532307 3.380566 7.233663e-04 +#> fitted(MEbinom1)vec9 5.215112 2.3949596 2.177537 2.944054e-02 anova(glmbase, glmME, test="Chisq") #> Analysis of Deviance Table #> @@ -183,10 +183,10 @@Examples
#> Model 2: c(hopkins_part) ~ 1 + fitted(MEbinom1) #> Resid. Df Resid. Dev Df Deviance Pr(>Chi) #> 1 255 292.23 -#> 2 253 266.27 2 25.956 2.311e-06 *** +#> 2 253 275.39 2 16.841 0.0002203 *** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -# \dontrun{ +if (FALSE) { require("sf", quietly=TRUE) columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE) #require("spdep", quietly=TRUE) @@ -195,108 +195,30 @@Examples
lmbase <- lm(CRIME ~ INC + HOVAL, data=columbus) lagcol <- SpatialFiltering(CRIME ~ 1, ~ INC + HOVAL, data=columbus, nb=col.gal.nb, style="W", alpha=0.1, verbose=TRUE) -#> Step 0 SelEvec 0 MinMi 0.2123742 ZMinMi 2.681 Pr(ZI) 0.007340246 -#> Step 1 SelEvec 6 MinMi 0.1178225 ZMinMi 1.84512 Pr(ZI) 0.06502014 -#> Step 2 SelEvec 4 MinMi 0.06242664 ZMinMi 1.494821 Pr(ZI) 0.1349611 lagcol -#> Step SelEvec Eval MinMi ZMinMi Pr(ZI) R2 gamma -#> 0 0 0 0.0000000 0.21237415 2.681000 0.007340246 0.5524040 0.00000 -#> 1 1 6 0.7161123 0.11782248 1.845120 0.065020139 0.6038801 25.46181 -#> 2 2 4 0.8682938 0.06242664 1.494821 0.134961136 0.6531288 26.68319 lmlag <- lm(CRIME ~ INC + HOVAL + fitted(lagcol), data=columbus) anova(lmbase, lmlag) -#> Analysis of Variance Table -#> -#> Model 1: CRIME ~ INC + HOVAL -#> Model 2: CRIME ~ INC + HOVAL + fitted(lagcol) -#> Res.Df RSS Df Sum of Sq F Pr(>F) -#> 1 46 6014.9 -#> 2 44 4661.3 2 1353.6 6.3884 0.003666 ** -#> --- -#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 set.seed(123) system.time(lagcol1 <- ME(CRIME ~ INC + HOVAL, data=columbus, family="gaussian", listw=lw, alpha=0.1, verbose=TRUE)) -#> eV[,6], I: 0.1178225 ZI: NA, pr(ZI): 0.08 -#> eV[,4], I: 0.06242664 ZI: NA, pr(ZI): 0.27 -#> user system elapsed -#> 0.783 0.000 0.786 lagcol1 -#> Eigenvector ZI pr(ZI) -#> 0 NA NA 0.01 -#> 1 6 NA 0.08 -#> 2 4 NA 0.27 lmlag1 <- lm(CRIME ~ INC + HOVAL + fitted(lagcol1), data=columbus) anova(lmbase, lmlag1) -#> Analysis of Variance Table -#> -#> Model 1: CRIME ~ INC + HOVAL -#> Model 2: CRIME ~ INC + HOVAL + fitted(lagcol1) -#> Res.Df RSS Df Sum of Sq F Pr(>F) -#> 1 46 6014.9 -#> 2 44 4661.3 2 1353.6 6.3884 0.003666 ** -#> --- -#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 set.seed(123) lagcol2 <- ME(CRIME ~ INC + HOVAL, data=columbus, family="gaussian", listw=lw, alpha=0.1, stdev=TRUE, verbose=TRUE) -#> eV[,6], I: 0.1178225 ZI: 1.5509, pr(ZI): 0.06046283 -#> eV[,4], I: 0.06242664 ZI: 0.681174, pr(ZI): 0.2478807 lagcol2 -#> Eigenvector ZI pr(ZI) -#> 0 NA 2.351591 0.009346653 -#> 1 6 1.550900 0.060462832 -#> 2 4 0.681174 0.247880696 lmlag2 <- lm(CRIME ~ INC + HOVAL + fitted(lagcol2), data=columbus) anova(lmbase, lmlag2) -#> Analysis of Variance Table -#> -#> Model 1: CRIME ~ INC + HOVAL -#> Model 2: CRIME ~ INC + HOVAL + fitted(lagcol2) -#> Res.Df RSS Df Sum of Sq F Pr(>F) -#> 1 46 6014.9 -#> 2 44 4661.3 2 1353.6 6.3884 0.003666 ** -#> --- -#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 NA.columbus <- columbus NA.columbus$CRIME[20:25] <- NA COL.ME.NA <- ME(CRIME ~ INC + HOVAL, data=NA.columbus, family="gaussian", listw=lw, alpha=0.1, stdev=TRUE, verbose=TRUE, na.action=na.exclude) -#> eV[,8], I: 0.1426723 ZI: 1.483169, pr(ZI): 0.06901474 -#> eV[,1], I: 0.09838877 ZI: 0.9862904, pr(ZI): 0.1619953 COL.ME.NA$na.action -#> 20 21 22 23 24 25 -#> 20 21 22 23 24 25 -#> attr(,"class") -#> [1] "exclude" summary(lm(CRIME ~ INC + HOVAL + fitted(COL.ME.NA), data=NA.columbus, na.action=na.exclude)) -#> -#> Call: -#> lm(formula = CRIME ~ INC + HOVAL + fitted(COL.ME.NA), data = NA.columbus, -#> na.action = na.exclude) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -30.1382 -6.0105 0.4095 7.1504 19.9399 -#> -#> Coefficients: -#> Estimate Std. Error t value Pr(>|t|) -#> (Intercept) 66.92248 5.28663 12.659 3.33e-15 *** -#> INC -1.40484 0.35678 -3.938 0.00034 *** -#> HOVAL -0.30446 0.09831 -3.097 0.00366 ** -#> fitted(COL.ME.NA)1 29.69422 10.58481 2.805 0.00788 ** -#> fitted(COL.ME.NA)2 26.61612 11.29187 2.357 0.02367 * -#> --- -#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -#> -#> Residual standard error: 10.48 on 38 degrees of freedom -#> (6 observations deleted due to missingness) -#> Multiple R-squared: 0.6294, Adjusted R-squared: 0.5904 -#> F-statistic: 16.13 on 4 and 38 DF, p-value: 8.353e-08 -#> nc.sids <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE) rn <- as.character(nc.sids$FIPS) ncCC89_nb <- spdep::read.gal(system.file("weights/ncCC89.gal", package="spData")[1], @@ -308,43 +230,12 @@Examples
set.seed(123) MEpois1 <- ME(SID74 ~ 1, data=nc.sids, offset=log(BIR74), family="poisson", listw=spdep::nb2listw(ncCR85_nb, style="B"), alpha=0.2, verbose=TRUE) -#> eV[,1], I: 0.1327384 ZI: NA, pr(ZI): 0.03 -#> eV[,8], I: 0.06936385 ZI: NA, pr(ZI): 0.12 -#> eV[,4], I: 0.03584503 ZI: NA, pr(ZI): 0.3 MEpois1 -#> Eigenvector ZI pr(ZI) -#> 0 NA NA 0.01 -#> 1 1 NA 0.03 -#> 2 8 NA 0.12 -#> 3 4 NA 0.30 glmME <- glm(SID74 ~ 1 + fitted(MEpois1), data=nc.sids, offset=log(BIR74), family="poisson") anova(glmME, test="Chisq") -#> Analysis of Deviance Table -#> -#> Model: poisson, link: log -#> -#> Response: SID74 -#> -#> Terms added sequentially (first to last) -#> -#> -#> Df Deviance Resid. Df Resid. Dev Pr(>Chi) -#> NULL 99 203.34 -#> fitted(MEpois1) 3 32.499 96 170.84 4.108e-07 *** -#> --- -#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 anova(glmbase, glmME, test="Chisq") -#> Analysis of Deviance Table -#> -#> Model 1: SID74 ~ 1 -#> Model 2: SID74 ~ 1 + fitted(MEpois1) -#> Resid. Df Resid. Dev Df Deviance Pr(>Chi) -#> 1 99 203.34 -#> 2 96 170.84 3 32.499 4.108e-07 *** -#> --- -#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -# } +}Examples
Examples
#> rho: 0.4310232 function value: -182.3904 #> rho: 0.4310232 function value: -182.3904 #> rho: 0.4310232 function value: -182.3904 +#> rho: 0.4310232 function value: -182.3904 (x <- summary(COL.lag.eig, correlation=TRUE)) #> #> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, @@ -398,7 +399,7 @@Examples
#> Type: lag #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.079251 7.177347 6.2808 3.369e-10 +#> (Intercept) 45.079249 7.177346 6.2808 3.369e-10 #> INC -1.031616 0.305143 -3.3808 0.0007229 #> HOVAL -0.265926 0.088499 -3.0049 0.0026570 #> @@ -411,9 +412,9 @@Examples
#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) #> Number of observations: 49 #> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) +#> AIC: NA (not available for weighted model), (AIC for lm: 382.75) #> LM test for residual autocorrelation -#> test value: 0.31955, p-value: 0.57188 +#> test value: 0.31954, p-value: 0.57188 #> #> Correlation of coefficients #> sigma rho (Intercept) INC @@ -424,150 +425,42 @@Examples
#> coef(x) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.0792505 7.17734654 6.280768 3.369041e-10 -#> INC -1.0316157 0.30514297 -3.380762 7.228517e-04 +#> (Intercept) 45.0792493 7.17734647 6.280768 3.369043e-10 +#> INC -1.0316157 0.30514297 -3.380762 7.228519e-04 #> HOVAL -0.2659263 0.08849862 -3.004863 2.657002e-03 -# \dontrun{ +if (FALSE) { COL.lag.eig$fdHess -#> [1] FALSE COL.lag.eig$resvar -#> sigma rho (Intercept) INC HOVAL -#> sigma 379.77510361 -0.3236306318 16.3015079 -0.29590801 -0.0202478459 -#> rho -0.32363063 0.0138487533 -0.6975717 0.01266245 0.0008664428 -#> (Intercept) 16.30150788 -0.6975716738 51.5143034 -1.32602704 -0.1616379856 -#> INC -0.29590801 0.0126624511 -1.3260270 0.09311223 -0.0117959715 -#> HOVAL -0.02024785 0.0008664428 -0.1616380 -0.01179597 0.0078320057 # using the apparent sign in Ord (1975, equation B.1) COL.lag.eigb <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw=listw, method="eigen", control=list(pre_eig=ev, OrdVsign=-1)) summary(COL.lag.eigb) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> method = "eigen", control = list(pre_eig = ev, OrdVsign = -1)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.68585 -5.35636 0.05421 6.02013 23.20555 -#> -#> Type: lag -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.079251 9.617835 4.6870 2.772e-06 -#> INC -1.031616 0.326524 -3.1594 0.001581 -#> HOVAL -0.265926 0.088855 -2.9928 0.002764 -#> -#> Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588 -#> Asymptotic standard error: 0.17322 -#> z-value: 2.4884, p-value: 0.012833 -#> Wald statistic: 6.1919, p-value: 0.012833 -#> -#> Log likelihood: -182.3904 for lag model -#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) -#> LM test for residual autocorrelation -#> test value: -0.93825, p-value: 1 -#> COL.lag.eigb$fdHess -#> [1] FALSE COL.lag.eigb$resvar -#> sigma rho (Intercept) INC HOVAL -#> sigma 388.59742931 -0.701154266 35.3176451 -0.64109248 -0.043867490 -#> rho -0.70115427 0.030003688 -1.5113074 0.02743353 0.001877171 -#> (Intercept) 35.31764510 -1.511307358 92.5027556 -2.07005707 -0.212549096 -#> INC -0.64109248 0.027433533 -2.0700571 0.10661800 -0.010871824 -#> HOVAL -0.04386749 0.001877171 -0.2125491 -0.01087182 0.007895242 # force numerical Hessian COL.lag.eig1 <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw=listw, method="Matrix", control=list(small=25)) summary(COL.lag.eig1) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> method = "Matrix", control = list(small = 25)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.68585 -5.35636 0.05421 6.02013 23.20555 -#> -#> Type: lag -#> Coefficients: (numerical Hessian approximate standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.07925 7.87089 5.7273 1.02e-08 -#> INC -1.03162 0.32842 -3.1411 0.001683 -#> HOVAL -0.26593 0.08823 -3.0140 0.002578 -#> -#> Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588 -#> Approximate (numerical Hessian) standard error: 0.12362 -#> z-value: 3.4868, p-value: 0.00048892 -#> Wald statistic: 12.157, p-value: 0.00048892 -#> -#> Log likelihood: -182.3904 for lag model -#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) -#> COL.lag.eig1$fdHess -#> rho (Intercept) INC HOVAL -#> rho 0.0152812237 -0.8345267 0.02005558 0.0002830225 -#> (Intercept) -0.8345266779 61.9509427 -1.78340176 -0.1334585288 -#> INC 0.0200555771 -1.7834018 0.10786140 -0.0122204417 -#> HOVAL 0.0002830225 -0.1334585 -0.01222044 0.0077846055 # force LeSage & Pace (2008, p. 57) approximation COL.lag.eig1a <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw=listw, method="Matrix", control=list(small=25), trs=trMatc) summary(COL.lag.eig1a) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> method = "Matrix", trs = trMatc, control = list(small = 25)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.68585 -5.35636 0.05421 6.02013 23.20555 -#> -#> Type: lag -#> Coefficients: (numerical Hessian approximate standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.079249 7.937572 5.6792 1.353e-08 -#> INC -1.031616 0.329337 -3.1324 0.001734 -#> HOVAL -0.265926 0.088222 -3.0143 0.002576 -#> -#> Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588 -#> Approximate (numerical Hessian) standard error: 0.12503 -#> z-value: 3.4473, p-value: 0.00056624 -#> Wald statistic: 11.884, p-value: 0.00056624 -#> -#> Log likelihood: -182.3904 for lag model -#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) -#> COL.lag.eig1a$fdHess -#> sigma2 rho (Intercept) INC HOVAL -#> sigma2 380.749548407 -0.3653290756 19.9519208 -0.47947507 -0.0067851124 -#> rho -0.365329076 0.0156331058 -0.8537795 0.02051762 0.0002903475 -#> (Intercept) 19.951920812 -0.8537795385 63.0050547 -1.80875071 -0.1338515483 -#> INC -0.479475072 0.0205176238 -1.8087507 0.10846276 -0.0122071279 -#> HOVAL -0.006785112 0.0002903475 -0.1338515 -0.01220713 0.0077831895 COL.lag.eig$resvar[2,2] -#> [1] 0.01384875 # using the apparent sign in Ord (1975, equation B.1) COL.lag.eigb$resvar[2,2] -#> [1] 0.03000369 # force numerical Hessian COL.lag.eig1$fdHess[1,1] -#> [1] 0.01528122 # force LeSage & Pace (2008, p. 57) approximation COL.lag.eig1a$fdHess[2,2] -#> [1] 0.01563311 -# } +} system.time(COL.lag.M <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix", quiet=FALSE)) #> #> Spatial lag model #> Jacobian calculated using sparse matrix Cholesky decomposition +#> Warning: the default value of argument 'sqrt' of method 'determinant(<CHMfactor>, <logical>)' may change from TRUE to FALSE as soon as the next release of Matrix; set 'sqrt' when programming #> rho: -0.2364499 function value: -192.9523 #> rho: 0.2354499 function value: -183.542 #> rho: 0.5271001 function value: -182.7039 @@ -583,7 +476,7 @@Examples
#> Computing eigenvalues ... #> #> user system elapsed -#> 0.196 0.000 0.197 +#> 0.145 0.001 0.146 summary(COL.lag.M) #> #> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, @@ -609,7 +502,7 @@Examples
#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) #> Number of observations: 49 #> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) +#> AIC: NA (not available for weighted model), (AIC for lm: 382.75) #> LM test for residual autocorrelation #> test value: 0.31954, p-value: 0.57188 #> @@ -618,455 +511,66 @@Examples
#> Direct Indirect Total #> INC -1.0860220 -0.7270848 -1.8131068 #> HOVAL -0.2799509 -0.1874254 -0.4673763 -# \dontrun{ +if (FALSE) { system.time(COL.lag.sp <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw=listw, method="spam", quiet=FALSE)) -#> -#> Spatial lag model -#> Jacobian calculated using sparse matrix Cholesky decomposition -#> rho: -0.2364499 function value: -192.9523 -#> rho: 0.2354499 function value: -183.542 -#> rho: 0.5271001 function value: -182.7039 -#> rho: 0.4455543 function value: -182.3974 -#> rho: 0.4267907 function value: -182.391 -#> rho: 0.4311986 function value: -182.3904 -#> rho: 0.4310114 function value: -182.3904 -#> rho: 0.4310231 function value: -182.3904 -#> rho: 0.4310232 function value: -182.3904 -#> rho: 0.4310232 function value: -182.3904 -#> rho: 0.4310232 function value: -182.3904 -#> rho: 0.4310232 function value: -182.3904 -#> Computing eigenvalues ... -#> -#> user system elapsed -#> 0.518 0.001 0.522 summary(COL.lag.sp) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> method = "spam", quiet = FALSE) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.68585 -5.35636 0.05421 6.02013 23.20555 -#> -#> Type: lag -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.079250 7.177347 6.2808 3.369e-10 -#> INC -1.031616 0.305143 -3.3808 0.0007229 -#> HOVAL -0.265926 0.088499 -3.0049 0.0026570 -#> -#> Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588 -#> Asymptotic standard error: 0.11768 -#> z-value: 3.6626, p-value: 0.00024962 -#> Wald statistic: 13.415, p-value: 0.00024962 -#> -#> Log likelihood: -182.3904 for lag model -#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) -#> LM test for residual autocorrelation -#> test value: 0.31954, p-value: 0.57188 -#> COL.lag.B <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, spdep::nb2listw(COL.nb, style="B"), control=list(pre_eig=ev)) summary(COL.lag.B) -#> -#> Call: -#> lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = spdep::nb2listw(COL.nb, -#> style = "B"), control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -33.6620 -4.8615 -1.3576 5.1567 25.7563 -#> -#> Type: lag -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 51.604815 6.075285 8.4942 < 2.2e-16 -#> INC -1.154463 0.301808 -3.8252 0.0001307 -#> HOVAL -0.251633 0.087612 -2.8721 0.0040773 -#> -#> Rho: 0.054543, LR test value: 13.453, p-value: 0.00024461 -#> Asymptotic standard error: 0.014836 -#> z-value: 3.6763, p-value: 0.00023662 -#> Wald statistic: 13.515, p-value: 0.00023662 -#> -#> Log likelihood: -180.6507 for lag model -#> ML residual variance (sigma squared): 93.22, (sigma: 9.655) -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 371.3, (AIC for lm: 382.75) -#> LM test for residual autocorrelation -#> test value: 0.006827, p-value: 0.93415 -#> COL.mixed.B <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, spdep::nb2listw(COL.nb, style="B"), type="mixed", tol.solve=1e-9, control=list(pre_eig=ev)) summary(COL.mixed.B) -#> -#> Call: -#> lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = spdep::nb2listw(COL.nb, -#> style = "B"), type = "mixed", tol.solve = 1e-09, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -34.8460 -4.2057 -0.1195 4.6525 21.6112 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 5.4215e+01 5.6639e+00 9.5719 < 2.2e-16 -#> INC -8.2386e-01 3.1643e-01 -2.6036 0.0092262 -#> HOVAL -3.0085e-01 8.5629e-02 -3.5134 0.0004424 -#> lag.(Intercept) -7.5493e+00 1.7307e+00 -4.3620 1.289e-05 -#> lag.INC 2.1544e-05 1.2216e-01 0.0002 0.9998593 -#> lag.HOVAL 7.2458e-02 3.9007e-02 1.8576 0.0632281 -#> -#> Rho: 0.15212, LR test value: 7.435, p-value: 0.0063967 -#> Asymptotic standard error: 0.015565 -#> z-value: 9.7735, p-value: < 2.22e-16 -#> Wald statistic: 95.522, p-value: < 2.22e-16 -#> -#> Log likelihood: -177.7722 for mixed model -#> ML residual variance (sigma squared): 82.502, (sigma: 9.0831) -#> Number of observations: 49 -#> Number of parameters estimated: 8 -#> AIC: 371.54, (AIC for lm: 376.98) -#> LM test for residual autocorrelation -#> test value: -26.79, p-value: 1 -#> COL.mixed.W <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, type="mixed", control=list(pre_eig=ev)) summary(COL.mixed.W) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> type = "mixed", control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.47829 -6.46731 -0.33835 6.05200 22.62969 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 42.822414 12.667204 3.3806 0.0007233 -#> INC -0.914223 0.331094 -2.7612 0.0057586 -#> HOVAL -0.293738 0.089212 -3.2926 0.0009927 -#> lag.INC -0.520284 0.565129 -0.9206 0.3572355 -#> lag.HOVAL 0.245640 0.178917 1.3729 0.1697756 -#> -#> Rho: 0.42634, LR test value: 5.3693, p-value: 0.020494 -#> Asymptotic standard error: 0.15623 -#> z-value: 2.7288, p-value: 0.0063561 -#> Wald statistic: 7.4465, p-value: 0.0063561 -#> -#> Log likelihood: -181.3935 for mixed model -#> ML residual variance (sigma squared): 91.791, (sigma: 9.5808) -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 376.79, (AIC for lm: 380.16) -#> LM test for residual autocorrelation -#> test value: 0.28919, p-value: 0.59074 -#> COL.mixed.D00 <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin=TRUE, control=list(pre_eig=ev)) summary(COL.mixed.D00) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> Durbin = TRUE, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.47829 -6.46731 -0.33835 6.05200 22.62969 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 42.822414 12.667204 3.3806 0.0007233 -#> INC -0.914223 0.331094 -2.7612 0.0057586 -#> HOVAL -0.293738 0.089212 -3.2926 0.0009927 -#> lag.INC -0.520284 0.565129 -0.9206 0.3572355 -#> lag.HOVAL 0.245640 0.178917 1.3729 0.1697756 -#> -#> Rho: 0.42634, LR test value: 5.3693, p-value: 0.020494 -#> Asymptotic standard error: 0.15623 -#> z-value: 2.7288, p-value: 0.0063561 -#> Wald statistic: 7.4465, p-value: 0.0063561 -#> -#> Log likelihood: -181.3935 for mixed model -#> ML residual variance (sigma squared): 91.791, (sigma: 9.5808) -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 376.79, (AIC for lm: 380.16) -#> LM test for residual autocorrelation -#> test value: 0.28919, p-value: 0.59074 -#> COL.mixed.D01 <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin=FALSE, control=list(pre_eig=ev)) summary(COL.mixed.D01) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> Durbin = FALSE, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.68585 -5.35636 0.05421 6.02013 23.20555 -#> -#> Type: lag -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 45.079251 7.177347 6.2808 3.369e-10 -#> INC -1.031616 0.305143 -3.3808 0.0007229 -#> HOVAL -0.265926 0.088499 -3.0049 0.0026570 -#> -#> Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588 -#> Asymptotic standard error: 0.11768 -#> z-value: 3.6626, p-value: 0.00024962 -#> Wald statistic: 13.415, p-value: 0.00024962 -#> -#> Log likelihood: -182.3904 for lag model -#> ML residual variance (sigma squared): 95.494, (sigma: 9.7721) -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 374.78, (AIC for lm: 382.75) -#> LM test for residual autocorrelation -#> test value: 0.31955, p-value: 0.57188 -#> COL.mixed.D1 <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin= ~ INC + HOVAL, control=list(pre_eig=ev)) summary(COL.mixed.D1) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> Durbin = ~INC + HOVAL, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.47829 -6.46731 -0.33835 6.05200 22.62969 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 42.822414 12.667204 3.3806 0.0007233 -#> INC -0.914223 0.331094 -2.7612 0.0057586 -#> HOVAL -0.293738 0.089212 -3.2926 0.0009927 -#> lag.INC -0.520284 0.565129 -0.9206 0.3572355 -#> lag.HOVAL 0.245640 0.178917 1.3729 0.1697756 -#> -#> Rho: 0.42634, LR test value: 5.3693, p-value: 0.020494 -#> Asymptotic standard error: 0.15623 -#> z-value: 2.7288, p-value: 0.0063561 -#> Wald statistic: 7.4465, p-value: 0.0063561 -#> -#> Log likelihood: -181.3935 for mixed model -#> ML residual variance (sigma squared): 91.791, (sigma: 9.5808) -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 376.79, (AIC for lm: 380.16) -#> LM test for residual autocorrelation -#> test value: 0.28919, p-value: 0.59074 -#> f <- CRIME ~ INC + HOVAL COL.mixed.D2 <- lagsarlm(f, data=COL.OLD, listw, Durbin=as.formula(delete.response(terms(f))), control=list(pre_eig=ev)) summary(COL.mixed.D2) -#> -#> Call: -#> lagsarlm(formula = f, data = COL.OLD, listw = listw, Durbin = as.formula(delete.response(terms(f))), -#> control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.47829 -6.46731 -0.33835 6.05200 22.62969 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 42.822414 12.667204 3.3806 0.0007233 -#> INC -0.914223 0.331094 -2.7612 0.0057586 -#> HOVAL -0.293738 0.089212 -3.2926 0.0009927 -#> lag.INC -0.520284 0.565129 -0.9206 0.3572355 -#> lag.HOVAL 0.245640 0.178917 1.3729 0.1697756 -#> -#> Rho: 0.42634, LR test value: 5.3693, p-value: 0.020494 -#> Asymptotic standard error: 0.15623 -#> z-value: 2.7288, p-value: 0.0063561 -#> Wald statistic: 7.4465, p-value: 0.0063561 -#> -#> Log likelihood: -181.3935 for mixed model -#> ML residual variance (sigma squared): 91.791, (sigma: 9.5808) -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 376.79, (AIC for lm: 380.16) -#> LM test for residual autocorrelation -#> test value: 0.28919, p-value: 0.59074 -#> COL.mixed.D1a <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin= ~ INC, control=list(pre_eig=ev)) summary(COL.mixed.D1a) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> Durbin = ~INC, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.82800 -5.85207 0.12047 6.00137 23.19963 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 48.814687 12.198232 4.0018 6.287e-05 -#> INC -1.006620 0.330641 -3.0445 0.002331 -#> HOVAL -0.265514 0.088768 -2.9911 0.002780 -#> lag.INC -0.186684 0.530550 -0.3519 0.724936 -#> -#> Rho: 0.39229, LR test value: 4.5007, p-value: 0.033881 -#> Asymptotic standard error: 0.1561 -#> z-value: 2.513, p-value: 0.011971 -#> Wald statistic: 6.3151, p-value: 0.011971 -#> -#> Log likelihood: -182.3328 for mixed model -#> ML residual variance (sigma squared): 96.122, (sigma: 9.8042) -#> Number of observations: 49 -#> Number of parameters estimated: 6 -#> AIC: 376.67, (AIC for lm: 379.17) -#> LM test for residual autocorrelation -#> test value: 2.6134, p-value: 0.10596 -#> try(COL.mixed.D1 <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin= ~ inc + HOVAL, control=list(pre_eig=ev))) -#> Error in eval(predvars, data, env) : object 'inc' not found try(COL.mixed.D1 <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin= ~ DISCBD + HOVAL, control=list(pre_eig=ev))) -#> Error in lagsarlm(CRIME ~ INC + HOVAL, data = COL.OLD, listw, Durbin = ~DISCBD + : -#> WX variables not in X: DISCBD NA.COL.OLD <- COL.OLD NA.COL.OLD$CRIME[20:25] <- NA COL.lag.NA <- lagsarlm(CRIME ~ INC + HOVAL, data=NA.COL.OLD, listw, na.action=na.exclude) COL.lag.NA$na.action -#> 1020 1021 1022 1023 1024 1025 -#> 20 21 22 23 24 25 -#> attr(,"class") -#> [1] "exclude" COL.lag.NA -#> -#> Call: -#> lagsarlm(formula = CRIME ~ INC + HOVAL, data = NA.COL.OLD, listw = listw, -#> na.action = na.exclude) -#> Type: lag -#> -#> Coefficients: -#> rho (Intercept) INC HOVAL -#> 0.4537820 43.1054975 -0.9267352 -0.2715541 -#> -#> Log likelihood: -160.8867 resid(COL.lag.NA) -#> 1001 1002 1003 1004 1005 1006 -#> -4.4352945 -13.2750948 -2.9233555 -37.3067542 1.3568011 -3.8828027 -#> 1007 1008 1009 1010 1011 1012 -#> 5.9980436 -12.8113318 -4.0482558 18.1813323 5.9714203 1.0035599 -#> 1013 1014 1015 1016 1017 1018 -#> -1.7490666 -1.7490651 5.8456328 10.1074158 -4.3706895 3.3713771 -#> 1019 1020 1021 1022 1023 1024 -#> -4.0823022 NA NA NA NA NA -#> 1025 1026 1027 1028 1029 1030 -#> NA -6.0553296 -11.5813311 -8.0011389 -1.9047478 3.9744243 -#> 1031 1032 1033 1034 1035 1036 -#> 4.8148614 6.0673483 10.2059847 22.7419913 -2.0830998 0.1422777 -#> 1037 1038 1039 1040 1041 1042 -#> 8.6436388 8.8150864 6.4974685 15.4121789 9.4182148 5.6427603 -#> 1043 1044 1045 1046 1047 1048 -#> -7.5727123 -7.5983733 -9.7792620 -10.0870764 -3.1242393 3.4921796 -#> 1049 -#> 0.7173251 COL.lag.NA1 <- lagsarlm(CRIME ~ INC + HOVAL, data=NA.COL.OLD, listw, Durbin=~INC) # https://github.com/r-spatial/spatialreg/issues/10 COL.lag.NA1$na.action -#> 1020 1021 1022 1023 1024 1025 -#> 20 21 22 23 24 25 -#> attr(,"class") -#> [1] "omit" COL.lag.NA2 <- lagsarlm(CRIME ~ INC + HOVAL, data=NA.COL.OLD, listw, Durbin=~INC, na.action=na.exclude) COL.lag.NA2$na.action -#> 1020 1021 1022 1023 1024 1025 -#> 20 21 22 23 24 25 -#> attr(,"class") -#> [1] "exclude" # https://github.com/r-spatial/spatialreg/issues/11 COL.lag.NA3 <- lagsarlm(CRIME ~ INC + HOVAL, data=NA.COL.OLD, listw, control=list(pre_eig=ev)) -#> Warning: NAs found, precomputed eigenvalues ignored COL.lag.NA3$na.action -#> 1020 1021 1022 1023 1024 1025 -#> 20 21 22 23 24 25 -#> attr(,"class") -#> [1] "omit" -# } +} -# \dontrun{ +if (FALSE) { data(boston, package="spData") gp2mM <- lagsarlm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT), data=boston.c, spdep::nb2listw(boston.soi), type="mixed", method="Matrix") summary(gp2mM) -#> -#> Call:lagsarlm(formula = log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + -#> I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + -#> log(LSTAT), data = boston.c, listw = spdep::nb2listw(boston.soi), -#> type = "mixed", method = "Matrix") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -0.6316833 -0.0629790 -0.0090776 0.0682421 0.6991072 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 1.89816218 0.22759605 8.3400 < 2.2e-16 -#> CRIM -0.00571021 0.00093504 -6.1069 1.016e-09 -#> ZN 0.00069091 0.00051869 1.3320 0.182851 -#> INDUS -0.00111343 0.00307354 -0.3623 0.717155 -#> CHAS1 -0.04163225 0.02739364 -1.5198 0.128567 -#> I(NOX^2) -0.01034950 0.19360214 -0.0535 0.957367 -#> I(RM^2) 0.00794979 0.00102063 7.7891 6.661e-15 -#> AGE -0.00128789 0.00048920 -2.6326 0.008473 -#> log(DIS) -0.12404108 0.09510940 -1.3042 0.192168 -#> log(RAD) 0.05863502 0.02257078 2.5978 0.009382 -#> TAX -0.00049084 0.00012145 -4.0416 5.308e-05 -#> PTRATIO -0.01319853 0.00595352 -2.2169 0.026628 -#> B 0.00056383 0.00011089 5.0847 3.682e-07 -#> log(LSTAT) -0.24724454 0.02262033 -10.9302 < 2.2e-16 -#> lag.CRIM -0.00464215 0.00172935 -2.6843 0.007267 -#> lag.ZN -0.00037937 0.00070584 -0.5375 0.590940 -#> lag.INDUS 0.00025064 0.00385901 0.0649 0.948215 -#> lag.CHAS1 0.12518252 0.04071559 3.0746 0.002108 -#> lag.I(NOX^2) -0.38640401 0.22157523 -1.7439 0.081177 -#> lag.I(RM^2) -0.00451252 0.00153180 -2.9459 0.003220 -#> lag.AGE 0.00149678 0.00068418 2.1877 0.028693 -#> lag.log(DIS) -0.00453785 0.10046478 -0.0452 0.963973 -#> lag.log(RAD) -0.00940702 0.03104930 -0.3030 0.761912 -#> lag.TAX 0.00041083 0.00017867 2.2994 0.021481 -#> lag.PTRATIO 0.00060355 0.00788837 0.0765 0.939012 -#> lag.B -0.00050781 0.00014155 -3.5874 0.000334 -#> lag.log(LSTAT) 0.09846781 0.03399423 2.8966 0.003772 -#> -#> Rho: 0.59578, LR test value: 181.68, p-value: < 2.22e-16 -#> Asymptotic standard error: 0.038445 -#> z-value: 15.497, p-value: < 2.22e-16 -#> Wald statistic: 240.16, p-value: < 2.22e-16 -#> -#> Log likelihood: 300.6131 for mixed model -#> ML residual variance (sigma squared): 0.016011, (sigma: 0.12654) -#> Number of observations: 506 -#> Number of parameters estimated: 29 -#> AIC: -543.23, (AIC for lm: -363.55) -#> LM test for residual autocorrelation -#> test value: 29.772, p-value: 4.8604e-08 -#> W <- as(spdep::nb2listw(boston.soi), "CsparseMatrix") trMatb <- trW(W, type="mult") gp2mMi <- lagsarlm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + @@ -1074,61 +578,7 @@Examples
data=boston.c, spdep::nb2listw(boston.soi), type="mixed", method="Matrix", trs=trMatb) summary(gp2mMi) -#> -#> Call:lagsarlm(formula = log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) + -#> I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + -#> log(LSTAT), data = boston.c, listw = spdep::nb2listw(boston.soi), -#> type = "mixed", method = "Matrix", trs = trMatb) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -0.6316833 -0.0629790 -0.0090776 0.0682421 0.6991072 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 1.89816218 0.22759605 8.3400 < 2.2e-16 -#> CRIM -0.00571021 0.00093504 -6.1069 1.016e-09 -#> ZN 0.00069091 0.00051869 1.3320 0.182851 -#> INDUS -0.00111343 0.00307354 -0.3623 0.717155 -#> CHAS1 -0.04163225 0.02739364 -1.5198 0.128567 -#> I(NOX^2) -0.01034950 0.19360214 -0.0535 0.957367 -#> I(RM^2) 0.00794979 0.00102063 7.7891 6.661e-15 -#> AGE -0.00128789 0.00048920 -2.6326 0.008473 -#> log(DIS) -0.12404108 0.09510940 -1.3042 0.192168 -#> log(RAD) 0.05863502 0.02257078 2.5978 0.009382 -#> TAX -0.00049084 0.00012145 -4.0416 5.308e-05 -#> PTRATIO -0.01319853 0.00595352 -2.2169 0.026628 -#> B 0.00056383 0.00011089 5.0847 3.682e-07 -#> log(LSTAT) -0.24724454 0.02262033 -10.9302 < 2.2e-16 -#> lag.CRIM -0.00464215 0.00172935 -2.6843 0.007267 -#> lag.ZN -0.00037937 0.00070584 -0.5375 0.590940 -#> lag.INDUS 0.00025064 0.00385901 0.0649 0.948215 -#> lag.CHAS1 0.12518252 0.04071559 3.0746 0.002108 -#> lag.I(NOX^2) -0.38640401 0.22157523 -1.7439 0.081177 -#> lag.I(RM^2) -0.00451252 0.00153180 -2.9459 0.003220 -#> lag.AGE 0.00149678 0.00068418 2.1877 0.028693 -#> lag.log(DIS) -0.00453785 0.10046478 -0.0452 0.963973 -#> lag.log(RAD) -0.00940702 0.03104930 -0.3030 0.761912 -#> lag.TAX 0.00041083 0.00017867 2.2994 0.021481 -#> lag.PTRATIO 0.00060355 0.00788837 0.0765 0.939012 -#> lag.B -0.00050781 0.00014155 -3.5874 0.000334 -#> lag.log(LSTAT) 0.09846781 0.03399423 2.8966 0.003772 -#> -#> Rho: 0.59578, LR test value: 181.68, p-value: < 2.22e-16 -#> Asymptotic standard error: 0.038445 -#> z-value: 15.497, p-value: < 2.22e-16 -#> Wald statistic: 240.16, p-value: < 2.22e-16 -#> -#> Log likelihood: 300.6131 for mixed model -#> ML residual variance (sigma squared): 0.016011, (sigma: 0.12654) -#> Number of observations: 506 -#> Number of parameters estimated: 29 -#> AIC: -543.23, (AIC for lm: -363.55) -#> LM test for residual autocorrelation -#> test value: 29.772, p-value: 4.8604e-08 -#> -# } +} COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, quiet=FALSE, control=list(pre_eig=ev)) #> @@ -1147,8 +597,8 @@Examples
#> lambda: 0.5617888 function: -183.3805 Jacobian: -2.134769 SSE: 4683.153 #> lambda: 0.5617902 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 #> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 -#> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134783 SSE: 4683.151 -#> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 +#> lambda: 0.5617902 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 +#> lambda: 0.5617902 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 summary(COL.errW.eig) #> #> Call:errorsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, @@ -1161,7 +611,7 @@Examples
#> Type: error #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 59.893219 5.366163 11.1613 < 2.2e-16 +#> (Intercept) 59.893220 5.366162 11.1613 < 2.2e-16 #> INC -0.941312 0.330569 -2.8476 0.0044057 #> HOVAL -0.302250 0.090476 -3.3407 0.0008358 #> @@ -1228,10 +678,10 @@Examples
#> lambda: 0.5617866 function: -183.3805 Jacobian: -2.134749 SSE: 4683.157 #> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134783 SSE: 4683.15 #> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 -#> lambda: 0.5617902 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 #> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 #> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134783 SSE: 4683.151 -#> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134782 SSE: 4683.151 +#> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134783 SSE: 4683.151 +#> lambda: 0.5617903 function: -183.3805 Jacobian: -2.134783 SSE: 4683.151 summary(COL.errW.M) #> #> Call:errorsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, @@ -1290,215 +740,62 @@Examples
#> Number of parameters estimated: 7 #> AIC: 377.17, (AIC for lm: 380.16) #> -# \dontrun{ +if (FALSE) { COL.SDEM.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin=TRUE, control=list(pre_eig=ev)) summary(COL.SDEM.eig) -#> -#> Call:errorsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> Durbin = TRUE, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.31635 -6.54376 -0.22212 6.44591 23.15801 -#> -#> Type: error -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 73.545133 8.783543 8.3731 < 2.2e-16 -#> INC -1.051673 0.319514 -3.2915 0.0009966 -#> HOVAL -0.275608 0.091151 -3.0236 0.0024976 -#> lag.INC -1.156711 0.578629 -1.9991 0.0456024 -#> lag.HOVAL 0.111691 0.198993 0.5613 0.5746048 -#> -#> Lambda: 0.4254, LR test value: 4.9871, p-value: 0.025537 -#> Asymptotic standard error: 0.15842 -#> z-value: 2.6852, p-value: 0.0072485 -#> Wald statistic: 7.2103, p-value: 0.0072485 -#> -#> Log likelihood: -181.5846 for error model -#> ML residual variance (sigma squared): 92.531, (sigma: 9.6193) -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 377.17, (AIC for lm: 380.16) -#> COL.SDEM.eig <- errorsarlm(CRIME ~ DISCBD + INC + HOVAL, data=COL.OLD, listw, Durbin=~INC, control=list(pre_eig=ev)) summary(COL.SDEM.eig) -#> -#> Call:errorsarlm(formula = CRIME ~ DISCBD + INC + HOVAL, data = COL.OLD, -#> listw = listw, Durbin = ~INC, control = list(pre_eig = ev)) -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -34.61867 -7.31993 0.82879 5.92877 17.82211 -#> -#> Type: error -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 68.961912 6.985785 9.8717 < 2.2e-16 -#> DISCBD -5.412936 2.009281 -2.6940 0.007061 -#> INC -0.899425 0.315174 -2.8537 0.004321 -#> HOVAL -0.202846 0.090125 -2.2507 0.024403 -#> lag.INC 0.161858 0.603859 0.2680 0.788669 -#> -#> Lambda: 0.24524, LR test value: 1.2719, p-value: 0.25942 -#> Asymptotic standard error: 0.18393 -#> z-value: 1.3334, p-value: 0.18241 -#> Wald statistic: 1.7778, p-value: 0.18241 -#> -#> Log likelihood: -178.9895 for error model -#> ML residual variance (sigma squared): 85.935, (sigma: 9.2701) -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 371.98, (AIC for lm: 371.25) -#> summary(impacts(COL.SDEM.eig)) -#> Impact measures (SDEM, estimable, n): -#> Direct Indirect Total -#> DISCBD -5.4129364 NA -5.4129364 -#> INC -0.8994250 0.161858 -0.7375670 -#> HOVAL -0.2028457 NA -0.2028457 -#> ======================================================== -#> Standard errors: -#> Direct Indirect Total -#> DISCBD 2.00928149 NA 2.00928149 -#> INC 0.31517363 0.6038588 0.67559404 -#> HOVAL 0.09012493 NA 0.09012493 -#> ======================================================== -#> Z-values: -#> Direct Indirect Total -#> DISCBD -2.693966 NA -2.693966 -#> INC -2.853745 0.2680395 -1.091731 -#> HOVAL -2.250717 NA -2.250717 -#> -#> p-values: -#> Direct Indirect Total -#> DISCBD 0.0070607 NA 0.0070607 -#> INC 0.0043207 0.78867 0.2749513 -#> HOVAL 0.0244035 NA 0.0244035 -#> NA.COL.OLD <- COL.OLD NA.COL.OLD$CRIME[20:25] <- NA COL.err.NA <- errorsarlm(CRIME ~ INC + HOVAL, data=NA.COL.OLD, listw, na.action=na.exclude) COL.err.NA$na.action -#> 1020 1021 1022 1023 1024 1025 -#> 20 21 22 23 24 25 -#> attr(,"class") -#> [1] "exclude" COL.err.NA -#> -#> Call: -#> errorsarlm(formula = CRIME ~ INC + HOVAL, data = NA.COL.OLD, -#> listw = listw, na.action = na.exclude) -#> Type: error -#> -#> Coefficients: -#> lambda (Intercept) INC HOVAL -#> 0.5748430 58.2460528 -0.8473028 -0.3024909 -#> -#> Log likelihood: -161.8763 resid(COL.err.NA) -#> 1001 1002 1003 1004 1005 1006 -#> -4.18270830 -11.44133842 0.31874928 -34.47163074 2.42244758 -4.32095072 -#> 1007 1008 1009 1010 1011 1012 -#> 8.66744165 -13.38669935 -1.92276585 17.85753951 -1.11484596 -2.30434792 -#> 1013 1014 1015 1016 1017 1018 -#> -8.16935116 -5.80500231 0.14973721 5.93191444 -7.03028271 2.39112829 -#> 1019 1020 1021 1022 1023 1024 -#> -8.95099917 NA NA NA NA NA -#> 1025 1026 1027 1028 1029 1030 -#> NA -2.52940711 -9.60025384 -6.95635586 -0.43630579 5.98493664 -#> 1031 1032 1033 1034 1035 1036 -#> 6.25882669 7.75527031 10.83413235 23.23927250 -0.05594958 1.43808573 -#> 1037 1038 1039 1040 1041 1042 -#> 9.51995265 12.18295373 8.31031724 17.06834503 7.04418392 7.50088924 -#> 1043 1044 1045 1046 1047 1048 -#> -7.78485976 -6.79207446 -7.94977534 -11.25362117 -5.68994630 5.04837399 -#> 1049 -#> 2.22497381 print(system.time(ev <- eigenw(similar.listw(listw)))) -#> user system elapsed -#> 0.002 0.000 0.002 print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="eigen", control=list(pre_eig=ev)))) -#> user system elapsed -#> 0.198 0.000 0.199 ocoef <- coefficients(COL.errW.eig) print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="eigen", control=list(pre_eig=ev, LAPACK=FALSE)))) -#> user system elapsed -#> 0.195 0.000 0.196 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="eigen", control=list(pre_eig=ev, compiled_sse=TRUE)))) -#> user system elapsed -#> 0.192 0.000 0.193 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix_J", control=list(super=TRUE)))) -#> user system elapsed -#> 0.243 0.000 0.245 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix_J", control=list(super=FALSE)))) -#> user system elapsed -#> 0.238 0.000 0.239 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix_J", control=list(super=as.logical(NA))))) -#> user system elapsed -#> 0.238 0.000 0.239 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix", control=list(super=TRUE)))) -#> user system elapsed -#> 0.203 0.000 0.204 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix", control=list(super=FALSE)))) -#> user system elapsed -#> 0.202 0.000 0.203 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="Matrix", control=list(super=as.logical(NA))))) -#> user system elapsed -#> 0.200 0.000 0.202 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="spam", control=list(spamPivot="MMD")))) -#> user system elapsed -#> 0.201 0.000 0.201 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="spam", control=list(spamPivot="RCM")))) -#> user system elapsed -#> 0.202 0.000 0.204 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="spam_update", control=list(spamPivot="MMD")))) -#> user system elapsed -#> 0.201 0.000 0.201 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE print(system.time(COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, method="spam_update", control=list(spamPivot="RCM")))) -#> user system elapsed -#> 0.202 0.000 0.202 print(all.equal(ocoef, coefficients(COL.errW.eig))) -#> [1] TRUE -# } +} COL.sacW.eig <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, control=list(pre_eig1=ev, pre_eig2=ev)) summary(COL.sacW.eig) @@ -1513,7 +810,7 @@Examples
#> Type: sac #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 47.783766 9.902659 4.8253 1.398e-06 +#> (Intercept) 47.783764 9.902659 4.8253 1.398e-06 #> INC -1.025894 0.326326 -3.1438 0.001668 #> HOVAL -0.281651 0.090033 -3.1283 0.001758 #> @@ -1530,14 +827,14 @@Examples
#> ML residual variance (sigma squared): 95.604, (sigma: 9.7777) #> Number of observations: 49 #> Number of parameters estimated: 6 -#> AIC: 376.47, (AIC for lm: 382.75) +#> AIC: NA (not available for weighted model), (AIC for lm: 382.75) #> set.seed(1) summary(impacts(COL.sacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE) #> Impact measures (sac, trace): #> Direct Indirect Total -#> INC -1.0632723 -0.5601501 -1.6234223 -#> HOVAL -0.2919129 -0.1537847 -0.4456977 +#> INC -1.0632723 -0.5601502 -1.6234225 +#> HOVAL -0.2919129 -0.1537848 -0.4456977 #> ======================================================== #> Simulation results ( variance matrix): #> ======================================================== @@ -1548,13 +845,13 @@Examples
#> #> Simulated z-values: #> Direct Indirect Total -#> INC -3.376570 -0.8371024 -1.861677 -#> HOVAL -3.159909 -0.7431009 -1.553156 +#> INC -3.376570 -0.8371025 -1.861677 +#> HOVAL -3.159909 -0.7431010 -1.553156 #> #> Simulated p-values: #> Direct Indirect Total -#> INC 0.00073396 0.40254 0.062649 -#> HOVAL 0.00157818 0.45742 0.120386 +#> INC 0.00073396 0.40253 0.062649 +#> HOVAL 0.00157819 0.45742 0.120386 COL.msacW.eig <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, type="sacmixed", control=list(pre_eig1=ev, pre_eig2=ev)) summary(COL.msacW.eig) @@ -1569,7 +866,7 @@Examples
#> Type: sacmixed #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 50.92025 68.25721 0.7460 0.455664 +#> (Intercept) 50.92026 68.25721 0.7460 0.455664 #> INC -0.95072 0.44033 -2.1591 0.030841 #> HOVAL -0.28650 0.09994 -2.8667 0.004148 #> lag.INC -0.69261 1.69113 -0.4096 0.682132 @@ -1588,7 +885,7 @@Examples
#> ML residual variance (sigma squared): 93.149, (sigma: 9.6514) #> Number of observations: 49 #> Number of parameters estimated: 8 -#> AIC: 378.68, (AIC for lm: 382.75) +#> AIC: NA (not available for weighted model), (AIC for lm: 382.75) #> set.seed(1) summary(impacts(COL.msacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE) @@ -1601,13 +898,13 @@Examples
#> ======================================================== #> Simulated standard errors #> Direct Indirect Total -#> INC 0.3778064 2.1709292 2.3099451 -#> HOVAL 0.1144714 0.7736217 0.8389227 +#> INC 0.3778064 2.1709289 2.3099448 +#> HOVAL 0.1144714 0.7736216 0.8389226 #> #> Simulated z-values: #> Direct Indirect Total -#> INC -2.737182 -0.7297008 -1.13347001 -#> HOVAL -2.364128 0.2526688 -0.08958504 +#> INC -2.737182 -0.7297008 -1.13347010 +#> HOVAL -2.364129 0.2526689 -0.08958505 #> #> Simulated p-values: #> Direct Indirect Total @@ -1627,7 +924,7 @@Examples
#> Type: sacmixed #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 50.92025 68.25721 0.7460 0.455664 +#> (Intercept) 50.92026 68.25721 0.7460 0.455664 #> INC -0.95072 0.44033 -2.1591 0.030841 #> HOVAL -0.28650 0.09994 -2.8667 0.004148 #> lag.INC -0.69261 1.69113 -0.4096 0.682132 @@ -1646,7 +943,7 @@Examples
#> ML residual variance (sigma squared): 93.149, (sigma: 9.6514) #> Number of observations: 49 #> Number of parameters estimated: 8 -#> AIC: 378.68, (AIC for lm: 382.75) +#> AIC: NA (not available for weighted model), (AIC for lm: 382.75) #> set.seed(1) summary(impacts(COL.msacW1.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE) @@ -1659,13 +956,13 @@Examples
#> ======================================================== #> Simulated standard errors #> Direct Indirect Total -#> INC 0.3778064 2.1709292 2.3099451 -#> HOVAL 0.1144714 0.7736217 0.8389227 +#> INC 0.3778064 2.1709289 2.3099448 +#> HOVAL 0.1144714 0.7736216 0.8389226 #> #> Simulated z-values: #> Direct Indirect Total -#> INC -2.737182 -0.7297008 -1.13347001 -#> HOVAL -2.364128 0.2526688 -0.08958504 +#> INC -2.737182 -0.7297008 -1.13347010 +#> HOVAL -2.364129 0.2526689 -0.08958505 #> #> Simulated p-values: #> Direct Indirect Total @@ -1686,7 +983,7 @@Examples
#> Type: sacmixed #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 74.064506 37.738937 1.9625 0.049699 +#> (Intercept) 74.064502 37.738940 1.9625 0.049699 #> DISCBD -5.707678 3.248629 -1.7569 0.078926 #> INC -0.900975 0.314996 -2.8603 0.004233 #> HOVAL -0.203399 0.092982 -2.1875 0.028705 @@ -1705,161 +1002,45 @@Examples
#> ML residual variance (sigma squared): 85.135, (sigma: 9.2269) #> Number of observations: 49 #> Number of parameters estimated: 8 -#> AIC: 373.93, (AIC for lm: 369.42) +#> AIC: NA (not available for weighted model), (AIC for lm: 369.42) #> summary(impacts(COL.msacW2.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE) #> Impact measures (sacmixed, trace): #> Direct Indirect Total -#> DISCBD -5.7150801 0.41841708 -5.2966631 -#> INC -0.9031729 0.12421194 -0.7789609 -#> HOVAL -0.2036631 0.01491075 -0.1887524 +#> DISCBD -5.7150799 0.41841671 -5.2966632 +#> INC -0.9031729 0.12421198 -0.7789609 +#> HOVAL -0.2036631 0.01491074 -0.1887524 #> ======================================================== #> Simulation results ( variance matrix): #> ======================================================== #> Simulated standard errors -#> Direct Indirect Total -#> DISCBD 3.1420102 5.8423996 5.8501653 -#> INC 0.3651443 1.8447846 2.0077003 -#> HOVAL 0.1088636 0.6010839 0.6599894 +#> Direct Indirect Total +#> DISCBD 3.1446963 5.7201656 5.757459 +#> INC 0.3621137 1.7150155 1.876353 +#> HOVAL 0.1090172 0.6360488 0.698836 #> #> Simulated z-values: #> Direct Indirect Total -#> DISCBD -1.854977 0.0193511 -0.9769469 -#> INC -2.433297 0.1815217 -0.2757565 -#> HOVAL -2.034012 -0.1886907 -0.5073549 +#> DISCBD -1.852002 0.0280018 -0.9837341 +#> INC -2.461864 0.1861184 -0.3049951 +#> HOVAL -2.021743 -0.1841543 -0.4829973 #> #> Simulated p-values: #> Direct Indirect Total -#> DISCBD 0.063599 0.98456 0.32860 -#> INC 0.014962 0.85596 0.78274 -#> HOVAL 0.041950 0.85034 0.61191 -# \dontrun{ +#> DISCBD 0.064026 0.97766 0.32525 +#> INC 0.013822 0.85235 0.76037 +#> HOVAL 0.043203 0.85389 0.62910 +if (FALSE) { COL.mix.eig <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, type="mixed", method="eigen") summary(COL.mix.eig, correlation=TRUE, Nagelkerke=TRUE) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> type = "mixed", method = "eigen") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.47829 -6.46731 -0.33835 6.05200 22.62969 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 42.822415 12.667205 3.3806 0.0007233 -#> INC -0.914223 0.331094 -2.7612 0.0057586 -#> HOVAL -0.293738 0.089212 -3.2926 0.0009927 -#> lag.INC -0.520284 0.565129 -0.9206 0.3572355 -#> lag.HOVAL 0.245640 0.178917 1.3729 0.1697756 -#> -#> Rho: 0.42634, LR test value: 5.3693, p-value: 0.020494 -#> Asymptotic standard error: 0.15623 -#> z-value: 2.7288, p-value: 0.0063561 -#> Wald statistic: 7.4465, p-value: 0.0063561 -#> -#> Log likelihood: -181.3935 for mixed model -#> ML residual variance (sigma squared): 91.791, (sigma: 9.5808) -#> Nagelkerke pseudo-R-squared: 0.6494 -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 376.79, (AIC for lm: 380.16) -#> LM test for residual autocorrelation -#> test value: 0.28919, p-value: 0.59074 -#> -#> Correlation of coefficients -#> sigma rho (Intercept) INC HOVAL lag.INC -#> rho -0.18 -#> (Intercept) 0.16 -0.89 -#> INC -0.03 0.14 -0.19 -#> HOVAL 0.02 -0.09 0.03 -0.45 -#> lag.INC -0.09 0.49 -0.53 -0.36 0.05 -#> lag.HOVAL -0.04 0.19 -0.36 0.19 -0.24 -0.41 -#> COL.mix.M <- lagsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, type="mixed", method="Matrix") summary(COL.mix.M, correlation=TRUE, Nagelkerke=TRUE) -#> -#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, -#> type = "mixed", method = "Matrix") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -37.47829 -6.46731 -0.33835 6.05200 22.62969 -#> -#> Type: mixed -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 42.822413 12.667204 3.3806 0.0007233 -#> INC -0.914223 0.331094 -2.7612 0.0057586 -#> HOVAL -0.293738 0.089212 -3.2926 0.0009927 -#> lag.INC -0.520283 0.565129 -0.9206 0.3572355 -#> lag.HOVAL 0.245640 0.178917 1.3729 0.1697756 -#> -#> Rho: 0.42634, LR test value: 5.3693, p-value: 0.020494 -#> Asymptotic standard error: 0.15623 -#> z-value: 2.7288, p-value: 0.0063561 -#> Wald statistic: 7.4465, p-value: 0.0063561 -#> -#> Log likelihood: -181.3935 for mixed model -#> ML residual variance (sigma squared): 91.791, (sigma: 9.5808) -#> Nagelkerke pseudo-R-squared: 0.6494 -#> Number of observations: 49 -#> Number of parameters estimated: 7 -#> AIC: 376.79, (AIC for lm: 380.16) -#> LM test for residual autocorrelation -#> test value: 0.28919, p-value: 0.59074 -#> -#> Correlation of coefficients -#> sigma rho (Intercept) INC HOVAL lag.INC -#> rho -0.18 -#> (Intercept) 0.16 -0.89 -#> INC -0.03 0.14 -0.19 -#> HOVAL 0.02 -0.09 0.03 -0.45 -#> lag.INC -0.09 0.49 -0.53 -0.36 0.05 -#> lag.HOVAL -0.04 0.19 -0.36 0.19 -0.24 -0.41 -#> COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, spdep::nb2listw(COL.nb, style="W"), method="eigen") summary(COL.errW.eig, correlation=TRUE, Nagelkerke=TRUE, Hausman=TRUE) -#> -#> Call: -#> errorsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = spdep::nb2listw(COL.nb, -#> style = "W"), method = "eigen") -#> -#> Residuals: -#> Min 1Q Median 3Q Max -#> -34.81174 -6.44031 -0.72142 7.61476 23.33626 -#> -#> Type: error -#> Coefficients: (asymptotic standard errors) -#> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 59.893219 5.366163 11.1613 < 2.2e-16 -#> INC -0.941312 0.330569 -2.8476 0.0044057 -#> HOVAL -0.302250 0.090476 -3.3407 0.0008358 -#> -#> Lambda: 0.56179, LR test value: 7.9935, p-value: 0.0046945 -#> Asymptotic standard error: 0.13387 -#> z-value: 4.1966, p-value: 2.7098e-05 -#> Wald statistic: 17.611, p-value: 2.7098e-05 -#> -#> Log likelihood: -183.3805 for error model -#> ML residual variance (sigma squared): 95.575, (sigma: 9.7762) -#> Nagelkerke pseudo-R-squared: 0.61978 -#> Number of observations: 49 -#> Number of parameters estimated: 5 -#> AIC: 376.76, (AIC for lm: 382.75) -#> Hausman test: 4.902, df: 3, p-value: 0.17911 -#> -#> Correlation of coefficients -#> sigma lambda (Intercept) INC -#> lambda -0.24 -#> (Intercept) 0.00 0.00 -#> INC 0.00 0.00 -0.56 -#> HOVAL 0.00 0.00 -0.26 -0.45 -#> -# } +}Examples
Examples
#require("spdep", quietly=TRUE) data(oldcol, package="spdep") lw <- spdep::nb2listw(COL.nb, style="W") -ev <- eigenw(lw) -W <- as(lw, "CsparseMatrix") -trMatc <- trW(W, type="mult") require("coda", quietly=TRUE) set.seed(1) COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw) @@ -336,621 +333,65 @@
Examples
#> lambda 3 1010 937 1.080 #> sige 2 930 937 0.993 #> -# \dontrun{ +if (FALSE) { +ev <- eigenw(lw) +W <- as(lw, "CsparseMatrix") +trMatc <- trW(W, type="mult") set.seed(1) COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, control=list(prior=list(lambdaMH=TRUE))) print(summary(COL.err.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 59.5304 8.9916 0.201059 0.275885 -#> INC -0.9142 0.3969 0.008874 0.011266 -#> HOVAL -0.3027 0.1021 0.002282 0.002282 -#> lambda 0.6011 0.1502 0.003358 0.007723 -#> sige 116.6773 27.9707 0.625443 0.660902 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 44.0607 54.802 59.6796 64.2473 72.9189 -#> INC -1.6866 -1.184 -0.9177 -0.6422 -0.1429 -#> HOVAL -0.5102 -0.370 -0.3005 -0.2359 -0.1019 -#> lambda 0.2826 0.511 0.6135 0.7023 0.8676 -#> sige 74.1788 97.167 113.2955 130.2498 182.6024 -#> print(raftery.diag(COL.err.Bayes, r=0.01)) -#> -#> Quantile (q) = 0.025 -#> Accuracy (r) = +/- 0.01 -#> Probability (s) = 0.95 -#> -#> Burn-in Total Lower bound Dependence -#> (M) (N) (Nmin) factor (I) -#> (Intercept) 3 1143 937 1.220 -#> INC 2 969 937 1.030 -#> HOVAL 2 892 937 0.952 -#> lambda 17 4454 937 4.750 -#> sige 2 930 937 0.993 -#> set.seed(1) COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=TRUE) print(summary(COL.err.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 72.66455 11.7701 0.263187 0.263187 -#> INC -1.01543 0.3763 0.008414 0.008414 -#> HOVAL -0.28060 0.1050 0.002347 0.002347 -#> lag.INC -1.06754 0.7388 0.016519 0.016519 -#> lag.HOVAL 0.08961 0.2383 0.005328 0.005328 -#> lambda 0.48878 0.1732 0.003873 0.003873 -#> sige 114.66234 26.8527 0.600444 0.659903 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 48.2182 65.57218 73.00283 80.0588 95.67986 -#> INC -1.7264 -1.27249 -1.01848 -0.7658 -0.26325 -#> HOVAL -0.4835 -0.35298 -0.27998 -0.2105 -0.07685 -#> lag.INC -2.4553 -1.54934 -1.08873 -0.6092 0.39927 -#> lag.HOVAL -0.3926 -0.06243 0.09177 0.2443 0.55929 -#> lambda 0.1165 0.37669 0.49775 0.6118 0.79690 -#> sige 73.4827 95.70364 110.68368 129.1459 177.97468 -#> print(summary(impacts(COL.err.Bayes))) -#> Impact measures (SDEM, MCMC, n): -#> Direct Indirect Total -#> INC -1.0154294 -1.06753527 -2.0829647 -#> HOVAL -0.2806034 0.08960515 -0.1909982 -#> ======================================================== -#> Standard errors: -#> Direct Indirect Total -#> INC 0.35656492 0.7000652 0.8257127 -#> HOVAL 0.09946104 0.2258118 0.2716818 -#> ======================================================== -#> Z-values: -#> Direct Indirect Total -#> INC -2.847811 -1.5249083 -2.5226265 -#> HOVAL -2.821239 0.3968134 -0.7030219 -#> -#> p-values: -#> Direct Indirect Total -#> INC 0.0044021 0.12728 0.011648 -#> HOVAL 0.0047839 0.69151 0.482042 -#> print(raftery.diag(COL.err.Bayes, r=0.01)) -#> -#> Quantile (q) = 0.025 -#> Accuracy (r) = +/- 0.01 -#> Probability (s) = 0.95 -#> -#> Burn-in Total Lower bound Dependence -#> (M) (N) (Nmin) factor (I) -#> (Intercept) 2 930 937 0.993 -#> INC 3 1010 937 1.080 -#> HOVAL 2 930 937 0.993 -#> lag.INC 2 892 937 0.952 -#> lag.HOVAL 2 892 937 0.952 -#> lambda 2 930 937 0.993 -#> sige 2 930 937 0.993 -#> set.seed(1) COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=TRUE, control=list(prior=list(lambdaMH=TRUE))) print(summary(COL.err.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 70.40019 70.4345 1.574962 1.385430 -#> INC -0.97037 0.4110 0.009190 0.009823 -#> HOVAL -0.28455 0.1145 0.002561 0.002561 -#> lag.INC -0.95215 0.8426 0.018842 0.027001 -#> lag.HOVAL 0.06555 0.2679 0.005990 0.007555 -#> lambda 0.60702 0.1781 0.003983 0.014276 -#> sige 120.45188 29.1908 0.652727 1.279486 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 37.9440 62.5478 71.80756 80.9018 104.95272 -#> INC -1.7687 -1.2346 -0.98524 -0.7063 -0.14242 -#> HOVAL -0.5141 -0.3592 -0.28460 -0.2088 -0.06554 -#> lag.INC -2.6086 -1.5191 -0.95619 -0.4294 0.72268 -#> lag.HOVAL -0.4869 -0.1027 0.06799 0.2443 0.56913 -#> lambda 0.2366 0.5043 0.62080 0.7333 0.92570 -#> sige 76.7686 99.6498 116.80421 135.9061 189.12857 -#> print(summary(impacts(COL.err.Bayes))) -#> Impact measures (SDEM, MCMC, n): -#> Direct Indirect Total -#> INC -0.9703690 -0.95214894 -1.9225179 -#> HOVAL -0.2845537 0.06555386 -0.2189999 -#> ======================================================== -#> Standard errors: -#> Direct Indirect Total -#> INC 0.3894378 0.7984944 0.9933065 -#> HOVAL 0.1085247 0.2538476 0.3175629 -#> ======================================================== -#> Z-values: -#> Direct Indirect Total -#> INC -2.491717 -1.192430 -1.9354730 -#> HOVAL -2.622018 0.258241 -0.6896269 -#> -#> p-values: -#> Direct Indirect Total -#> INC 0.0127127 0.23309 0.052932 -#> HOVAL 0.0087411 0.79622 0.490429 -#> print(raftery.diag(COL.err.Bayes, r=0.01)) -#> -#> Quantile (q) = 0.025 -#> Accuracy (r) = +/- 0.01 -#> Probability (s) = 0.95 -#> -#> Burn-in Total Lower bound Dependence -#> (M) (N) (Nmin) factor (I) -#> (Intercept) 5 1472 937 1.570 -#> INC 2 892 937 0.952 -#> HOVAL 3 1010 937 1.080 -#> lag.INC 3 1010 937 1.080 -#> lag.HOVAL 2 930 937 0.993 -#> lambda 14 3645 937 3.890 -#> sige 2 930 937 0.993 -#> set.seed(1) COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=~INC) print(summary(COL.err.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 74.5423 10.98931 0.245728 0.245728 -#> INC -1.0339 0.38293 0.008563 0.008987 -#> HOVAL -0.2860 0.09857 0.002204 0.002204 -#> lag.INC -0.9190 0.59170 0.013231 0.012507 -#> lambda 0.4879 0.17196 0.003845 0.003672 -#> sige 112.0658 25.47022 0.569531 0.567228 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 52.0905 67.8363 74.9216 81.7016 94.97558 -#> INC -1.7796 -1.2965 -1.0425 -0.7816 -0.23925 -#> HOVAL -0.4762 -0.3494 -0.2861 -0.2205 -0.08885 -#> lag.INC -2.0736 -1.3191 -0.9217 -0.5457 0.27871 -#> lambda 0.1355 0.3737 0.4987 0.6088 0.80100 -#> sige 73.0813 94.1656 108.2980 125.8857 173.13082 -#> print(summary(impacts(COL.err.Bayes))) -#> Impact measures (SDEM, MCMC, n): -#> Direct Indirect Total -#> INC -1.0339443 -0.9190497 -1.9529940 -#> HOVAL -0.2860029 NA -0.2860029 -#> ======================================================== -#> Standard errors: -#> Direct Indirect Total -#> INC 0.36696535 0.5670367 0.69945670 -#> HOVAL 0.09446444 NA 0.09446444 -#> ======================================================== -#> Z-values: -#> Direct Indirect Total -#> INC -2.817553 -1.620794 -2.792159 -#> HOVAL -3.027625 NA -3.027625 -#> -#> p-values: -#> Direct Indirect Total -#> INC 0.0048391 0.10506 0.0052358 -#> HOVAL 0.0024648 NA 0.0024648 -#> print(raftery.diag(COL.err.Bayes, r=0.01)) -#> -#> Quantile (q) = 0.025 -#> Accuracy (r) = +/- 0.01 -#> Probability (s) = 0.95 -#> -#> Burn-in Total Lower bound Dependence -#> (M) (N) (Nmin) factor (I) -#> (Intercept) 3 1010 937 1.080 -#> INC 2 892 937 0.952 -#> HOVAL 2 930 937 0.993 -#> lag.INC 2 892 937 0.952 -#> lambda 2 892 937 0.952 -#> sige 3 1010 937 1.080 -#> set.seed(1) COL.err.Bayes <- spBreg_err(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=~INC, control=list(prior=list(lambdaMH=TRUE))) print(summary(COL.err.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 72.8864 13.2070 0.295318 0.332323 -#> INC -0.9843 0.3831 0.008566 0.009768 -#> HOVAL -0.2889 0.1008 0.002253 0.002253 -#> lag.INC -0.8606 0.6640 0.014847 0.014847 -#> lambda 0.5779 0.1662 0.003715 0.010526 -#> sige 114.9801 26.2240 0.586387 0.665850 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 45.5280 65.3717 73.3929 81.2895 97.2255 -#> INC -1.7170 -1.2373 -0.9898 -0.7261 -0.2197 -#> HOVAL -0.4870 -0.3565 -0.2850 -0.2177 -0.1032 -#> lag.INC -2.1579 -1.2914 -0.8847 -0.4386 0.4898 -#> lambda 0.2436 0.4725 0.5835 0.6981 0.8672 -#> sige 73.5737 96.1297 111.5039 130.5606 176.0332 -#> print(summary(impacts(COL.err.Bayes))) -#> Impact measures (SDEM, MCMC, n): -#> Direct Indirect Total -#> INC -0.9843411 -0.8605551 -1.8448962 -#> HOVAL -0.2889068 NA -0.2889068 -#> ======================================================== -#> Standard errors: -#> Direct Indirect Total -#> INC 0.36712699 0.6362974 0.80457985 -#> HOVAL 0.09657666 NA 0.09657666 -#> ======================================================== -#> Z-values: -#> Direct Indirect Total -#> INC -2.681201 -1.352442 -2.292993 -#> HOVAL -2.991476 NA -2.991476 -#> -#> p-values: -#> Direct Indirect Total -#> INC 0.0073359 0.17623 0.0218484 -#> HOVAL 0.0027763 NA 0.0027763 -#> print(raftery.diag(COL.err.Bayes, r=0.01)) -#> -#> Quantile (q) = 0.025 -#> Accuracy (r) = +/- 0.01 -#> Probability (s) = 0.95 -#> -#> Burn-in Total Lower bound Dependence -#> (M) (N) (Nmin) factor (I) -#> (Intercept) 4 1192 937 1.270 -#> INC 2 892 937 0.952 -#> HOVAL 2 930 937 0.993 -#> lag.INC 2 892 937 0.952 -#> lambda 10 2853 937 3.040 -#> sige 2 930 937 0.993 -#> set.seed(1) COL.sacW.B0 <- spBreg_sac(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=FALSE, control=list(ndraw=1500L, nomit=500L)) print(summary(COL.sacW.B0)) -#> -#> Iterations = 501:1500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 1000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 49.3672 9.23725 0.292108 0.847100 -#> INC -0.9925 0.35777 0.011314 0.012624 -#> HOVAL -0.2847 0.09798 0.003099 0.003307 -#> rho 0.3117 0.19633 0.006208 0.020392 -#> lambda 0.1927 0.27130 0.008579 0.027162 -#> sige 105.0462 23.96805 0.757936 0.757936 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 31.8783 43.12583 48.9398 55.3573 68.1858 -#> INC -1.6817 -1.24459 -0.9872 -0.7516 -0.3168 -#> HOVAL -0.4952 -0.34433 -0.2827 -0.2201 -0.1028 -#> rho -0.1548 0.21285 0.3406 0.4438 0.6176 -#> lambda -0.3788 0.01455 0.2171 0.3794 0.6988 -#> sige 68.5025 88.13909 101.6888 118.2545 160.2792 -#> print(summary(impacts(COL.sacW.B0, tr=trMatc), zstats=TRUE, short=TRUE)) -#> Impact measures (sac, trace): -#> Direct Indirect Total -#> INC -1.017325 -0.4246675 -1.4419930 -#> HOVAL -0.291816 -0.1218143 -0.4136303 -#> ======================================================== -#> Simulation results ( variance matrix): -#> ======================================================== -#> Simulated standard errors -#> Direct Indirect Total -#> INC 0.3668822 0.4180525 0.6533395 -#> HOVAL 0.1016518 0.1341689 0.1983041 -#> -#> Simulated z-values: -#> Direct Indirect Total -#> INC -2.808257 -1.219897 -2.357550 -#> HOVAL -2.910868 -1.106171 -2.240542 -#> -#> Simulated p-values: -#> Direct Indirect Total -#> INC 0.0049810 0.22250 0.018396 -#> HOVAL 0.0036043 0.26865 0.025056 set.seed(1) COL.sacW.B1 <- spBreg_sac(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=TRUE, control=list(ndraw=1500L, nomit=500L)) print(summary(COL.sacW.B1)) -#> -#> Iterations = 501:1500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 1000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 60.1833 24.14827 0.763635 3.169760 -#> INC -0.9880 0.36251 0.011464 0.015819 -#> HOVAL -0.2839 0.09749 0.003083 0.003433 -#> lag.INC -0.8740 0.77272 0.024435 0.061565 -#> lag.HOVAL 0.1704 0.22184 0.007015 0.013821 -#> rho 0.1808 0.32007 0.010121 0.044327 -#> lambda 0.2106 0.32071 0.010142 0.041657 -#> sige 100.9463 24.32264 0.769149 0.946662 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 22.2042 42.277313 56.1007 76.4887 113.5316 -#> INC -1.7045 -1.223372 -0.9842 -0.7329 -0.3299 -#> HOVAL -0.4803 -0.349479 -0.2841 -0.2152 -0.1085 -#> lag.INC -2.5306 -1.336778 -0.8386 -0.3451 0.5156 -#> lag.HOVAL -0.3031 0.026900 0.1851 0.3249 0.5752 -#> rho -0.5076 -0.021227 0.2386 0.4320 0.6655 -#> lambda -0.5106 -0.007594 0.2129 0.4482 0.7389 -#> sige 64.6867 84.395156 96.9599 114.4101 159.7252 -#> print(summary(impacts(COL.sacW.B1, tr=trMatc), zstats=TRUE, short=TRUE)) -#> Impact measures (sacmixed, trace): -#> Direct Indirect Total -#> INC -1.0332779 -1.239705 -2.2729830 -#> HOVAL -0.2787547 0.140259 -0.1384957 -#> ======================================================== -#> Simulation results ( variance matrix): -#> ======================================================== -#> Simulated standard errors -#> Direct Indirect Total -#> INC 0.3510687 0.9383466 1.0337787 -#> HOVAL 0.1002385 0.2776460 0.3184971 -#> -#> Simulated z-values: -#> Direct Indirect Total -#> INC -2.962199 -1.4351185 -2.3085930 -#> HOVAL -2.775887 0.5415921 -0.4015108 -#> -#> Simulated p-values: -#> Direct Indirect Total -#> INC 0.0030545 0.15125 0.020966 -#> HOVAL 0.0055051 0.58810 0.688044 set.seed(1) COL.lag.Bayes <- spBreg_lag(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw) print(summary(COL.lag.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 45.8256 8.04911 0.179984 0.179984 -#> INC -1.0459 0.34746 0.007770 0.007770 -#> HOVAL -0.2662 0.09366 0.002094 0.002094 -#> rho 0.4146 0.12434 0.002780 0.002780 -#> sige 107.8827 23.75095 0.531087 0.561152 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 30.3334 40.3140 45.9376 51.0903 61.51008 -#> INC -1.7252 -1.2774 -1.0477 -0.8158 -0.35423 -#> HOVAL -0.4533 -0.3312 -0.2686 -0.2019 -0.08107 -#> rho 0.1645 0.3327 0.4157 0.5008 0.64982 -#> sige 70.5778 90.6316 104.8556 120.4706 164.75421 -#> print(summary(impacts(COL.lag.Bayes, tr=trMatc), short=TRUE, zstats=TRUE)) -#> Impact measures (lag, trace): -#> Direct Indirect Total -#> INC -1.0962113 -0.6903627 -1.7865741 -#> HOVAL -0.2790451 -0.1757347 -0.4547797 -#> ======================================================== -#> Simulation results ( variance matrix): -#> ======================================================== -#> Simulated standard errors -#> Direct Indirect Total -#> INC 0.3680189 0.4927553 0.7674941 -#> HOVAL 0.1002735 0.1393959 0.2171077 -#> -#> Simulated z-values: -#> Direct Indirect Total -#> INC -3.003178 -1.554252 -2.437925 -#> HOVAL -2.809370 -1.425055 -2.212507 -#> -#> Simulated p-values: -#> Direct Indirect Total -#> INC 0.0026718 0.12012 0.014772 -#> HOVAL 0.0049639 0.15414 0.026932 print(summary(impacts(COL.lag.Bayes, evalues=ev), short=TRUE, zstats=TRUE)) -#> Impact measures (lag, evalues): -#> Direct Indirect Total -#> INC -1.0962113 -0.6903627 -1.7865741 -#> HOVAL -0.2790451 -0.1757347 -0.4547797 -#> ======================================================== -#> Simulation results ( variance matrix): -#> ======================================================== -#> Simulated standard errors -#> Direct Indirect Total -#> INC 0.3680193 0.4929524 0.7676307 -#> HOVAL 0.1002746 0.1395496 0.2172252 -#> -#> Simulated z-values: -#> Direct Indirect Total -#> INC -3.003178 -1.553686 -2.437527 -#> HOVAL -2.809345 -1.423571 -2.211367 -#> -#> Simulated p-values: -#> Direct Indirect Total -#> INC 0.0026718 0.12026 0.014788 -#> HOVAL 0.0049642 0.15457 0.027010 set.seed(1) COL.D0.Bayes <- spBreg_lag(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw, Durbin=TRUE) print(summary(COL.D0.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 45.3106 13.74660 0.307383 0.307383 -#> INC -0.9255 0.36779 0.008224 0.008224 -#> HOVAL -0.2941 0.09628 0.002153 0.002153 -#> lag.INC -0.5834 0.63860 0.014279 0.014279 -#> lag.HOVAL 0.2396 0.19206 0.004295 0.004295 -#> rho 0.3932 0.16296 0.003644 0.003644 -#> sige 109.7428 24.97328 0.558419 0.607251 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 19.18523 35.9952 45.0347 54.5253 72.8545 -#> INC -1.63995 -1.1761 -0.9206 -0.6651 -0.2148 -#> HOVAL -0.48230 -0.3566 -0.2931 -0.2308 -0.1021 -#> lag.INC -1.82283 -1.0066 -0.5712 -0.1763 0.6692 -#> lag.HOVAL -0.14391 0.1068 0.2393 0.3670 0.6168 -#> rho 0.06053 0.2826 0.3977 0.5078 0.6949 -#> sige 71.19421 91.4872 106.1690 123.7826 170.9503 -#> print(summary(impacts(COL.D0.Bayes, tr=trMatc), short=TRUE, zstats=TRUE)) -#> Impact measures (mixed, trace): -#> Direct Indirect Total -#> INC -1.027773 -1.4589793 -2.48675192 -#> HOVAL -0.280732 0.1908138 -0.08991828 -#> ======================================================== -#> Simulation results ( variance matrix): -#> ======================================================== -#> Simulated standard errors -#> Direct Indirect Total -#> INC 0.38264314 1.3114491 1.4534522 -#> HOVAL 0.09937135 0.3347428 0.3698437 -#> -#> Simulated z-values: -#> Direct Indirect Total -#> INC -2.739321 -1.2582210 -1.8564594 -#> HOVAL -2.838369 0.5610536 -0.2548209 -#> -#> Simulated p-values: -#> Direct Indirect Total -#> INC 0.0061566 0.20831 0.063388 -#> HOVAL 0.0045345 0.57476 0.798861 set.seed(1) COL.D1.Bayes <- spBreg_lag(CRIME ~ DISCBD + INC + HOVAL, data=COL.OLD, listw=lw, Durbin= ~ INC) print(summary(COL.D1.Bayes)) -#> -#> Iterations = 501:2500 -#> Thinning interval = 1 -#> Number of chains = 1 -#> Sample size per chain = 2000 -#> -#> 1. Empirical mean and standard deviation for each variable, -#> plus standard error of the mean: -#> -#> Mean SD Naive SE Time-series SE -#> (Intercept) 56.4447 13.43635 0.300446 0.300446 -#> DISCBD -4.7748 2.14395 0.047940 0.047940 -#> INC -0.9493 0.34681 0.007755 0.007755 -#> HOVAL -0.1798 0.09983 0.002232 0.002232 -#> lag.INC 0.4262 0.63352 0.014166 0.014166 -#> rho 0.1873 0.18174 0.004064 0.004064 -#> sige 103.8701 23.29384 0.520866 0.566802 -#> -#> 2. Quantiles for each variable: -#> -#> 2.5% 25% 50% 75% 97.5% -#> (Intercept) 30.5093 47.769573 56.4789 65.0877 84.86254 -#> DISCBD -9.2416 -6.153288 -4.7453 -3.3591 -0.67219 -#> INC -1.6655 -1.174948 -0.9451 -0.7102 -0.28365 -#> HOVAL -0.3744 -0.247270 -0.1775 -0.1128 0.01524 -#> lag.INC -0.8028 0.004752 0.4224 0.8433 1.66618 -#> rho -0.2116 0.072536 0.1836 0.3067 0.53579 -#> sige 67.8946 87.054945 100.9815 116.4528 161.25826 -#> print(summary(impacts(COL.D1.Bayes, tr=trMatc), short=TRUE, zstats=TRUE)) -#> Impact measures (mixed, trace): -#> Direct Indirect Total -#> DISCBD -4.8148501 -1.06067837 -5.8755285 -#> INC -0.9381442 0.29454370 -0.6436005 -#> HOVAL -0.1813164 -0.03994275 -0.2212591 -#> ======================================================== -#> Simulation results ( variance matrix): -#> ======================================================== -#> Simulated standard errors -#> Direct Indirect Total -#> DISCBD 2.1926383 1.82145831 3.3798560 -#> INC 0.3502903 0.81693522 0.8950370 -#> HOVAL 0.1015612 0.06692703 0.1429623 -#> -#> Simulated z-values: -#> Direct Indirect Total -#> DISCBD -2.219068 -0.7449248 -1.8410439 -#> INC -2.696796 0.3333552 -0.7511779 -#> HOVAL -1.802805 -0.7412143 -1.6277185 -#> -#> Simulated p-values: -#> Direct Indirect Total -#> DISCBD 0.026482 0.45632 0.065615 -#> INC 0.007001 0.73887 0.452546 -#> HOVAL 0.071419 0.45856 0.103585 #data(elect80, package="spData") #lw <- spdep::nb2listw(e80_queen, zero.policy=TRUE) #el_ml <- lagsarlm(log(pc_turnout) ~ log(pc_college) + log(pc_homeownership) @@ -962,7 +403,7 @@Examples
#print(summary(el_B)) #print(el_ml$timings) #print(attr(el_B, "timings")) -# } +}Examples
lmSLX(formula, data = list(), listw, na.action, weights=NULL,
- Durbin=TRUE, zero.policy=NULL)
+ lmSLX(formula, data = list(), listw, na.action, weights=NULL, Durbin=TRUE,
+ zero.policy=NULL)
create_WX(x, listw, zero.policy=NULL, prefix="")
# S3 method for SlX
impacts(obj, ...)
@@ -106,8 +106,7 @@ Arguments
default TRUE for lmSLX
(Durbin model including WX); if TRUE, full spatial Durbin model; if a formula object, the subset of explanatory variables to lag
zero.policy
-default NULL, use global option value; if TRUE assign zero to the lagged value of zones without
-neighbours, if FALSE assign NA
+default NULL, use global option value; if TRUE assign zero to the lagged value of zones without neighbours, if FALSE assign NA
obj
A spatial regression object created by lmSLX
@@ -174,7 +173,7 @@ Examples
#> F-statistic: 17.12 on 4 and 44 DF, p-value: 1.553e-08
#>
summary(impacts(COL.SLX))
-#> Impact measures (SlX, estimable, n-k):
+#> Impact measures (SlX, glht, n-k):
#> Direct Indirect Total
#> INC -1.1089293 -1.3709725 -2.47990173
#> HOVAL -0.2897283 0.1917608 -0.09796753
@@ -196,7 +195,7 @@ Examples
#>
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL + I(HOVAL^2), data=COL.OLD, listw=lw, Durbin=TRUE)
summary(impacts(COL.SLX))
-#> Impact measures (SlX, estimable, n-k):
+#> Impact measures (SlX, glht, n-k):
#> Direct Indirect Total
#> INC -0.947594274 -1.275338647 -2.22293292
#> HOVAL -0.777427839 -0.355048446 -1.13247628
@@ -248,7 +247,7 @@ Examples
#>
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL + I(HOVAL^2), data=COL.OLD, listw=lw, Durbin=~INC)
summary(impacts(COL.SLX))
-#> Impact measures (SlX, estimable, n-k):
+#> Impact measures (SlX, glht, n-k):
#> Direct Indirect Total
#> INC -1.079064628 -1.010896 -2.089960575
#> HOVAL -0.634518755 NA -0.634518755
@@ -320,7 +319,7 @@ Examples
#> F-statistic: 26.49 on 2 and 46 DF, p-value: 2.214e-08
#>
summary(impacts(COL.SLX))
-#> Impact measures (SlX, estimable, n-k):
+#> Impact measures (SlX, glht, n-k):
#> Direct Indirect Total
#> INC -1.588901 -1.085867 -2.674768
#> ========================================================
@@ -336,7 +335,7 @@ Examples
#> Direct Indirect Total
#> INC 8.2671e-06 0.024029 5.1387e-11
#>
-# \dontrun{
+if (FALSE) {
crds <- cbind(COL.OLD$X, COL.OLD$Y)
mdist <- sqrt(sum(diff(apply(crds, 2, range))^2))
dnb <- spdep::dnearneigh(crds, 0, mdist)
@@ -350,66 +349,12 @@ Examples
}
opt <- optimize(f, interval=c(0.1, 4), form=CRIME ~ INC + HOVAL,
data=COL.OLD, dnb=dnb, dists=dists, verbose=TRUE, maximum=TRUE)
-#> power: 1.589667 logLik: -172.6864
-#> power: 2.510333 logLik: -177.741
-#> power: 1.020665 logLik: -171.5379
-#> power: 0.8721475 logLik: -171.7979
-#> power: 1.11302 logLik: -171.4973
-#> power: 1.107329 logLik: -171.497
-#> power: 1.105705 logLik: -171.497
-#> power: 1.105746 logLik: -171.497
-#> power: 1.105664 logLik: -171.497
-#> power: 1.105705 logLik: -171.497
glst <- lapply(dists, function(d) 1/(d^opt$maximum))
lw <- spdep::nb2listw(dnb, glist=glst, style="B")
SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw)
summary(SLX)
-#>
-#> Call:
-#> lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))),
-#> data = as.data.frame(x), weights = weights)
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -25.1535 -4.7529 -0.9819 4.8053 22.2274
-#>
-#> Coefficients:
-#> Estimate Std. Error t value Pr(>|t|)
-#> (Intercept) 18.78498 11.23795 1.672 0.1019
-#> INC -0.65428 0.29463 -2.221 0.0317 *
-#> HOVAL -0.18244 0.08138 -2.242 0.0302 *
-#> lag..Intercept. 6.27790 4.44284 1.413 0.1648
-#> lag.INC -0.12958 0.28794 -0.450 0.6550
-#> lag.HOVAL 0.02668 0.11273 0.237 0.8141
-#> ---
-#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-#>
-#> Residual standard error: 8.553 on 43 degrees of freedom
-#> Multiple R-squared: 0.7659, Adjusted R-squared: 0.7387
-#> F-statistic: 28.14 on 5 and 43 DF, p-value: 1.544e-12
-#>
summary(impacts(SLX))
-#> Impact measures (SlX, estimable, n-k):
-#> Direct Indirect Total
-#> INC -0.6542760 -0.12957739 -0.7838534
-#> HOVAL -0.1824383 0.02667561 -0.1557627
-#> ========================================================
-#> Standard errors:
-#> Direct Indirect Total
-#> INC 0.29463342 0.2879383 0.3742251
-#> HOVAL 0.08138257 0.1127279 0.1404591
-#> ========================================================
-#> Z-values:
-#> Direct Indirect Total
-#> INC -2.220644 -0.4500179 -2.094604
-#> HOVAL -2.241736 0.2366371 -1.108954
-#>
-#> p-values:
-#> Direct Indirect Total
-#> INC 0.026375 0.65270 0.036206
-#> HOVAL 0.024978 0.81294 0.267450
-#>
-# }
+}
COL.SLX <- lmSLX(CRIME ~ INC + HOVAL, data=COL.OLD, listw=lw)
pslx0 <- predict(COL.SLX)
pslx1 <- predict(COL.SLX, newdata=COL.OLD, listw=lw)
@@ -419,9 +364,9 @@ Examples
COL.OLD1$INC <- COL.OLD1$INC + 1
pslx2 <- predict(COL.SLX, newdata=COL.OLD1, listw=lw)
sum(coef(COL.SLX)[c(2,4)])
-#> [1] 5.623621
+#> [1] -2.479902
mean(pslx2-pslx1)
-#> [1] -1.425045
+#> [1] -2.479902
wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
library(spdep)
-#> Loading required package: sp
#>
#> Attaching package: ‘spdep’
#> The following objects are masked from ‘package:spatialreg’:
@@ -148,20 +147,9 @@ Examples
#>
aple(as.vector(scale(wheat$yield_detrend, scale=FALSE)), spdep::nb2listw(nbr12, style="W"))
#> [1] 0.6601805
-# \dontrun{
+if (FALSE) {
errorsarlm(yield_detrend ~ 1, wheat, spdep::nb2listw(nbr12, style="W"))
-#>
-#> Call:
-#> errorsarlm(formula = yield_detrend ~ 1, data = wheat, listw = spdep::nb2listw(nbr12,
-#> style = "W"))
-#> Type: error
-#>
-#> Coefficients:
-#> lambda (Intercept)
-#> 0.60189686 -0.00251772
-#>
-#> Log likelihood: -192.9519
-# }
+}
# \dontrun{
+ if (FALSE) {
wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
nbr1 <- spdep::poly2nb(wheat, queen=FALSE)
nbrl <- spdep::nblag(nbr1, 2)
@@ -127,23 +127,10 @@ Examples
boot_out_ser <- aple.mc(as.vector(scale(wheat$yield_detrend, scale=FALSE)),
spdep::nb2listw(nbr12, style="W"), nsim=500)
plot(boot_out_ser)
-
boot_out_ser
-#>
-#> DATA PERMUTATION
-#>
-#>
-#> Call:
-#> boot(data = x, statistic = aple.boot, R = nsim, sim = "permutation",
-#> pre = pre, parallel = parallel, ncpus = ncpus, cl = cl)
-#>
-#>
-#> Bootstrap Statistics :
-#> original bias std. error
-#> t1* 0.6601805 -0.6708415 0.1092922
library(parallel)
oldCores <- set.coresOption(NULL)
-nc <- detectCores(logical=FALSE)
+nc <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L
# set nc to 1L here
if (nc > 1L) nc <- 1L
invisible(set.coresOption(nc))
@@ -161,21 +148,9 @@ Examples
stopCluster(cl)
}
boot_out_par
-#>
-#> DATA PERMUTATION
-#>
-#>
-#> Call:
-#> boot(data = x, statistic = aple.boot, R = nsim, sim = "permutation",
-#> pre = pre, parallel = parallel, ncpus = ncpus, cl = cl)
-#>
-#>
-#> Bootstrap Statistics :
-#> original bias std. error
-#> t1* 0.6601805 -0.6708415 0.1092922
invisible(set.coresOption(oldCores))
RNGkind(oldRNG[1], oldRNG[2])
-# }
+}
# \dontrun{
+ if (FALSE) {
wheat <- st_read(system.file("shapes/wheat.shp", package="spData")[1], quiet=TRUE)
nbr1 <- spdep::poly2nb(wheat, queen=FALSE)
nbrl <- spdep::nblag(nbr1, 2)
@@ -141,72 +141,18 @@ Examples
abline(lm_obj)
abline(v=0, h=0, lty=2)
zz <- summary(influence.measures(lm_obj))
-#> Potentially influential observations of
-#> lm(formula = Y ~ X, data = plt_out) :
-#>
-#> dfb.1_ dfb.X dffit cov.r cook.d hat
-#> 34 -0.12 0.01 -0.12 0.98_* 0.01 0.00
-#> 50 0.10 0.01 0.11 0.98_* 0.01 0.00
-#> 60 0.00 0.00 0.00 1.01_* 0.00 0.01
-#> 118 0.01 0.01 0.02 1.01_* 0.00 0.01
-#> 137 -0.10 -0.01 -0.10 0.98_* 0.01 0.00
-#> 143 -0.02 -0.04 -0.05 1.01_* 0.00 0.01
-#> 157 -0.10 -0.05 -0.11 0.99_* 0.01 0.00
-#> 166 0.00 0.00 0.00 1.02_* 0.00 0.02_*
-#> 168 0.01 0.02 0.03 1.01_* 0.00 0.01
-#> 176 -0.10 0.17 -0.20_* 0.99 0.02 0.01
-#> 177 0.03 -0.07 0.08 1.01_* 0.00 0.01
-#> 191 0.01 0.04 0.04 1.02_* 0.00 0.02_*
-#> 192 0.10 0.18 0.20_* 0.99 0.02 0.01
-#> 201 -0.10 0.07 -0.12 0.99_* 0.01 0.00
-#> 216 0.02 0.05 0.05 1.02_* 0.00 0.01_*
-#> 217 0.03 0.08 0.09 1.02_* 0.00 0.01_*
-#> 225 -0.11 -0.07 -0.13 0.98_* 0.01 0.00
-#> 237 -0.10 0.02 -0.10 0.99_* 0.01 0.00
-#> 242 -0.02 -0.04 -0.04 1.01_* 0.00 0.01
-#> 287 0.14 -0.23 0.27_* 0.97_* 0.04 0.01
-#> 290 -0.18 -0.35 -0.40_* 0.95_* 0.08 0.01
-#> 295 -0.01 -0.01 -0.01 1.01_* 0.00 0.01
-#> 322 0.00 0.00 0.00 1.01_* 0.00 0.01
-#> 325 -0.10 -0.07 -0.12 0.99_* 0.01 0.00
-#> 351 0.19 -0.04 0.20_* 0.94_* 0.02 0.00
-#> 369 0.01 -0.03 0.03 1.01_* 0.00 0.01
-#> 376 -0.05 -0.13 -0.14 1.02_* 0.01 0.02_*
-#> 392 -0.04 0.08 -0.09 1.01_* 0.00 0.01
-#> 393 -0.03 0.06 -0.07 1.01_* 0.00 0.01
-#> 394 -0.01 0.03 -0.03 1.01_* 0.00 0.01
-#> 402 0.10 0.02 0.10 0.99_* 0.01 0.00
-#> 429 0.13 -0.10 0.16 0.98_* 0.01 0.00
-#> 430 -0.11 -0.23 -0.25_* 0.99 0.03 0.01
-#> 438 0.00 -0.01 -0.01 1.01_* 0.00 0.01
-#> 442 -0.01 0.04 -0.04 1.02_* 0.00 0.02_*
-#> 443 -0.01 0.02 -0.02 1.02_* 0.00 0.01_*
-#> 461 0.02 0.04 0.04 1.01_* 0.00 0.01
-#> 462 0.01 0.03 0.03 1.03_* 0.00 0.02_*
-#> 466 0.01 -0.02 0.02 1.02_* 0.00 0.01_*
-#> 467 -0.03 0.08 -0.08 1.02_* 0.00 0.01_*
-#> 468 0.02 -0.04 0.04 1.01_* 0.00 0.01
-#> 480 0.13 0.05 0.14 0.97_* 0.01 0.00
-#> 488 0.16 0.09 0.18 0.96_* 0.02 0.00
-#> 492 -0.13 0.12 -0.17 0.98_* 0.01 0.00
infl <- as.integer(rownames(zz))
points(plt_out$X[infl], plt_out$Y[infl], pch=3, cex=0.6, col="red")
-
crossprod(plt_out$Y, plt_out$X)/crossprod(plt_out$X)
-#> [,1]
-#> [1,] 0.6601805
wheat$localAple <- localAple(as.vector(scale(wheat$yield_detrend, scale=FALSE)),
spdep::nb2listw(nbr12, style="W"))
mean(wheat$localAple)
-#> [1] 0.6601805
hist(wheat$localAple)
-
opar <- par(no.readonly=TRUE)
plot(wheat[,"localAple"], reset=FALSE)
text(st_coordinates(st_centroid(st_geometry(wheat)))[infl,], labels=rep("*", length(infl)))
-
par(opar)
-# }
+}
default TRUE: truncate Smirnov correction term, see trW
default TRUE
use equation 7 in Smirnov and Anselin (2009), if FALSE no unit root correction
-default TRUE; use equation 7 in Smirnov and Anselin (2009), if FALSE no unit root correction
default “LU”, alternatively “MC”; underlying lndet method to use for generating SE toolbox emulation grid
gstsls(formula, data = list(), listw, listw2 = NULL, na.action = na.fail,
- zero.policy = NULL, pars=NULL, scaleU=FALSE, control = list(),
+ zero.policy = attr(listw, "zero.policy"), pars=NULL, scaleU=FALSE, control = list(),
verbose=NULL, method="nlminb", robust=FALSE, legacy=FALSE, W2X=TRUE)
# S3 method for Gmsar
impacts(obj, ..., n = NULL, tr = NULL, R = NULL,
@@ -327,7 +327,7 @@ Examples
a listw
object created for example by nb2listw
a listw
object created for example by spdep::nb2listw()
default NULL, use global option value; if TRUE assign zero to the lagged value of zones without @@ -118,9 +118,9 @@
The underlying spatial lag model:
$$y = \rho W y + X \beta + \varepsilon$$
-where \(\rho\) is the spatial parameter may be fitted by maximum likelihood. In that case, the log likelihood function includes the logartithm of cumbersome Jacobian term \(|I - \rho W|\). If we rewrite the model as:
+where \(\rho\) is the spatial parameter may be fitted by maximum likelihood. In that case, the log likelihood function includes the logarithm of cumbersome Jacobian term \(|I - \rho W|\). If we rewrite the model as:
$$S y = X \beta + \varepsilon$$
-we see that in the ML case \(S y = (I - \rho W) y\). If W is row-stochastic, S may be expressed as a linear combination of row-stochastic matrices. By pre-computing the matrix \([y Wy, W^2y, ..., W^{q-1}y]\), the term \(S y (\alpha)\) can readily be found by numerical optimization using the matrix exponential approach. \(\alpha\) and \(\rho\) are related as \(\rho = 1 - \exp{\alpha}\), conditional on the number of matrix power terms taken q
.
we see that in the ML case \(S y = (I - \rho W) y\). If W is row-stochastic, S may be expressed as a linear combination of row-stochastic matrices. By pre-computing the matrix \([y, Wy, W^2y, ..., W^{q-1}y]\), the term \(S y (\alpha)\) can readily be found by numerical optimization using the matrix exponential approach. \(\alpha\) and \(\rho\) are related as \(\rho = 1 - \exp{\alpha}\), conditional on the number of matrix power terms taken q
.
# \dontrun{
+ if (FALSE) {
require(sf, quietly=TRUE)
columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
#require(spdep, quietly=TRUE)
@@ -125,24 +125,6 @@ Examples
col.sp <- as.spam.listw(col.listw)
str(col.sp)
}
-#> Spam version 2.9-1 (2022-08-07) is loaded.
-#> Type 'help( Spam)' or 'demo( spam)' for a short introduction
-#> and overview of this package.
-#> Help for individual functions is also obtained by adding the
-#> suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
-#>
-#> Attaching package: ‘spam’
-#> The following object is masked from ‘package:Matrix’:
-#>
-#> det
-#> The following objects are masked from ‘package:base’:
-#>
-#> backsolve, forwardsolve
-#> Formal class 'spam' [package "spam"] with 4 slots
-#> ..@ entries : num [1:230] 0.5 0.5 0.333 0.333 0.333 ...
-#> ..@ colindices : int [1:230] 2 3 1 3 4 1 2 4 5 2 ...
-#> ..@ rowpointers: int [1:50] 1 3 6 10 14 21 23 27 33 41 ...
-#> ..@ dimension : int [1:2] 49 49
suppressMessages(nyadjmat <- as.matrix(foreign::read.dbf(system.file(
"misc/nyadjwts.dbf", package="spData")[1])[-1]))
nyadjlw <- spdep::mat2listw(nyadjmat)
@@ -154,14 +136,11 @@ Examples
I <- Diagonal(n)
rho <- 0.1
c(determinant(I - rho * W_S, logarithm=TRUE)$modulus)
-#> [1] -9.587255
sum(log(1 - rho * eigenw(listw_NY)))
-#> [1] -9.587255
nW <- - W_S
nChol <- Cholesky(nW, Imult=8)
n * log(rho) + (2 * c(determinant(update(nChol, nW, 1/rho))$modulus))
-#> [1] 99.8069
-# }
+}
nb7rt <- spdep::cell2nb(7, 7, torus=TRUE)
x <- matrix(sample(rnorm(500*length(nb7rt))), nrow=length(nb7rt))
lw <- spdep::nb2listw(nb7rt)
@@ -179,46 +158,33 @@ Examples
W <- as(lw, "CsparseMatrix")
system.time(e <- invIrM(nb7rt, rho=0.98, method="solve", feasible=NULL) %*% x)
#> user system elapsed
-#> 0.004 0.000 0.004
+#> 0.003 0.000 0.003
system.time(ee <- powerWeights(W, rho=0.98, X=x))
#> Warning: not converged within order iterations
#> user system elapsed
-#> 0.212 0.003 0.224
+#> 0.185 0.004 0.189
str(attr(ee, "internal"))
#> List of 5
-#> $ series: num [1:250] 0.287 0.234 0.201 0.178 0.16 ...
+#> $ series: num [1:250] 0.286 0.233 0.199 0.175 0.157 ...
#> $ order : num 250
#> $ tol : num 4.05e-10
#> $ iter : num 250
#> $ conv : logi FALSE
-all.equal(e, as(ee, "matrix"), check.attributes=FALSE)
-#> [1] "Mean relative difference: 0.0060747"
-# \dontrun{
+all.equal(e, as(ee, "matrix"), check.attributes=FALSE)
+#> [1] "Mean relative difference: 0.00604604"
+if (FALSE) {
system.time(ee <- powerWeights(W, rho=0.9, X=x))
-#> user system elapsed
-#> 0.167 0.002 0.179
system.time(ee <- powerWeights(W, rho=0.98, order=1000, X=x))
-#> user system elapsed
-#> 0.806 0.005 0.856
-all.equal(e, as(ee, "matrix"), check.attributes=FALSE)
-#> [1] TRUE
+all.equal(e, as(ee, "matrix"), check.attributes=FALSE)
nb60rt <- spdep::cell2nb(60, 60, torus=TRUE)
W <- as(spdep::nb2listw(nb60rt), "CsparseMatrix")
set.seed(1)
x <- matrix(rnorm(dim(W)[1]), ncol=1)
system.time(ee <- powerWeights(W, rho=0.3, X=x))
-#> user system elapsed
-#> 0.021 0.000 0.022
str(as(ee, "matrix"))
-#> num [1:3600, 1] -0.383 0.207 -0.731 1.552 0.32 ...
-#> - attr(*, "dimnames")=List of 2
-#> ..$ : chr [1:3600] "1:1" "2:1" "3:1" "4:1" ...
-#> ..$ : NULL
obj <- errorsarlm(as(ee, "matrix")[,1] ~ 1, listw=spdep::nb2listw(nb60rt), method="Matrix")
coefficients(obj)
-#> lambda (Intercept)
-#> 0.30639880 0.01380415
-# }
+}
require("sf", quietly=TRUE)
nydata <- st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData")[1], quiet=TRUE)
-# \dontrun{
+if (FALSE) {
lm0 <- lm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata)
summary(lm0)
-#>
-#> Call:
-#> lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata)
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.7417 -0.3957 -0.0326 0.3353 4.1398
-#>
-#> Coefficients:
-#> Estimate Std. Error t value Pr(>|t|)
-#> (Intercept) -0.51728 0.15856 -3.262 0.00124 **
-#> PEXPOSURE 0.04884 0.03506 1.393 0.16480
-#> PCTAGE65P 3.95089 0.60550 6.525 3.22e-10 ***
-#> PCTOWNHOME -0.56004 0.17031 -3.288 0.00114 **
-#> ---
-#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-#>
-#> Residual standard error: 0.6571 on 277 degrees of freedom
-#> Multiple R-squared: 0.1932, Adjusted R-squared: 0.1844
-#> F-statistic: 22.1 on 3 and 277 DF, p-value: 7.306e-13
-#>
lm0w <- lm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, weights=POP8)
summary(lm0w)
-#>
-#> Call:
-#> lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> weights = POP8)
-#>
-#> Weighted Residuals:
-#> Min 1Q Median 3Q Max
-#> -129.067 -14.714 5.817 25.624 70.723
-#>
-#> Coefficients:
-#> Estimate Std. Error t value Pr(>|t|)
-#> (Intercept) -0.77837 0.14116 -5.514 8.03e-08 ***
-#> PEXPOSURE 0.07626 0.02731 2.792 0.00560 **
-#> PCTAGE65P 3.85656 0.57126 6.751 8.60e-11 ***
-#> PCTOWNHOME -0.39869 0.15305 -2.605 0.00968 **
-#> ---
-#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
-#>
-#> Residual standard error: 33.5 on 277 degrees of freedom
-#> Multiple R-squared: 0.1977, Adjusted R-squared: 0.189
-#> F-statistic: 22.75 on 3 and 277 DF, p-value: 3.382e-13
-#>
-# }
+}
suppressMessages(nyadjmat <- as.matrix(foreign::read.dbf(system.file(
"misc/nyadjwts.dbf", package="spData")[1])[-1]))
suppressMessages(ID <- as.character(names(foreign::read.dbf(system.file(
@@ -386,263 +343,43 @@ Examples
#> [1] TRUE
#require("spdep", quietly=TRUE)
nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY))
+#> Warning: style is M (missing); style should be set to a valid value
listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B")
eigs <- eigenw(listw_NY)
-# \dontrun{
+if (FALSE) {
esar0 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY)
summary(esar0)
-#>
-#> Call:errorsarlm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
-#> data = nydata, listw = listw_NY)
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.56754 -0.38239 -0.02643 0.33109 4.01219
-#>
-#> Type: error
-#> Coefficients: (asymptotic standard errors)
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.618193 0.176784 -3.4969 0.0004707
-#> PEXPOSURE 0.071014 0.042051 1.6888 0.0912635
-#> PCTAGE65P 3.754200 0.624722 6.0094 1.862e-09
-#> PCTOWNHOME -0.419890 0.191329 -2.1946 0.0281930
-#>
-#> Lambda: 0.040487, LR test value: 5.2438, p-value: 0.022026
-#> Asymptotic standard error: 0.016214
-#> z-value: 2.4971, p-value: 0.01252
-#> Wald statistic: 6.2356, p-value: 0.01252
-#>
-#> Log likelihood: -276.1069 for error model
-#> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 564.21, (AIC for lm: 567.46)
-#>
system.time(esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data=nydata, listw=listw_NY, family="SAR", method="eigen",
control=list(pre_eig=eigs)))
-#> user system elapsed
-#> 0.272 0.000 0.274
res <- summary(esar1f)
print(res)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, family = "SAR", method = "eigen", control = list(pre_eig = eigs))
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.56754 -0.38239 -0.02643 0.33109 4.01219
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.618193 0.176784 -3.4969 0.0004707
-#> PEXPOSURE 0.071014 0.042051 1.6888 0.0912635
-#> PCTAGE65P 3.754200 0.624722 6.0094 1.862e-09
-#> PCTOWNHOME -0.419890 0.191329 -2.1946 0.0281930
-#>
-#> Lambda: 0.040487 LR test value: 5.2438 p-value: 0.022026
-#> Numerical Hessian standard error of lambda: 0.017209
-#>
-#> Log likelihood: -276.1069
-#> ML residual variance (sigma squared): 0.41388, (sigma: 0.64333)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 564.21
-#>
coef(res)
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.61819272 0.17678351 -3.496891 4.707136e-04
-#> PEXPOSURE 0.07101384 0.04205063 1.688770 9.126351e-02
-#> PCTAGE65P 3.75419997 0.62472153 6.009397 1.862141e-09
-#> PCTOWNHOME -0.41988961 0.19132936 -2.194591 2.819298e-02
sqrt(diag(res$resvar))
-#> (Intercept) PEXPOSURE PCTAGE65P PCTOWNHOME
-#> 0.17678351 0.04205063 0.62472153 0.19132936
sqrt(diag(esar1f$fit$imat)*esar1f$fit$s2)
-#> (Intercept) PEXPOSURE PCTAGE65P PCTOWNHOME
-#> 0.17678351 0.04205063 0.62472153 0.19132936
sqrt(diag(esar1f$fdHess))
-#> [1] 0.01720868 0.18535631 0.04389387 0.63003835 0.20373413
system.time(esar1M <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data=nydata, listw=listw_NY, family="SAR", method="Matrix"))
-#> user system elapsed
-#> 0.323 0.000 0.325
summary(esar1M)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, family = "SAR", method = "Matrix")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -3.406132 -0.561646 -0.092662 0.474796 5.384405
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.414826 0.102166 -4.0603 4.901e-05
-#> PEXPOSURE 0.015081 0.017772 0.8486 0.3961
-#> PCTAGE65P 5.159749 0.476498 10.8285 < 2.2e-16
-#> PCTOWNHOME -0.892387 0.099241 -8.9921 < 2.2e-16
-#>
-#> Lambda: -0.38889 LR test value: 254.73 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: 0.044857
-#>
-#> Log likelihood: -151.3662
-#> ML residual variance (sigma squared): 0.95111, (sigma: 0.97525)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 314.73
-#>
system.time(esar1M <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data=nydata, listw=listw_NY, family="SAR", method="Matrix",
control=list(super=TRUE)))
-#> user system elapsed
-#> 0.282 0.000 0.283
summary(esar1M)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, family = "SAR", method = "Matrix", control = list(super = TRUE))
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -3.178535 -0.521860 -0.074421 0.414212 5.184465
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.411041 0.104426 -3.9362 8.279e-05
-#> PEXPOSURE 0.015768 0.018366 0.8585 0.3906
-#> PCTAGE65P 5.070130 0.483332 10.4900 < 2.2e-16
-#> PCTOWNHOME -0.880883 0.102093 -8.6282 < 2.2e-16
-#>
-#> Lambda: -0.33146 LR test value: 247.44 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: 0.030659
-#>
-#> Log likelihood: -155.0108
-#> ML residual variance (sigma squared): 0.82436, (sigma: 0.90795)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 322.02
-#>
esar1wf <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="SAR", method="eigen",
control=list(pre_eig=eigs))
summary(esar1wf)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, weights = POP8, family = "SAR", method = "eigen",
-#> control = list(pre_eig = eigs))
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.48488 -0.26823 0.09489 0.46552 4.28343
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.797063 0.144054 -5.5331 3.146e-08
-#> PEXPOSURE 0.080545 0.028334 2.8428 0.004473
-#> PCTAGE65P 3.816731 0.576037 6.6258 3.453e-11
-#> PCTOWNHOME -0.380778 0.156507 -2.4330 0.014975
-#>
-#> Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764
-#> Numerical Hessian standard error of lambda: 0.016466
-#>
-#> Log likelihood: -251.6017
-#> ML residual variance (sigma squared): 1104.1, (sigma: 33.229)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 515.2
-#>
system.time(esar1wM <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data=nydata, listw=listw_NY, weights=POP8, family="SAR", method="Matrix"))
-#> user system elapsed
-#> 0.320 0.000 0.322
summary(esar1wM)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, weights = POP8, family = "SAR", method = "Matrix")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.561100 -0.374524 0.057405 0.591094 5.700142
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.578546 0.090006 -6.4279 1.294e-10
-#> PEXPOSURE 0.035402 0.013959 2.5361 0.01121
-#> PCTAGE65P 4.651137 0.421285 11.0404 < 2.2e-16
-#> PCTOWNHOME -0.666898 0.091443 -7.2931 3.029e-13
-#>
-#> Lambda: -0.34423 LR test value: 264.24 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: 0.036776
-#>
-#> Log likelihood: -119.6468
-#> ML residual variance (sigma squared): 2129.1, (sigma: 46.142)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 251.29
-#>
esar1wlu <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="SAR", method="LU")
summary(esar1wlu)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, weights = POP8, family = "SAR", method = "LU")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.48488 -0.26823 0.09489 0.46552 4.28343
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.797063 0.144054 -5.5331 3.146e-08
-#> PEXPOSURE 0.080545 0.028334 2.8428 0.004473
-#> PCTAGE65P 3.816731 0.576037 6.6258 3.453e-11
-#> PCTOWNHOME -0.380778 0.156507 -2.4330 0.014975
-#>
-#> Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764
-#> Numerical Hessian standard error of lambda: 0.016522
-#>
-#> Log likelihood: -251.6017
-#> ML residual variance (sigma squared): 1104.1, (sigma: 33.229)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 515.2
-#>
esar1wch <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="SAR", method="Chebyshev")
summary(esar1wch)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, weights = POP8, family = "SAR", method = "Chebyshev")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -4.39831 -0.86303 0.12117 1.09320 8.78570
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.538555 0.087573 -6.1498 7.76e-10
-#> PEXPOSURE 0.029782 0.013074 2.2780 0.02273
-#> PCTAGE65P 4.896346 0.419359 11.6758 < 2.2e-16
-#> PCTOWNHOME -0.737829 0.087569 -8.4257 < 2.2e-16
-#>
-#> Lambda: -1 LR test value: 236970 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: NaN
-#>
-#> Log likelihood: 118232.5
-#> ML residual variance (sigma squared): 9336.4, (sigma: 96.625)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: -236450
-#>
-# }
+}
ecar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, family="CAR", method="eigen",
control=list(pre_eig=eigs))
@@ -664,7 +401,7 @@ Examples
#> PCTOWNHOME -0.382789 0.195564 -1.9574 0.0503053
#>
#> Lambda: 0.084123 LR test value: 5.8009 p-value: 0.016018
-#> Numerical Hessian standard error of lambda: 0.030868
+#> Numerical Hessian standard error of lambda: 0.030872
#>
#> Log likelihood: -275.8283
#> ML residual variance (sigma squared): 0.40758, (sigma: 0.63842)
@@ -672,38 +409,11 @@ Examples
#> Number of parameters estimated: 6
#> AIC: 563.66
#>
-# \dontrun{
+if (FALSE) {
system.time(ecar1M <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data=nydata, listw=listw_NY, family="CAR", method="Matrix"))
-#> user system elapsed
-#> 0.349 0.000 0.350
summary(ecar1M)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, family = "CAR", method = "Matrix")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -3.449951 -0.633777 -0.072436 0.550248 6.039594
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.402788 0.112612 -3.5768 0.0003479
-#> PEXPOSURE 0.020131 0.020783 0.9686 0.3327222
-#> PCTAGE65P 4.632644 0.500526 9.2555 < 2.2e-16
-#> PCTOWNHOME -0.812981 0.112867 -7.2030 5.891e-13
-#>
-#> Lambda: -0.5 LR test value: 228.48 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: NaN
-#>
-#> Log likelihood: -164.4897
-#> ML residual variance (sigma squared): 0.50386, (sigma: 0.70983)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 340.98
-#>
-# }
+}
ecar1wf <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
listw=listw_NY, weights=POP8, family="CAR", method="eigen",
control=list(pre_eig=eigs))
@@ -726,47 +436,20 @@ Examples
#> PCTOWNHOME -0.386820 0.157436 -2.4570 0.014010
#>
#> Lambda: 0.022419 LR test value: 0.38785 p-value: 0.53343
-#> Numerical Hessian standard error of lambda: 0.038543
+#> Numerical Hessian standard error of lambda: 0.038977
#>
#> Log likelihood: -251.5711
#> ML residual variance (sigma squared): 1102.9, (sigma: 33.21)
#> Number of observations: 281
#> Number of parameters estimated: 6
-#> AIC: 515.14
+#> AIC: NA (not available for weighted model)
#>
-# \dontrun{
+if (FALSE) {
system.time(ecar1wM <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data=nydata, listw=listw_NY, weights=POP8, family="CAR", method="Matrix"))
-#> user system elapsed
-#> 0.342 0.000 0.344
summary(ecar1wM)
-#>
-#> Call:
-#> spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = nydata,
-#> listw = listw_NY, weights = POP8, family = "CAR", method = "Matrix")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.98144 -0.15716 0.37342 1.08857 7.10495
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) -0.714952 0.093905 -7.6136 2.665e-14
-#> PEXPOSURE 0.041467 0.015279 2.7139 0.00665
-#> PCTAGE65P 4.149207 0.425446 9.7526 < 2.2e-16
-#> PCTOWNHOME -0.478396 0.097030 -4.9304 8.205e-07
-#>
-#> Lambda: -0.5 LR test value: 262.65 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: NaN
-#>
-#> Log likelihood: -120.4395
-#> ML residual variance (sigma squared): 1159.2, (sigma: 34.046)
-#> Number of observations: 281
-#> Number of parameters estimated: 6
-#> AIC: 252.88
-#>
-# }
-# \dontrun{
+}
+if (FALSE) {
require("sf", quietly=TRUE)
nc.sids <- st_read(system.file("shapes/sids.shp", package="spData")[1], quiet=TRUE)
ft.SID74 <- sqrt(1000)*(sqrt(nc.sids$SID74/nc.sids$BIR74) +
@@ -777,21 +460,18 @@ Examples
sids.nhbr30.dist <- spdep::nbdists(sids.nhbr30, cbind(nc.sids$east, nc.sids$north))
sids.nhbr <- spdep::listw2sn(spdep::nb2listw(sids.nhbr30,
glist=sids.nhbr30.dist, style="B", zero.policy=TRUE))
-#> Warning: zero sum general weights
dij <- sids.nhbr[,3]
n <- nc.sids$BIR74
el1 <- min(dij)/dij
el2 <- sqrt(n[sids.nhbr$to]/n[sids.nhbr$from])
sids.nhbr$weights <- el1*el2
sids.nhbr.listw <- spdep::sn2listw(sids.nhbr)
-#> Warning: 56, 87 are not origins
both <- factor(paste(nc.sids$L_id, nc.sids$M_id, sep=":"))
ft.NWBIR74 <- sqrt(1000)*(sqrt(nc.sids$NWBIR74/nc.sids$BIR74) +
sqrt((nc.sids$NWBIR74+1)/nc.sids$BIR74))
mdata <- data.frame(both, ft.NWBIR74, ft.SID74, BIR74=nc.sids$BIR74)
outl <- which.max(rstandard(lm_nc))
as.character(nc.sids$NAME[outl])
-#> [1] "Anson"
mdata.4 <- mdata[-outl,]
W <- spdep::listw2mat(sids.nhbr.listw)
W.4 <- W[-outl, -outl]
@@ -799,395 +479,50 @@ Examples
esarI <- errorsarlm(ft.SID74 ~ 1, data=mdata, listw=sids.nhbr.listw,
zero.policy=TRUE)
summary(esarI)
-#>
-#> Call:errorsarlm(formula = ft.SID74 ~ 1, data = mdata, listw = sids.nhbr.listw,
-#> zero.policy = TRUE)
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.887117 -0.636573 -0.043429 0.448767 3.406724
-#>
-#> Type: error
-#> Regions with no neighbours included:
-#> 56 87
-#> Coefficients: (asymptotic standard errors)
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 2.97463 0.13011 22.862 < 2.2e-16
-#>
-#> Lambda: 0.66864, LR test value: 10.214, p-value: 0.0013939
-#> Asymptotic standard error: 0.11473
-#> z-value: 5.8278, p-value: 5.6147e-09
-#> Wald statistic: 33.964, p-value: 5.6147e-09
-#>
-#> Log likelihood: -133.8616 for error model
-#> ML residual variance (sigma squared): 0.81932, (sigma: 0.90516)
-#> Number of observations: 100
-#> Number of parameters estimated: 3
-#> AIC: 273.72, (AIC for lm: 281.94)
-#>
esarIa <- spautolm(ft.SID74 ~ 1, data=mdata, listw=sids.nhbr.listw,
family="SAR")
summary(esarIa)
-#>
-#> Call: spautolm(formula = ft.SID74 ~ 1, data = mdata, listw = sids.nhbr.listw,
-#> family = "SAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.887117 -0.636573 -0.043429 0.448767 3.406724
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 2.97463 0.13011 22.862 < 2.2e-16
-#>
-#> Lambda: 0.66864 LR test value: 10.214 p-value: 0.0013939
-#> Numerical Hessian standard error of lambda: 0.16506
-#>
-#> Log likelihood: -133.8616
-#> ML residual variance (sigma squared): 0.81932, (sigma: 0.90516)
-#> Number of observations: 100
-#> Number of parameters estimated: 3
-#> AIC: 273.72
-#>
esarIV <- errorsarlm(ft.SID74 ~ ft.NWBIR74, data=mdata, listw=sids.nhbr.listw,
zero.policy=TRUE)
summary(esarIV)
-#>
-#> Call:
-#> errorsarlm(formula = ft.SID74 ~ ft.NWBIR74, data = mdata, listw = sids.nhbr.listw,
-#> zero.policy = TRUE)
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.123648 -0.573163 0.017859 0.468022 2.693604
-#>
-#> Type: error
-#> Regions with no neighbours included:
-#> 56 87
-#> Coefficients: (asymptotic standard errors)
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.549443 0.219230 7.0677 1.576e-12
-#> ft.NWBIR74 0.041974 0.006171 6.8018 1.033e-11
-#>
-#> Lambda: 0.18465, LR test value: 0.50496, p-value: 0.47733
-#> Asymptotic standard error: 0.20648
-#> z-value: 0.89424, p-value: 0.37119
-#> Wald statistic: 0.79967, p-value: 0.37119
-#>
-#> Log likelihood: -117.7464 for error model
-#> ML residual variance (sigma squared): 0.61546, (sigma: 0.78451)
-#> Number of observations: 100
-#> Number of parameters estimated: 4
-#> AIC: 243.49, (AIC for lm: 242)
-#>
esarIVa <- spautolm(ft.SID74 ~ ft.NWBIR74, data=mdata, listw=sids.nhbr.listw,
family="SAR")
summary(esarIVa)
-#>
-#> Call:
-#> spautolm(formula = ft.SID74 ~ ft.NWBIR74, data = mdata, listw = sids.nhbr.listw,
-#> family = "SAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.123648 -0.573163 0.017859 0.468022 2.693604
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.549443 0.219230 7.0677 1.576e-12
-#> ft.NWBIR74 0.041974 0.006171 6.8018 1.033e-11
-#>
-#> Lambda: 0.18465 LR test value: 0.50496 p-value: 0.47733
-#> Numerical Hessian standard error of lambda: 0.25591
-#>
-#> Log likelihood: -117.7464
-#> ML residual variance (sigma squared): 0.61546, (sigma: 0.78451)
-#> Number of observations: 100
-#> Number of parameters estimated: 4
-#> AIC: 243.49
-#>
esarIaw <- spautolm(ft.SID74 ~ 1, data=mdata, listw=sids.nhbr.listw,
weights=BIR74, family="SAR")
summary(esarIaw)
-#>
-#> Call: spautolm(formula = ft.SID74 ~ 1, data = mdata, listw = sids.nhbr.listw,
-#> weights = BIR74, family = "SAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.867485 -0.568644 0.019717 0.502197 3.498013
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 2.852052 0.090271 31.594 < 2.2e-16
-#>
-#> Lambda: 0.7338 LR test value: 12.917 p-value: 0.00032554
-#> Numerical Hessian standard error of lambda: 0.13887
-#>
-#> Log likelihood: -130.0975
-#> ML residual variance (sigma squared): 1539.4, (sigma: 39.236)
-#> Number of observations: 100
-#> Number of parameters estimated: 3
-#> AIC: 266.19
-#>
esarIIaw <- spautolm(ft.SID74 ~ both - 1, data=mdata, listw=sids.nhbr.listw,
weights=BIR74, family="SAR")
summary(esarIIaw)
-#>
-#> Call:
-#> spautolm(formula = ft.SID74 ~ both - 1, data = mdata, listw = sids.nhbr.listw,
-#> weights = BIR74, family = "SAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.590809 -0.432976 0.016736 0.357284 3.536718
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> both1:2 2.05545 0.22184 9.2654 < 2.2e-16
-#> both1:3 2.87260 0.16181 17.7531 < 2.2e-16
-#> both1:4 4.16365 0.34330 12.1283 < 2.2e-16
-#> both2:1 2.47255 0.29757 8.3090 < 2.2e-16
-#> both2:2 2.15307 0.21172 10.1692 < 2.2e-16
-#> both2:3 2.64235 0.17296 15.2770 < 2.2e-16
-#> both2:4 3.26604 0.28287 11.5459 < 2.2e-16
-#> both3:1 3.11277 0.34166 9.1107 < 2.2e-16
-#> both3:2 2.76541 0.15667 17.6508 < 2.2e-16
-#> both3:3 2.86582 0.18593 15.4134 < 2.2e-16
-#> both3:4 3.18142 0.21617 14.7169 < 2.2e-16
-#> both4:3 3.69333 0.23348 15.8188 < 2.2e-16
-#>
-#> Lambda: 0.32136 LR test value: 1.4004 p-value: 0.23666
-#> Numerical Hessian standard error of lambda: 0.25465
-#>
-#> Log likelihood: -109.8922
-#> ML residual variance (sigma squared): 1071.6, (sigma: 32.735)
-#> Number of observations: 100
-#> Number of parameters estimated: 14
-#> AIC: 247.78
-#>
esarIVaw <- spautolm(ft.SID74 ~ ft.NWBIR74, data=mdata,
listw=sids.nhbr.listw, weights=BIR74, family="SAR")
summary(esarIVaw)
-#>
-#> Call:
-#> spautolm(formula = ft.SID74 ~ ft.NWBIR74, data = mdata, listw = sids.nhbr.listw,
-#> weights = BIR74, family = "SAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.00956 -0.45229 0.12547 0.55952 2.92223
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.5769279 0.2501334 6.3043 2.894e-10
-#> ft.NWBIR74 0.0368573 0.0069413 5.3099 1.097e-07
-#>
-#> Lambda: 0.3839 LR test value: 1.9983 p-value: 0.15747
-#> Numerical Hessian standard error of lambda: 0.25769
-#>
-#> Log likelihood: -119.5648
-#> ML residual variance (sigma squared): 1295.8, (sigma: 35.997)
-#> Number of observations: 100
-#> Number of parameters estimated: 4
-#> AIC: 247.13
-#>
ecarIaw <- spautolm(ft.SID74 ~ 1, data=mdata.4, listw=sids.nhbr.listw.4,
weights=BIR74, family="CAR")
-#> Warning: Non-symmetric spatial weights in CAR model
summary(ecarIaw)
-#>
-#> Call:
-#> spautolm(formula = ft.SID74 ~ 1, data = mdata.4, listw = sids.nhbr.listw.4,
-#> weights = BIR74, family = "CAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.009350 -0.638915 -0.060761 0.428526 2.019409
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 2.942864 0.095304 30.879 < 2.2e-16
-#>
-#> Lambda: 0.86832 LR test value: 23.003 p-value: 1.6172e-06
-#> Numerical Hessian standard error of lambda: 0.048102
-#>
-#> Log likelihood: -118.7564
-#> ML residual variance (sigma squared): 1264, (sigma: 35.553)
-#> Number of observations: 99
-#> Number of parameters estimated: 3
-#> AIC: 243.51
-#>
ecarIIaw <- spautolm(ft.SID74 ~ both - 1, data=mdata.4,
listw=sids.nhbr.listw.4, weights=BIR74, family="CAR")
-#> Warning: Non-symmetric spatial weights in CAR model
-#> Warning: NaNs produced
summary(ecarIIaw)
-#>
-#> Call:
-#> spautolm(formula = ft.SID74 ~ both - 1, data = mdata.4, listw = sids.nhbr.listw.4,
-#> weights = BIR74, family = "CAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -2.564067 -0.461531 -0.020982 0.384458 2.054255
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> both1:2 2.06282 0.20065 10.2806 < 2.2e-16
-#> both1:3 2.91982 0.14171 20.6048 < 2.2e-16
-#> both1:4 4.12159 0.30076 13.7037 < 2.2e-16
-#> both2:1 2.58281 0.27014 9.5611 < 2.2e-16
-#> both2:2 2.17549 0.18265 11.9104 < 2.2e-16
-#> both2:3 2.67030 0.15355 17.3910 < 2.2e-16
-#> both2:4 3.10806 0.24748 12.5588 < 2.2e-16
-#> both3:1 2.93237 0.30007 9.7724 < 2.2e-16
-#> both3:2 2.65317 0.14139 18.7646 < 2.2e-16
-#> both3:3 2.91685 0.17134 17.0234 < 2.2e-16
-#> both3:4 3.20447 0.20402 15.7063 < 2.2e-16
-#> both4:3 3.80672 0.20831 18.2742 < 2.2e-16
-#>
-#> Lambda: 0.22163 LR test value: 1.3827 p-value: 0.23964
-#> Numerical Hessian standard error of lambda: NaN
-#>
-#> Log likelihood: -99.2181
-#> ML residual variance (sigma squared): 890.66, (sigma: 29.844)
-#> Number of observations: 99
-#> Number of parameters estimated: 14
-#> AIC: 226.44
-#>
ecarIVaw <- spautolm(ft.SID74 ~ ft.NWBIR74, data=mdata.4,
listw=sids.nhbr.listw.4, weights=BIR74, family="CAR")
-#> Warning: Non-symmetric spatial weights in CAR model
summary(ecarIVaw)
-#>
-#> Call:
-#> spautolm(formula = ft.SID74 ~ ft.NWBIR74, data = mdata.4, listw = sids.nhbr.listw.4,
-#> weights = BIR74, family = "CAR")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -1.99259 -0.44794 0.15464 0.60748 1.95751
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 1.434705 0.225521 6.3618 1.995e-10
-#> ft.NWBIR74 0.040903 0.006299 6.4936 8.382e-11
-#>
-#> Lambda: 0.22724 LR test value: 1.1936 p-value: 0.2746
-#> Numerical Hessian standard error of lambda: 0.55494
-#>
-#> Log likelihood: -114.0196
-#> ML residual variance (sigma squared): 1201, (sigma: 34.655)
-#> Number of observations: 99
-#> Number of parameters estimated: 4
-#> AIC: 236.04
-#>
nc.sids$fitIV <- append(fitted.values(ecarIVaw), NA, outl-1)
plot(nc.sids[,"fitIV"], nbreaks=12) # Cressie 1993, p. 565
-
-# }
-# \dontrun{
+}
+if (FALSE) {
data(oldcol, package="spdep")
COL.errW.eig <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD,
spdep::nb2listw(COL.nb, style="W"))
summary(COL.errW.eig)
-#>
-#> Call:
-#> errorsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = spdep::nb2listw(COL.nb,
-#> style = "W"))
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -34.81174 -6.44031 -0.72142 7.61476 23.33626
-#>
-#> Type: error
-#> Coefficients: (asymptotic standard errors)
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 59.893219 5.366163 11.1613 < 2.2e-16
-#> INC -0.941312 0.330569 -2.8476 0.0044057
-#> HOVAL -0.302250 0.090476 -3.3407 0.0008358
-#>
-#> Lambda: 0.56179, LR test value: 7.9935, p-value: 0.0046945
-#> Asymptotic standard error: 0.13387
-#> z-value: 4.1966, p-value: 2.7098e-05
-#> Wald statistic: 17.611, p-value: 2.7098e-05
-#>
-#> Log likelihood: -183.3805 for error model
-#> ML residual variance (sigma squared): 95.575, (sigma: 9.7762)
-#> Number of observations: 49
-#> Number of parameters estimated: 5
-#> AIC: 376.76, (AIC for lm: 382.75)
-#>
COL.errW.sar <- spautolm(CRIME ~ INC + HOVAL, data=COL.OLD,
spdep::nb2listw(COL.nb, style="W"))
summary(COL.errW.sar)
-#>
-#> Call:
-#> spautolm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = spdep::nb2listw(COL.nb,
-#> style = "W"))
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -34.81174 -6.44031 -0.72142 7.61476 23.33626
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 59.893219 5.366163 11.1613 < 2.2e-16
-#> INC -0.941312 0.330569 -2.8476 0.0044057
-#> HOVAL -0.302250 0.090476 -3.3407 0.0008358
-#>
-#> Lambda: 0.56179 LR test value: 7.9935 p-value: 0.0046945
-#> Numerical Hessian standard error of lambda: 0.15242
-#>
-#> Log likelihood: -183.3805
-#> ML residual variance (sigma squared): 95.575, (sigma: 9.7762)
-#> Number of observations: 49
-#> Number of parameters estimated: 5
-#> AIC: 376.76
-#>
data(boston, package="spData")
gp1 <- spautolm(log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2)
+ I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B + log(LSTAT),
data=boston.c, spdep::nb2listw(boston.soi), family="SMA")
summary(gp1)
-#>
-#> Call: spautolm(formula = log(CMEDV) ~ CRIM + ZN + INDUS + CHAS + I(NOX^2) +
-#> I(RM^2) + AGE + log(DIS) + log(RAD) + TAX + PTRATIO + B +
-#> log(LSTAT), data = boston.c, listw = spdep::nb2listw(boston.soi),
-#> family = "SMA")
-#>
-#> Residuals:
-#> Min 1Q Median 3Q Max
-#> -0.5847694 -0.0713881 0.0012284 0.0827517 0.6071219
-#>
-#> Coefficients:
-#> Estimate Std. Error z value Pr(>|z|)
-#> (Intercept) 4.28501607 0.15367176 27.8842 < 2.2e-16
-#> CRIM -0.00718807 0.00106298 -6.7622 1.359e-11
-#> ZN 0.00023008 0.00051897 0.4433 0.6575185
-#> INDUS 0.00047300 0.00263339 0.1796 0.8574551
-#> CHAS1 0.01020698 0.02872047 0.3554 0.7222970
-#> I(NOX^2) -0.44885530 0.13675913 -3.2821 0.0010304
-#> I(RM^2) 0.00638094 0.00110330 5.7835 7.316e-09
-#> AGE -0.00043973 0.00051336 -0.8566 0.3916862
-#> log(DIS) -0.15650578 0.03856337 -4.0584 4.941e-05
-#> log(RAD) 0.07583760 0.02016468 3.7609 0.0001693
-#> TAX -0.00049364 0.00012162 -4.0588 4.933e-05
-#> PTRATIO -0.02494959 0.00538791 -4.6307 3.645e-06
-#> B 0.00048517 0.00010944 4.4334 9.277e-06
-#> log(LSTAT) -0.32961379 0.02353891 -14.0029 < 2.2e-16
-#>
-#> Lambda: 0.61991 LR test value: 144.28 p-value: < 2.22e-16
-#> Numerical Hessian standard error of lambda: 0.044359
-#>
-#> Log likelihood: 229.1208
-#> ML residual variance (sigma squared): 0.02596, (sigma: 0.16112)
-#> Number of observations: 506
-#> Number of parameters estimated: 16
-#> AIC: -426.24
-#>
-# }
+}