diff --git a/DESCRIPTION b/DESCRIPTION index f249883..982e91b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: spatialreg -Version: 1.3-2 -Date: 2024-02-06 +Version: 1.3-3 +Date: 2024-04-29 Title: Spatial Regression Analysis Encoding: UTF-8 Authors@R: c(person("Roger", "Bivand", role = c("cre", "aut"), email = "Roger.Bivand@nhh.no", comment=c(ORCID="0000-0003-2392-6140")), @@ -29,7 +29,7 @@ Depends: R (>= 3.3.0), spData, Matrix, sf Imports: spdep (>= 1.3-1), coda, methods, MASS, boot, splines, LearnBayes, nlme, multcomp Suggests: parallel, RSpectra, tmap, foreign, spam, knitr, lmtest, expm, - sandwich, rmarkdown, igraph, tinytest + sandwich, rmarkdown, igraph (>= 2.0.0), tinytest Description: A collection of all the estimation functions for spatial cross-sectional models (on lattice/areal data using spatial weights matrices) contained up to now in 'spdep'. These model fitting functions include maximum likelihood methods for cross-sectional models proposed by 'Cliff' and 'Ord' (1973, ISBN:0850860369) and (1981, ISBN:0850860814), fitting methods initially described by 'Ord' (1975) . The models are further described by 'Anselin' (1988) . Spatial two stage least squares and spatial general method of moment models initially proposed by 'Kelejian' and 'Prucha' (1998) and (1999) are provided. Impact methods and MCMC fitting methods proposed by 'LeSage' and 'Pace' (2009) are implemented for the family of cross-sectional spatial regression models. Methods for fitting the log determinant term in maximum likelihood and MCMC fitting are compared by 'Bivand et al.' (2013) , and model fitting methods by 'Bivand' and 'Piras' (2015) ; both of these articles include extensive lists of references. A recent review is provided by 'Bivand', 'Millo' and 'Piras' (2021) . 'spatialreg' >= 1.1-* corresponded to 'spdep' >= 1.1-1, in which the model fitting functions were deprecated and passed through to 'spatialreg', but masked those in 'spatialreg'. From versions 1.2-*, the functions have been made defunct in 'spdep'. License: GPL-2 URL: https://github.com/r-spatial/spatialreg/, https://r-spatial.github.io/spatialreg/ diff --git a/NEWS.md b/NEWS.md index a41eae0..7ec3a50 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,11 @@ # Version 1.3-2 (development) +* add `sqrt=TRUE` in calls to Matrix `determinant` methods in `matrix_ldet` + +* add `igraph (>= 2.0.0)` in DESCRIPTION for re-named `igraph` functions + +# Version 1.3-2 (2024-02-06) + * pass through SlX formula in call * re-corrected #19 because the fitted model weights component may be NULL diff --git a/R/jacobian.R b/R/jacobian.R index 8dd1400..bc3c97d 100644 --- a/R/jacobian.R +++ b/R/jacobian.R @@ -376,10 +376,10 @@ Matrix_ldet <- function(coef, env, which=1) { .f <- get(".f", envir=env) Jacobian <- ifelse(coef > b, n * log(coef) + - (.f * c(determinant(update(nChol, nW, 1/coef))$modulus)), + (.f * c(determinant(update(nChol, nW, 1/coef), sqrt=TRUE)$modulus)), ifelse(coef < a, n* log(-(coef)) + - (.f * c(determinant(update(pChol, csrw, 1/(-coef)))$modulus)), - 0.0)) + (.f * c(determinant(update(pChol, csrw, 1/(-coef)), + sqrt=TRUE)$modulus)), 0.0)) Jacobian } diff --git a/docs/404.html b/docs/404.html index 8faea44..1d5dc7c 100644 --- a/docs/404.html +++ b/docs/404.html @@ -6,7 +6,7 @@ Page not found (404) • spatialreg - + @@ -32,7 +32,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -49,14 +49,14 @@
  • @@ -106,7 +106,7 @@

    Page not found (404)

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/articles/SpatialFiltering.html b/docs/articles/SpatialFiltering.html index 484e162..1f4e6bd 100644 --- a/docs/articles/SpatialFiltering.html +++ b/docs/articles/SpatialFiltering.html @@ -6,7 +6,7 @@ Moran Eigenvectors • spatialreg - + @@ -33,7 +33,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -50,14 +50,14 @@
  • @@ -325,7 +325,7 @@

    References

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/articles/index.html b/docs/articles/index.html index 2053313..c9de137 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -1,5 +1,5 @@ -Articles • spatialregArticles • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@

  • Changelog @@ -66,12 +66,12 @@

    Articles

    All vignettes

    -
    Moran Eigenvectors
    -
    -
    Spatial weights objects as sparse matrices and graphs
    +
    Spatial weights objects as sparse matrices and graphs
    North Carolina SIDS data set (models)
    +
    Moran Eigenvectors
    +
    @@ -82,7 +82,7 @@

    All vignettes

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/articles/nb_igraph.html b/docs/articles/nb_igraph.html index c560d81..7ab39dd 100644 --- a/docs/articles/nb_igraph.html +++ b/docs/articles/nb_igraph.html @@ -6,7 +6,7 @@ Spatial weights objects as sparse matrices and graphs • spatialreg - + @@ -33,7 +33,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -50,14 +50,14 @@
  • @@ -162,7 +162,7 @@ 

    ## Loading required package: spData
    ## Loading required package: Matrix
    ## Loading required package: sf
    -
    ## Linking to GEOS 3.12.1, GDAL 3.8.3, PROJ 9.3.1; sf_use_s2() is TRUE
    +
    ## Linking to GEOS 3.12.1, GDAL 3.9.0beta1, PROJ 9.4.0; sf_use_s2() is TRUE

    Getting some data @@ -181,9 +181,9 @@

    Getting some datasf_extSoftVersion() }

    ##           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
-##       "3.12.1"        "3.8.3"        "9.3.1"         "true"         "true" 
+##       "3.12.1"   "3.9.0beta1"        "9.4.0"         "true"         "true" 
 ##           PROJ 
-##        "9.3.1"
    +## "9.4.0" 
     library(sf)
 columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1])
    @@ -311,25 +311,14 @@

    Symmetric sparse matrices"M" to indicate conversion from a matrix. The neighbour links are retreived correctly, as are the IDs:

    -nb_B1 <- spdep::mat2listw(as(as(B, "TsparseMatrix"), "CsparseMatrix"))
    -
    ## 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"
    -
    +nb_B1 <- spdep::mat2listw(as(as(B, "TsparseMatrix"), "CsparseMatrix"),
+    style="B", zero.policy=TRUE)
+nb_B1$style
    +
    ## [1] "B"
    +
     all.equal(nb_B1$neighbours, col2, check.attributes=FALSE)
    ## [1] TRUE
    -
    +
     all.equal(attr(nb_B1$neighbours, "region.id"), attr(nb_B$neighbours, "region.id"))
    
 
    ## [1] TRUE
    @@ -342,7 +331,7 @@

    Log deter 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
@@ -357,14 +346,14 @@ 
    Log deter
 "dsTMatrix", from a "ddiMatrix". The value of
 the log determinant follows, calling a sparse Cholesky decomposition
 internally for suitable input matrices.
    
-
    +
     n <- nrow(B)
 I <- Diagonal(n)
 class(I - rho * B)
    
 
    ## [1] "dgCMatrix"
 ## attr(,"package")
 ## [1] "Matrix"
    
-
    +
     c(determinant(I - rho * B, logarithm=TRUE)$modulus)
    
 
    ## [1] -1.44787
    
 
    The computation of a sparse Cholesky decomposition for each value of
@@ -372,13 +361,10 @@ 
    Log deter
 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
    
+nChol <- Cholesky(nW, Imult=8)
+n * log(rho) + (2 * c(determinant(update(nChol, nW, 1/rho), sqrt=TRUE)$modulus))
    
 
    ## [1] 15.40218
    
 
    
 
    
@@ -387,7 +373,7 @@ 
    Asymmetric sparse matricesThe 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):
    
-
    +
     nb_W <- spdep::nb2listw(col2, style="W", zero.policy=TRUE)
 W <- as(nb_W, "CsparseMatrix")
 str(W)
    
@@ -400,7 +386,7 @@ 
    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()
    
-
    +
     all(W == t(W))
    
 
    ## [1] FALSE
    
 
    The lag method for listw objects is often
@@ -409,17 +395,17 @@ 
    Asymmetric sparse matriceszero.policy to
 TRUE, the spatial lag of entity 21, which has no
 neighbours, is zero by construction:
    
-
    +
     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
    
-
    +
     plot(x, r1, ylim=c(0,1))
    
 
    
-
    +
     c(x[21], r1[21])
    
 
    ## [1] 0.9347052 0.0000000
    
 
    
@@ -433,21 +419,21 @@ 
    Log dete
 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:
    
-
    +
     rho <- 0.5
 sum(log(1 - rho * eigenw(nb_W)))
    
 
    ## [1] -1.594376
    
-
    +
     class(I - rho * W)
    
 
    ## [1] "dgCMatrix"
 ## attr(,"package")
 ## [1] "Matrix"
    
-
    +
     c(determinant(I - rho * W, logarithm=TRUE)$modulus)
    
 
    ## [1] -1.594376
    
 
    We can show the internal workings of the method as:
    
-
    -LU <- lu(I - rho * W)
+
    +LU <- lu(I - rho * W)
 sum(log(abs(diag(slot(LU, "U")))))
    
 
    ## [1] -1.594376
    
 
    
@@ -460,7 +446,7 @@ 
    Log determinants: symmetric by
 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:
    
-
    +
     d <- attr(nb_W$weights, "comp")$d
 all.equal(d, spdep::card(col2))
    
 
    ## [1] TRUE
    
@@ -468,24 +454,24 @@ 
    Log determinants: symmetric by
 (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:
    
-
    +
     dW <- Diagonal(n, d) %*% W
 all(dW == t(dW))
    
 
    ## [1] TRUE
    
-
    +
     isd <- Diagonal(n, 1/sqrt(d))
 isd[21,21]
    
 
    ## [1] Inf
    
-
    +
     Ws <- as(isd %*% dW %*% isd, "symmetricMatrix")
-rowSums(Ws)[21]
    
+rowSums(Ws)[21]
    
 
    ## [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
@@ -506,19 +492,19 @@ 
    Using eigsn is somewhat larger, use may be made of the
 eigs function in RSpectra:
    
-
    +
     1/range(spatialreg::eigenw(nb_B))
    
 
    ## [1] -0.3212551  0.1638329
    
-
    +
     if (!require("RSpectra", quietly=TRUE)) dothis <- FALSE
    
-
    +
     1/c(eigs(B, k=1, which="SR")$values, eigs(B, k=1, which="LR")$values)
    
 
    ## [1] -0.3212551  0.1638329
    
 
    In this case, the results are trivial with small k.
    
-
    +
     1/range(eigenw(nb_W))
    
 
    ## [1] -1.544645  1.000000
    
-
    +
     1/Re(c(eigs(W, k=1, which="SR")$values, eigs(W, k=1, which="LR")$values))
    
 
    ## [1] -1.544645  1.000000
    
 
    Using row-standardisation has the nice feature of setting the upper
@@ -534,15 +520,15 @@ 
    Converting from sym
 
    
 
    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"
    
-
    +
 
    ## 10824 bytes
    
-
    +
     if (!require("igraph", quietly=FALSE)) dothis <- FALSE
    
 
    ## Loading required package: igraph
    
 
    ## 
@@ -553,17 +539,11 @@ 
    Converting from sym
 
    ## The following object is masked from 'package:base':
 ## 
 ##     union
    
-
    -g1 <- graph.adjacency(B, mode="undirected")
    
-
    ## Warning: `graph.adjacency()` was deprecated in igraph 2.0.0.
-## ℹ Please use `graph_from_adjacency_matrix()` instead.
-## This warning is displayed once every 8 hours.
-## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
-## generated.
    
-
    -class(g1)
    
+
    +g1 <- graph_from_adjacency_matrix(B, mode="undirected")
+class(g1)
    
 
    ## [1] "igraph"
    
-
    +
 
    ## 6544 bytes
    
 
    
@@ -574,15 +554,23 @@ 
    Converting from gra
 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)
-if (!inherits(B1, "try-error")) print(class(B1))
-if (!inherits(B1, "try-error")) print(object.size(B1))
-if (!inherits(B1, "try-error")) print(all.equal(B, as(B1, "CsparseMatrix")))
-options("warn"=ow)
    
+B1 <- try(as_adjacency_matrix(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
    
+
    +if (!inherits(B1, "try-error")) print(all.equal(B, as(B1, "CsparseMatrix")))
    
+
    ## [1] TRUE
    
+
    +options("warn"=ow)
    
 
    
 
    
 
    Graph components in spdep
@@ -592,7 +580,7 @@ 
    Graph components in spdepn.comp.nb from the early days of the package. It is useful
 to know whether an nb object is divided up into separate
 subgraphs, and which entities are members of which such subgraph.
    
-
    +
     res <- spdep::n.comp.nb(col2)
 table(res$comp.id)
    
 
    ## 
@@ -604,27 +592,21 @@ 
    Graph components in igraph<
 
    
 
    The same result can be obtained using the clusters
 function in igraph:
    
-
    -c1 <- clusters(g1)
    
-
    ## Warning: `clusters()` was deprecated in igraph 2.0.0.
-## ℹ Please use `components()` instead.
-## This warning is displayed once every 8 hours.
-## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
-## generated.
    
-
    -c1$no == res$nc
    
+
    +c1 <- components(g1)
+c1$no == res$nc
    
 
    ## [1] TRUE
    
-
    +
     all.equal(c1$membership, res$comp.id)
    
 
    ## [1] "names for target but not for current"
    
-
    +
     all.equal(c1$csize, c(table(res$comp.id)), check.attributes=FALSE)
    
 
    ## [1] TRUE
    
 
    The same holds for the row-standardised variant:
    
-
    +
     W <- as(spdep::nb2listw(col2, style="W", zero.policy=TRUE), "CsparseMatrix")
-g1W <- graph.adjacency(W, mode="directed", weighted="W")
-c1W <- clusters(g1W)
+g1W <- graph_from_adjacency_matrix(W, mode="directed", weighted="W")
+c1W <- components(g1W)
 all.equal(c1W$membership, res$comp.id)
    
 
    ## [1] "names for target but not for current"
    
 
    
@@ -640,27 +622,16 @@ 
    Shortest paths in weights mat
 above. The diameter measure is then the diameter of the largest
 component subgraph. Note that this generates an n x
 n matrix:
    
-
-
    ## Warning: `is.connected()` was deprecated in igraph 2.0.0.
-## ℹ Please use `is_connected()` instead.
-## This warning is displayed once every 8 hours.
-## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
-## generated.
    
+
 
    ## [1] FALSE
    
-
    +
     dg1 <- diameter(g1)
 dg1
    
 
    ## [1] 7
    
-
    -sp_mat <- shortest.paths(g1)
    
-
    ## Warning: `shortest.paths()` was deprecated in igraph 2.0.0.
-## ℹ Please use `distances()` instead.
-## This warning is displayed once every 8 hours.
-## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
-## generated.
    
-
    -str(sp_mat)
    
+
    +sp_mat <- distances(g1)
+str(sp_mat)
    
 
    ##  num [1:49, 1:49] 0 1 1 2 2 3 4 3 3 4 ...
 ##  - attr(*, "dimnames")=List of 2
 ##   ..$ : chr [1:49] "1" "2" "3" "4" ...
@@ -674,7 +645,7 @@ 
    Shortest paths in weights matr
 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)
@@ -686,7 +657,7 @@ 
    Shortest paths in weights matr
 produced by shortest.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:
    
-
    +
     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]]
@@ -720,23 +691,23 @@ 
    Smirnov/Anselin (2009) cyclical m
 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:
    
-
    +
 
    ## [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:
    
-
    +
     spdep:::find_q1_q2(nb_W)
    
 
    ## [1] 2 0
    
-
    +
 
    ## [1] -1.544645  1.000000
    
 
    By construction, all two-node connected graph components also meet
@@ -763,7 +734,7 @@ 
    Smirnov/Anselin (2009) cyclical m
 
 
    
   
    
-
    Site built with pkgdown 2.0.7.
    
+
    Site built with pkgdown 2.0.9.
    
 
    
 
       
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 68226ef..889845c 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/nb_igraph_files/figure-html/unnamed-chunk-5-1.png b/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-5-1.png
index ee1309f..26e4120 100644
Binary files a/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-5-1.png and b/docs/articles/nb_igraph_files/figure-html/unnamed-chunk-5-1.png differ
diff --git a/docs/articles/sids_models.html b/docs/articles/sids_models.html
index ce56078..d4b78d5 100644
--- a/docs/articles/sids_models.html
+++ b/docs/articles/sids_models.html
@@ -6,7 +6,7 @@
 
 
 North Carolina SIDS data set (models) • spatialreg
-
+
 
 
 
@@ -33,7 +33,7 @@
       
       
         spatialreg
-        1.3-2
+        1.3-3
       
     
    
 
@@ -50,14 +50,14 @@

  • @@ -121,7 +121,7 @@

    Getting the data into R

    We will be using the spdep and spatialreg packages, here version: spdep, version -1.3-2, 2024-01-17, the sf package and the +1.3-4, 2024-03-27, 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 @@ -279,11 +279,11 @@

    CAR model fitting## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) -## (Intercept) 1.4371520 0.2252729 6.3796 1.775e-10 +## (Intercept) 1.4371519 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.55147 +## Numerical Hessian standard error of lambda: 0.5537 ## ## Log likelihood: -114.0376 ## ML residual variance (sigma squared): 1201.5, (sigma: 34.663) @@ -358,7 +358,7 @@

    References

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    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 79e8f99..9d48d8d 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 94021a1..de4505d 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 8631e1f..c2ce1d3 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -1,5 +1,5 @@ -Authors and Citation • spatialregAuthors and Citation • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@

  • Changelog @@ -60,7 +60,7 @@
    -

    Authors

    +

    Authors and Citation

    @@ -242,7 +242,7 @@

    Citation

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/index.html b/docs/index.html index aa31d2b..e72f3ce 100644 --- a/docs/index.html +++ b/docs/index.html @@ -6,7 +6,7 @@ Spatial Regression Analysis • spatialreg - + @@ -33,7 +33,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -50,14 +50,14 @@
  • @@ -143,7 +143,7 @@

    Developers

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/news/index.html b/docs/news/index.html index 0896f4c..bb534c0 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -1,5 +1,5 @@ -Changelog • spatialregChangelog • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -64,7 +64,12 @@

    Changelog

    -

    Version 1.3-2 (development)

    +

    Version 1.3-2 (development)2024-02-06

    +

    • add sqrt=TRUE in calls to Matrix determinant methods in matrix_ldet

      • +

    • add igraph (>= 2.0.0) in DESCRIPTION for re-named igraph functions

      • +
    +
    +

    Version 1.3-2 (2024-02-06)2024-02-06

    • pass through SlX formula in call

    • re-corrected #19 because the fitted model weights component may be NULL

    • suppress warning from multcomp::glht as the test which throws the warning is discarded

      • @@ -109,18 +114,18 @@

      Versi

    Version 1.1-8 (2021-05-03)2021-05-03

    -

    • #18 standardize use of coef() methods for (some) fitted model summary objects

      • +

      • #18 standardize use of coef() methods for (some) fitted model summary objects

      • https://github.com/tidymodels/broom/issues/1003#issuecomment-798694400 changing spatialreg model output class names: spdep sarlm -> spatialreg Sarlm, spautolm -> Spautolm, stsls -> Stsls, gmsar -> Gmsar, lagmess -> Lagmess, SLX -> , SlX, MCMC_s*_g -> MCMC_s*_G, SFResult -> SfResult, ME_res -> Me_res, lagImpact -> LagImpact, WXImpact -> WXimpact

        • -

      • #16 merged coordination of impacts methods (Gianfranco Piras)

        • -

      • #14 merged correction to SDEM and SLX impacts when a lagged intercept is present (Tobias Rüttenauer).

        • +

      • #16 merged coordination of impacts methods (Gianfranco Piras)

        • +

      • #14 merged correction to SDEM and SLX impacts when a lagged intercept is present (Tobias Rüttenauer).

    Version 1.1-5 (2019-12-01)2019-12-01

    -

    • #6, #11 na.action and precomputed eigenvalue bug

      • -

    • #9 Griddy Gibbs issue

      • -

    • #8 Predict method for SLX

      • -

    • #7, #13-14 Offset impacts for SDEM/SLX

      • -

    • #5, #10 Panel Durbin= argument

      • +

      • #6, #11 na.action and precomputed eigenvalue bug

        • +

      • #9 Griddy Gibbs issue

        • +

      • #8 Predict method for SLX

        • +

      • #7, #13-14 Offset impacts for SDEM/SLX

        • +

      • #5, #10 Panel Durbin= argument

    Version 1.1-3 (2019-04-01)2019-04-01

    @@ -140,7 +145,7 @@

    Versi

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 25030c2..4cda807 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,9 +1,9 @@ pandoc: 3.1.3 -pkgdown: 2.0.7 +pkgdown: 2.0.9 pkgdown_sha: ~ articles: - SpatialFiltering: SpatialFiltering.html nb_igraph: nb_igraph.html sids_models: sids_models.html -last_built: 2024-02-06T11:49Z + SpatialFiltering: SpatialFiltering.html +last_built: 2024-04-29T13:14Z diff --git a/docs/reference/GMerrorsar-1.png b/docs/reference/GMerrorsar-1.png index b35453b..1e5b37a 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 447012a..cf5d8a7 100644 --- a/docs/reference/GMerrorsar.html +++ b/docs/reference/GMerrorsar.html @@ -1,5 +1,5 @@ -Spatial simultaneous autoregressive error model estimation by GMM — GMerrorsar • spatialregSpatial simultaneous autoregressive error model estimation by GMM — GMerrorsar • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -370,9 +370,7 @@

    Examples

    "misc/nyadjwts.dbf", package="spData")[1]))[-1])) 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") +listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") summary(esar1f) @@ -465,7 +463,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/MCMCsamp.html b/docs/reference/MCMCsamp.html index b34c604..901b0e9 100644 --- a/docs/reference/MCMCsamp.html +++ b/docs/reference/MCMCsamp.html @@ -1,5 +1,5 @@ -MCMC sample from fitted spatial regression — MCMCsamp • spatialregMCMC sample from fitted spatial regression — MCMCsamp • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -137,9 +137,7 @@

    Examples

    identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10)) #> [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") +listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") summary(esar1f) @@ -195,28 +193,228 @@

    Examples

    #> PCTAGE65P 2.672002 3.33982 3.67753 4.04290 4.74392 #> PCTOWNHOME -0.845780 -0.56384 -0.41154 -0.21860 0.07637 #> -if (FALSE) { +# \dontrun{ 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.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: NA (not available for weighted model) +#> 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.01296 0.01568 0.0002218 0.0008193 +#> (Intercept) -0.79417 0.15064 0.0021303 0.0085373 +#> PEXPOSURE 0.07886 0.02974 0.0004206 0.0016661 +#> PCTAGE65P 3.79201 0.58160 0.0082250 0.0316924 +#> PCTOWNHOME -0.38114 0.16875 0.0023864 0.0100745 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> lambda -0.01785 0.00281 0.01326 0.02352 0.04254 +#> (Intercept) -1.09757 -0.89163 -0.79274 -0.69957 -0.47937 +#> PEXPOSURE 0.02125 0.05871 0.07958 0.09926 0.13413 +#> PCTAGE65P 2.68865 3.36618 3.80173 4.17535 4.89445 +#> PCTOWNHOME -0.70327 -0.49343 -0.37753 -0.27046 -0.03549 +#> 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.030872 +#> +#> 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.08485 0.03007 0.0004252 0.001841 +#> (Intercept) -0.66321 0.21541 0.0030463 0.014591 +#> PEXPOSURE 0.08242 0.04990 0.0007057 0.003038 +#> PCTAGE65P 3.65269 0.64308 0.0090945 0.039341 +#> PCTOWNHOME -0.35759 0.22515 0.0031840 0.014851 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> lambda 0.021075 0.06565 0.08492 0.1070 0.1398 +#> (Intercept) -1.173329 -0.78618 -0.63961 -0.5233 -0.2863 +#> PEXPOSURE -0.008489 0.04927 0.07836 0.1125 0.1967 +#> PCTAGE65P 2.400330 3.17755 3.65477 4.0908 4.9428 +#> PCTOWNHOME -0.761115 -0.51656 -0.37324 -0.2168 0.1507 +#> 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.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: NA (not available for weighted model) +#> 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.01421 0.01648 0.000233 0.0009767 +#> (Intercept) -0.79603 0.15920 0.002251 0.0097286 +#> PEXPOSURE 0.08092 0.03097 0.000438 0.0018774 +#> PCTAGE65P 3.77070 0.60576 0.008567 0.0381835 +#> PCTOWNHOME -0.37884 0.17038 0.002410 0.0102914 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> lambda -0.01732 0.002764 0.01370 0.02629 0.04646 +#> (Intercept) -1.12961 -0.899430 -0.80307 -0.69685 -0.48183 +#> PEXPOSURE 0.02262 0.059634 0.08016 0.10062 0.14212 +#> PCTAGE65P 2.58384 3.361796 3.81098 4.18747 4.90361 +#> PCTOWNHOME -0.69813 -0.495245 -0.38500 -0.26182 -0.03602 +#> 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.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: NA (not available for weighted model) +#> 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.03646 0.04020 0.0005686 0.002593 +#> (Intercept) -0.83219 0.15755 0.0022281 0.009562 +#> PEXPOSURE 0.09116 0.03438 0.0004863 0.002139 +#> PCTAGE65P 3.74091 0.59000 0.0083439 0.035836 +#> PCTOWNHOME -0.34906 0.17055 0.0024119 0.010630 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> lambda -0.04545 0.009909 0.04055 0.06473 0.108790 +#> (Intercept) -1.15043 -0.936965 -0.82876 -0.72797 -0.532623 +#> PEXPOSURE 0.02618 0.068120 0.08997 0.11245 0.165457 +#> PCTAGE65P 2.56456 3.350294 3.72433 4.13057 4.958447 +#> PCTOWNHOME -0.69327 -0.461490 -0.34685 -0.24507 -0.001185 +#> +# } esar0 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY) summary(esar0) @@ -259,53 +457,417 @@

    Examples

    #> plus standard error of the mean: #> #> Mean SD Naive SE Time-series SE -#> lambda 0.04086 0.01436 0.0004541 0.001639 -#> (Intercept) -0.62555 0.18275 0.0057791 0.024443 -#> PEXPOSURE 0.06900 0.04034 0.0012758 0.004965 -#> PCTAGE65P 3.68791 0.64881 0.0205172 0.077680 -#> PCTOWNHOME -0.38263 0.19769 0.0062513 0.025188 +#> lambda 0.04682 0.01895 0.0005993 0.002739 +#> (Intercept) -0.65775 0.21896 0.0069243 0.031633 +#> PEXPOSURE 0.08543 0.05260 0.0016634 0.007626 +#> PCTAGE65P 3.65044 0.57371 0.0181422 0.060983 +#> PCTOWNHOME -0.36116 0.23745 0.0075089 0.036472 #> #> 2. Quantiles for each variable: #> -#> 2.5% 25% 50% 75% 97.5% -#> lambda 0.014717 0.03100 0.04103 0.05272 0.065715 -#> (Intercept) -0.972641 -0.75486 -0.63637 -0.50400 -0.198074 -#> PEXPOSURE -0.004903 0.03891 0.06820 0.09560 0.152374 -#> PCTAGE65P 2.382993 3.31371 3.67406 4.12972 4.919868 -#> PCTOWNHOME -0.779951 -0.49064 -0.38080 -0.25234 -0.003732 +#> 2.5% 25% 50% 75% 97.5% +#> lambda 0.0109828 0.03340 0.04839 0.05958 0.08139 +#> (Intercept) -1.0771095 -0.81427 -0.62781 -0.51408 -0.23029 +#> PEXPOSURE -0.0004135 0.04738 0.07927 0.11564 0.20607 +#> PCTAGE65P 2.5911599 3.20769 3.63005 4.09097 4.74100 +#> PCTOWNHOME -0.8037868 -0.54136 -0.35322 -0.20277 0.17058 #> -if (FALSE) { +# \dontrun{ 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.01177 0.01571 0.0002221 0.0008781 +#> (Intercept) -0.79575 0.14274 0.0020186 0.0078442 +#> PEXPOSURE 0.08036 0.02929 0.0004143 0.0017545 +#> PCTAGE65P 3.81500 0.57129 0.0080793 0.0330363 +#> PCTOWNHOME -0.38426 0.15508 0.0021932 0.0087813 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> lambda -0.01899 0.00175 0.01220 0.02213 0.04328 +#> (Intercept) -1.07330 -0.89063 -0.79366 -0.70020 -0.51399 +#> PEXPOSURE 0.02200 0.06082 0.08123 0.10067 0.13675 +#> PCTAGE65P 2.70738 3.41885 3.80976 4.20150 4.94589 +#> PCTOWNHOME -0.70906 -0.47931 -0.38839 -0.28703 -0.06345 +#> 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.02791 0.01696 0.0002398 0.001257 +#> (Intercept) -1.11089 0.24295 0.0034359 0.019066 +#> PEXPOSURE 0.20998 0.08066 0.0011407 0.006616 +#> PCTAGE65P 3.44550 0.61861 0.0087484 0.045120 +#> PCTOWNHOME 0.03249 0.25058 0.0035437 0.022380 +#> lag.(Intercept) 0.11788 0.05850 0.0008273 0.004548 +#> lag.PEXPOSURE -0.03381 0.01626 0.0002299 0.001358 +#> lag.PCTAGE65P 0.25893 0.23503 0.0033238 0.019521 +#> lag.PCTOWNHOME -0.15337 0.05824 0.0008236 0.004496 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> lambda -0.004792 0.01604 0.02829 0.03980 0.0616941 +#> (Intercept) -1.549205 -1.28823 -1.12499 -0.91520 -0.6509996 +#> PEXPOSURE 0.063165 0.15363 0.20577 0.26455 0.3711420 +#> PCTAGE65P 2.335975 3.01533 3.44864 3.91005 4.6644758 +#> PCTOWNHOME -0.505045 -0.11695 0.05420 0.19038 0.5171339 +#> lag.(Intercept) 0.006027 0.07616 0.11689 0.16220 0.2310606 +#> lag.PEXPOSURE -0.064845 -0.04488 -0.03385 -0.02339 -0.0009086 +#> lag.PCTAGE65P -0.195440 0.09430 0.24778 0.43584 0.6912104 +#> lag.PCTOWNHOME -0.263176 -0.19213 -0.15386 -0.11592 -0.0317401 +#> 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.03924 0.01531 0.0002166 0.0009369 +#> (Intercept) -0.51482 0.16721 0.0023647 0.0103798 +#> PEXPOSURE 0.05057 0.03374 0.0004771 0.0019411 +#> PCTAGE65P 3.58543 0.66827 0.0094508 0.0442173 +#> PCTOWNHOME -0.41083 0.18638 0.0026358 0.0118147 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> rho 0.01129 0.02805 0.03879 0.05012 0.07010 +#> (Intercept) -0.84632 -0.62884 -0.50852 -0.39309 -0.20739 +#> PEXPOSURE -0.01768 0.02886 0.04981 0.07205 0.11808 +#> PCTAGE65P 2.28982 3.12634 3.60888 4.03227 4.88549 +#> PCTOWNHOME -0.74664 -0.54787 -0.40497 -0.28735 -0.03815 +#> 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.02471 0.01629 0.0002304 0.001219 +#> (Intercept) -1.12989 0.24291 0.0034353 0.017687 +#> PEXPOSURE 0.21704 0.08665 0.0012254 0.007356 +#> PCTAGE65P 3.40551 0.61299 0.0086689 0.043441 +#> PCTOWNHOME 0.04146 0.25683 0.0036321 0.019859 +#> lag.(Intercept) 0.12868 0.05650 0.0007991 0.004535 +#> lag.PEXPOSURE -0.03472 0.01732 0.0002449 0.001526 +#> lag.PCTAGE65P 0.17214 0.20763 0.0029364 0.015671 +#> lag.PCTOWNHOME -0.13618 0.06186 0.0008749 0.005326 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> rho -0.008938 0.01339 0.02554 0.03617 0.054650 +#> (Intercept) -1.604105 -1.28894 -1.12504 -0.97119 -0.660270 +#> PEXPOSURE 0.053725 0.15881 0.21725 0.27236 0.389018 +#> PCTAGE65P 2.213387 2.99625 3.42281 3.85452 4.551724 +#> PCTOWNHOME -0.451879 -0.14346 0.04721 0.22476 0.545593 +#> lag.(Intercept) 0.019705 0.09016 0.13030 0.16767 0.243720 +#> lag.PEXPOSURE -0.067460 -0.04719 -0.03527 -0.02388 0.003145 +#> lag.PCTAGE65P -0.238603 0.04557 0.17728 0.29777 0.562834 +#> lag.PCTOWNHOME -0.246111 -0.18197 -0.13460 -0.09357 -0.018730 +#> 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.04897 0.07268 0.001028 0.02171 +#> lambda 0.07158 0.07109 0.001005 0.01911 +#> (Intercept) -0.79253 0.30171 0.004267 0.05929 +#> PEXPOSURE 0.12518 0.08005 0.001132 0.01428 +#> PCTAGE65P 3.33373 0.62239 0.008802 0.04114 +#> PCTOWNHOME -0.24558 0.24700 0.003493 0.03231 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> rho -0.148205 -0.10220 -0.07484 -0.0003653 0.09879 +#> lambda -0.097739 0.03709 0.10543 0.1218063 0.13809 +#> (Intercept) -1.371833 -1.01584 -0.80342 -0.5327538 -0.26199 +#> PEXPOSURE -0.001867 0.05674 0.12102 0.1859304 0.27450 +#> PCTAGE65P 2.139989 2.91933 3.33752 3.7143648 4.58803 +#> PCTOWNHOME -0.743597 -0.40260 -0.25567 -0.0792356 0.24964 +#> 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.005935 0.06247 0.0008835 0.012391 +#> lambda 0.025149 0.06578 0.0009302 0.012564 +#> (Intercept) -1.104966 0.25946 0.0036694 0.020691 +#> PEXPOSURE 0.220056 0.08072 0.0011415 0.006564 +#> PCTAGE65P 3.448064 0.68477 0.0096841 0.053253 +#> PCTOWNHOME 0.002352 0.25423 0.0035954 0.018914 +#> lag.(Intercept) 0.107623 0.06651 0.0009406 0.007925 +#> lag.PEXPOSURE -0.034004 0.01609 0.0002276 0.001352 +#> lag.PCTAGE65P 0.279208 0.31780 0.0044943 0.037387 +#> lag.PCTOWNHOME -0.134012 0.06451 0.0009123 0.005533 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> rho -0.10755 -0.05711 -0.014110 0.05341 0.099104 +#> lambda -0.09950 -0.03206 0.037324 0.08249 0.118177 +#> (Intercept) -1.60962 -1.28203 -1.096564 -0.93087 -0.597130 +#> PEXPOSURE 0.05429 0.16818 0.217623 0.27874 0.368060 +#> PCTAGE65P 2.11661 3.00219 3.455553 3.93036 4.829667 +#> PCTOWNHOME -0.46367 -0.19124 -0.002169 0.17287 0.544393 +#> lag.(Intercept) -0.02153 0.05985 0.109887 0.15841 0.231289 +#> lag.PEXPOSURE -0.06423 -0.04610 -0.033598 -0.02331 -0.001144 +#> lag.PCTAGE65P -0.29223 0.02822 0.266936 0.50401 0.926973 +#> lag.PCTOWNHOME -0.26326 -0.17692 -0.130161 -0.08812 -0.009731 +#> +# }     @@ -320,7 +882,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/ME.html b/docs/reference/ME.html index 99540ec..48ae762 100644 --- a/docs/reference/ME.html +++ b/docs/reference/ME.html @@ -1,5 +1,5 @@ -Moran eigenvector GLM filtering — ME • spatialregMoran eigenvector GLM filtering — ME • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -168,7 +168,7 @@

    Examples

    #> 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.300 0.006 1.315 +#> 1.202 0.007 1.222 glmME <- glm(c(hopkins_part) ~ 1 + fitted(MEbinom1), family="binomial") #anova(glmME, test="Chisq") coef(summary(glmME)) @@ -186,7 +186,7 @@

    Examples

    #> 2 253 275.39 2 16.841 0.0002203 *** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 -if (FALSE) { +# \dontrun{ require("sf", quietly=TRUE) columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE) #require("spdep", quietly=TRUE) @@ -195,30 +195,109 @@

    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.589 0.003 0.606 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) +#> Warning: subsetting caused increase in subgraph count +#> 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], @@ -230,12 +309,43 @@

    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 +# }     @@ -250,7 +360,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/ML_models.html b/docs/reference/ML_models.html index 61c3403..4cca555 100644 --- a/docs/reference/ML_models.html +++ b/docs/reference/ML_models.html @@ -1,5 +1,5 @@ -Spatial simultaneous autoregressive model estimation by maximum likelihood — ML_models • spatialregSpatial simultaneous autoregressive model estimation by maximum likelihood — ML_models • spatialreg spatialreg - 1.3-2 + 1.3-3 @@ -44,14 +44,14 @@
  • Changelog @@ -428,39 +428,147 @@

    Examples

    #> (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 -if (FALSE) { +# \dontrun{ COL.lag.eig$fdHess +#> [1] FALSE COL.lag.eig$resvar +#> sigma rho (Intercept) INC HOVAL +#> sigma 379.77510023 -0.3236306420 16.3015085 -0.29590802 -0.0202478469 +#> rho -0.32363064 0.0138487528 -0.6975717 0.01266245 0.0008664428 +#> (Intercept) 16.30150847 -0.6975716519 51.5143024 -1.32602702 -0.1616379845 +#> INC -0.29590802 0.0126624508 -1.3260270 0.09311223 -0.0117959714 +#> 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.079249 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.59742708 -0.701154300 35.3176470 -0.64109252 -0.043867493 +#> rho -0.70115430 0.030003687 -1.5113073 0.02743353 0.001877171 +#> (Intercept) 35.31764699 -1.511307337 92.5027548 -2.07005706 -0.212549096 +#> INC -0.64109252 0.027433533 -2.0700571 0.10661800 -0.010871823 +#> 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.87142 5.727 1.022e-08 +#> INC -1.03162 0.32843 -3.141 0.001683 +#> HOVAL -0.26593 0.08823 -3.014 0.002578 +#> +#> Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588 +#> Approximate (numerical Hessian) standard error: 0.12363 +#> z-value: 3.4865, p-value: 0.00048934 +#> Wald statistic: 12.156, p-value: 0.00048934 +#> +#> 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.0152832589 -0.8346578 0.02005909 0.0002832041 +#> (Intercept) -0.8346577895 61.9591934 -1.78361966 -0.1334688151 +#> INC 0.0200590872 -1.7836197 0.10786711 -0.0122201879 +#> HOVAL 0.0002832041 -0.1334688 -0.01222019 0.0077846116 # 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.01528326 # 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 @@ -476,7 +584,7 @@

    Examples

    #> Computing eigenvalues ... #> #> user system elapsed -#> 0.180 0.000 0.182 +#> 0.143 0.000 0.145 summary(COL.lag.M) #> #> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = listw, @@ -511,66 +619,455 @@

    Examples

    #> Direct Indirect Total #> INC -1.0860220 -0.7270848 -1.8131068 #> HOVAL -0.2799509 -0.1874254 -0.4673763 -if (FALSE) { +# \dontrun{ 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.394 0.000 0.396 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.1531e-05 1.2216e-01 0.0002 0.9998594 +#> 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.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 +#> +#> 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.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" +# } -if (FALSE) { +# \dontrun{ 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.89816225 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.38640403 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) + @@ -578,7 +1075,61 @@

    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.89816225 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.38640403 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)) #> @@ -740,62 +1291,216 @@

    Examples

    #> Number of parameters estimated: 7 #> AIC: 377.17, (AIC for lm: 380.16) #> -if (FALSE) { +# \dontrun{ 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.985784 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, glht, n): +#> Direct Indirect Total +#> DISCBD -5.4129365 NA -5.4129365 +#> INC -0.8994251 0.1618581 -0.7375670 +#> HOVAL -0.2028457 NA -0.2028457 +#> ======================================================== +#> Standard errors: +#> Direct Indirect Total +#> DISCBD 2.00928140 NA 2.00928140 +#> INC 0.31517363 0.6038588 0.67559401 +#> HOVAL 0.09012493 NA 0.09012493 +#> ======================================================== +#> Z-values: +#> Direct Indirect Total +#> DISCBD -2.693966 NA -2.693966 +#> INC -2.853745 0.2680396 -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.44133843 0.31874928 -34.47163074 2.42244758 -4.32095072 +#> 1007 1008 1009 1010 1011 1012 +#> 8.66744165 -13.38669934 -1.92276585 17.85753950 -1.11484596 -2.30434792 +#> 1013 1014 1015 1016 1017 1018 +#> -8.16935116 -5.80500231 0.14973721 5.93191445 -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.52940712 -9.60025384 -6.95635586 -0.43630579 5.98493664 +#> 1031 1032 1033 1034 1035 1036 +#> 6.25882669 7.75527032 10.83413236 23.23927250 -0.05594957 1.43808573 +#> 1037 1038 1039 1040 1041 1042 +#> 9.51995266 12.18295373 8.31031724 17.06834503 7.04418393 7.50088924 +#> 1043 1044 1045 1046 1047 1048 +#> -7.78485976 -6.79207447 -7.94977534 -11.25362117 -5.68994630 5.04837399 +#> 1049 +#> 2.22497381 print(system.time(ev <- eigenw(similar.listw(listw)))) +#> user system elapsed +#> 0.001 0.000 0.001 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.143 0.000 0.143 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.139 0.000 0.140 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.151 0.000 0.153 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)))) +#> 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 +#> user system elapsed +#> 0.159 0.000 0.160 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.158 0.000 0.160 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.158 0.000 0.159 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.144 0.000 0.144 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.144 0.000 0.144 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.143 0.000 0.145 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.150 0.001 0.152 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.152 0.000 0.152 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.149 0.000 0.149 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.168 0.000 0.169 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) @@ -810,7 +1515,7 @@

    Examples

    #> Type: sac #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 47.783764 9.902659 4.8253 1.398e-06 +#> (Intercept) 47.783766 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 #> @@ -833,8 +1538,8 @@

    Examples

    summary(impacts(COL.sacW.eig, tr=trMatc, R=2000), zstats=TRUE, short=TRUE) #> Impact measures (sac, trace): #> Direct Indirect Total -#> INC -1.0632723 -0.5601502 -1.6234225 -#> HOVAL -0.2919129 -0.1537848 -0.4456977 +#> INC -1.0632723 -0.5601501 -1.6234223 +#> HOVAL -0.2919129 -0.1537847 -0.4456977 #> ======================================================== #> Simulation results ( variance matrix): #> ======================================================== @@ -845,13 +1550,13 @@

    Examples

    #> #> Simulated z-values: #> Direct Indirect Total -#> INC -3.376570 -0.8371025 -1.861677 -#> HOVAL -3.159909 -0.7431010 -1.553156 +#> INC -3.376570 -0.8371024 -1.861677 +#> HOVAL -3.159909 -0.7431009 -1.553156 #> #> Simulated p-values: #> Direct Indirect Total -#> INC 0.00073396 0.40253 0.062649 -#> HOVAL 0.00157819 0.45742 0.120386 +#> INC 0.00073396 0.40254 0.062649 +#> HOVAL 0.00157818 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) @@ -866,14 +1571,14 @@

    Examples

    #> Type: sacmixed #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 50.92026 68.25721 0.7460 0.455664 +#> (Intercept) 50.92026 68.25722 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 #> lag.HOVAL 0.20852 0.28702 0.7265 0.467546 #> #> Rho: 0.31557 -#> Asymptotic standard error: 0.9458 +#> Asymptotic standard error: 0.94581 #> z-value: 0.33365, p-value: 0.73864 #> Lambda: 0.15415 #> Asymptotic standard error: 1.0643 @@ -897,19 +1602,19 @@

    Examples

    #> Simulation results ( variance matrix): #> ======================================================== #> Simulated standard errors -#> Direct Indirect Total -#> INC 0.3778064 2.1709289 2.3099448 -#> HOVAL 0.1144714 0.7736216 0.8389226 +#> Direct Indirect Total +#> INC 0.3799971 2.356727 2.5113549 +#> HOVAL 0.1091241 0.761200 0.8220231 #> #> Simulated z-values: #> Direct Indirect Total -#> INC -2.737182 -0.7297008 -1.13347010 -#> HOVAL -2.364129 0.2526689 -0.08958505 +#> INC -2.706969 -0.6434228 -1.01340192 +#> HOVAL -2.478740 0.2512075 -0.09643407 #> #> Simulated p-values: #> Direct Indirect Total -#> INC 0.0061968 0.46557 0.25702 -#> HOVAL 0.0180725 0.80052 0.92862 +#> INC 0.0067901 0.51995 0.31087 +#> HOVAL 0.0131847 0.80165 0.92318 COL.msacW1.eig <- sacsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, listw, Durbin=TRUE, control=list(pre_eig1=ev, pre_eig2=ev)) summary(COL.msacW1.eig) @@ -924,14 +1629,14 @@

    Examples

    #> Type: sacmixed #> Coefficients: (asymptotic standard errors) #> Estimate Std. Error z value Pr(>|z|) -#> (Intercept) 50.92026 68.25721 0.7460 0.455664 +#> (Intercept) 50.92026 68.25722 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 #> lag.HOVAL 0.20852 0.28702 0.7265 0.467546 #> #> Rho: 0.31557 -#> Asymptotic standard error: 0.9458 +#> Asymptotic standard error: 0.94581 #> z-value: 0.33365, p-value: 0.73864 #> Lambda: 0.15415 #> Asymptotic standard error: 1.0643 @@ -955,19 +1660,19 @@

    Examples

    #> Simulation results ( variance matrix): #> ======================================================== #> Simulated standard errors -#> Direct Indirect Total -#> INC 0.3778064 2.1709289 2.3099448 -#> HOVAL 0.1144714 0.7736216 0.8389226 +#> Direct Indirect Total +#> INC 0.3799971 2.356727 2.5113549 +#> HOVAL 0.1091241 0.761200 0.8220231 #> #> Simulated z-values: #> Direct Indirect Total -#> INC -2.737182 -0.7297008 -1.13347010 -#> HOVAL -2.364129 0.2526689 -0.08958505 +#> INC -2.706969 -0.6434228 -1.01340192 +#> HOVAL -2.478740 0.2512075 -0.09643407 #> #> Simulated p-values: #> Direct Indirect Total -#> INC 0.0061968 0.46557 0.25702 -#> HOVAL 0.0180725 0.80052 0.92862 +#> INC 0.0067901 0.51995 0.31087 +#> HOVAL 0.0131847 0.80165 0.92318 COL.msacW2.eig <- sacsarlm(CRIME ~ DISCBD + INC + HOVAL, data=COL.OLD, listw, Durbin= ~ INC, control=list(pre_eig1=ev, pre_eig2=ev)) summary(COL.msacW2.eig) @@ -1015,7 +1720,7 @@

    Examples

    #> ======================================================== #> Simulated standard errors #> Direct Indirect Total -#> DISCBD 3.1446963 5.7201656 5.757459 +#> DISCBD 3.1446963 5.7201657 5.757459 #> INC 0.3621137 1.7150155 1.876353 #> HOVAL 0.1090172 0.6360488 0.698836 #> @@ -1030,17 +1735,133 @@

    Examples

    #> DISCBD 0.064026 0.97766 0.32525 #> INC 0.013822 0.85235 0.76037 #> HOVAL 0.043203 0.85389 0.62910 -if (FALSE) { +# \dontrun{ 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.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.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.822416 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.3572354 +#> 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 +#> +# }     @@ -1055,7 +1876,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/Rplot001.png b/docs/reference/Rplot001.png index 881792b..f4c4075 100644 Binary files a/docs/reference/Rplot001.png and b/docs/reference/Rplot001.png differ diff --git a/docs/reference/Rplot002.png b/docs/reference/Rplot002.png index 999a7c9..e54cd94 100644 Binary files a/docs/reference/Rplot002.png and b/docs/reference/Rplot002.png differ diff --git a/docs/reference/Rplot003.png b/docs/reference/Rplot003.png index ae0ec1c..26d3328 100644 Binary files a/docs/reference/Rplot003.png and b/docs/reference/Rplot003.png differ diff --git a/docs/reference/SET_MCMC.html b/docs/reference/SET_MCMC.html index 4f4e74a..fae47cb 100644 --- a/docs/reference/SET_MCMC.html +++ b/docs/reference/SET_MCMC.html @@ -1,5 +1,5 @@ -Bayesian MCMC spatial simultaneous autoregressive model estimation — spBreg_lag • spatialregBayesian MCMC spatial simultaneous autoregressive model estimation — spBreg_lag • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -333,7 +333,7 @@

    Examples

    #> lambda 3 1010 937 1.080 #> sige 2 930 937 0.993 #> -if (FALSE) { +# \dontrun{ ev <- eigenw(lw) W <- as(lw, "CsparseMatrix") trMatc <- trW(W, type="mult") @@ -341,57 +341,616 @@

    Examples

    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.66146 11.7776 0.263354 0.263354 +#> INC -1.01262 0.3770 0.008430 0.008430 +#> HOVAL -0.28032 0.1057 0.002363 0.002363 +#> lag.INC -1.07754 0.7304 0.016333 0.016333 +#> lag.HOVAL 0.09151 0.2395 0.005356 0.005702 +#> lambda 0.48878 0.1732 0.003873 0.003873 +#> sige 114.61924 26.8321 0.599984 0.653970 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> (Intercept) 48.2063 65.55434 73.00434 80.0130 95.69254 +#> INC -1.7385 -1.27266 -1.01806 -0.7591 -0.25126 +#> HOVAL -0.4883 -0.35221 -0.27730 -0.2111 -0.07194 +#> lag.INC -2.4291 -1.54241 -1.09764 -0.6156 0.38091 +#> lag.HOVAL -0.3892 -0.06269 0.09028 0.2468 0.55625 +#> lambda 0.1165 0.37669 0.49775 0.6118 0.79690 +#> sige 73.4822 95.76803 110.90096 128.2633 177.00970 +#> print(summary(impacts(COL.err.Bayes))) +#> Impact measures (SDEM, MCMC, n): +#> Direct Indirect Total +#> INC -1.0126232 -1.07753574 -2.090159 +#> HOVAL -0.2803178 0.09150584 -0.188812 +#> ======================================================== +#> Standard errors: +#> Direct Indirect Total +#> INC 0.3572678 0.6921465 0.8161704 +#> HOVAL 0.1001564 0.2269602 0.2729228 +#> ======================================================== +#> Z-values: +#> Direct Indirect Total +#> INC -2.834354 -1.5568031 -2.5609345 +#> HOVAL -2.798799 0.4031801 -0.6918146 +#> +#> p-values: +#> Direct Indirect Total +#> INC 0.0045918 0.11952 0.010439 +#> HOVAL 0.0051293 0.68682 0.489054 +#> 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 1096 937 1.170 +#> HOVAL 3 1052 937 1.120 +#> lag.INC 2 892 937 0.952 +#> lag.HOVAL 2 930 937 0.993 +#> 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) 71.32778 18.0269 0.403094 0.382389 +#> INC -0.99169 0.3978 0.008896 0.008896 +#> HOVAL -0.27840 0.1115 0.002494 0.002494 +#> lag.INC -0.95360 0.8293 0.018543 0.021986 +#> lag.HOVAL 0.06812 0.2557 0.005717 0.005904 +#> lambda 0.58149 0.1895 0.004238 0.010821 +#> sige 119.53101 28.0026 0.626158 0.839390 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> (Intercept) 37.9514 62.67210 72.07446 80.4143 99.84194 +#> INC -1.7844 -1.24565 -0.99600 -0.7380 -0.21146 +#> HOVAL -0.4977 -0.35256 -0.27741 -0.2077 -0.05449 +#> lag.INC -2.4414 -1.50144 -0.99500 -0.4736 0.79318 +#> lag.HOVAL -0.4571 -0.09065 0.07663 0.2398 0.55644 +#> lambda 0.1781 0.45916 0.59241 0.7231 0.90210 +#> sige 75.6501 99.39723 116.16515 135.1114 184.61449 +#> print(summary(impacts(COL.err.Bayes))) +#> Impact measures (SDEM, MCMC, n): +#> Direct Indirect Total +#> INC -0.9916937 -0.95359651 -1.945290 +#> HOVAL -0.2784029 0.06812087 -0.210282 +#> ======================================================== +#> Standard errors: +#> Direct Indirect Total +#> INC 0.3770053 0.7858279 0.9655825 +#> HOVAL 0.1056706 0.2422841 0.2998021 +#> ======================================================== +#> Z-values: +#> Direct Indirect Total +#> INC -2.63045 -1.2134929 -2.0146288 +#> HOVAL -2.63463 0.2811611 -0.7014028 +#> +#> p-values: +#> Direct Indirect Total +#> INC 0.0085272 0.22494 0.043944 +#> HOVAL 0.0084229 0.77859 0.483052 +#> 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 930 937 0.993 +#> HOVAL 2 930 937 0.993 +#> lag.INC 3 1013 937 1.080 +#> lag.HOVAL 2 969 937 1.030 +#> lambda 17 4647 937 4.960 +#> sige 3 1052 937 1.120 +#> 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.5397 10.99197 0.245788 0.245788 +#> INC -1.0322 0.38342 0.008574 0.008574 +#> HOVAL -0.2861 0.09986 0.002233 0.002233 +#> lag.INC -0.9210 0.59074 0.013209 0.012745 +#> lambda 0.4879 0.17196 0.003845 0.003672 +#> sige 112.0817 25.47154 0.569561 0.568218 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> (Intercept) 52.0865 67.8131 74.9358 81.6887 95.01447 +#> INC -1.7848 -1.2893 -1.0382 -0.7802 -0.24806 +#> HOVAL -0.4794 -0.3497 -0.2860 -0.2203 -0.08916 +#> lag.INC -2.0803 -1.3131 -0.9206 -0.5438 0.27869 +#> lambda 0.1355 0.3737 0.4987 0.6088 0.80100 +#> sige 73.0279 94.1912 108.3141 125.5058 173.13320 +#> print(summary(impacts(COL.err.Bayes))) +#> Impact measures (SDEM, MCMC, n): +#> Direct Indirect Total +#> INC -1.0321950 -0.9210334 -1.9532284 +#> HOVAL -0.2860547 NA -0.2860547 +#> ======================================================== +#> Standard errors: +#> Direct Indirect Total +#> INC 0.36744116 0.5661167 0.69838314 +#> HOVAL 0.09570117 NA 0.09570117 +#> ======================================================== +#> Z-values: +#> Direct Indirect Total +#> INC -2.809144 -1.626932 -2.796786 +#> HOVAL -2.989041 NA -2.989041 +#> +#> p-values: +#> Direct Indirect Total +#> INC 0.0049673 0.10375 0.0051614 +#> HOVAL 0.0027985 NA 0.0027985 +#> 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 930 937 0.993 +#> HOVAL 2 930 937 0.993 +#> lag.INC 2 930 937 0.993 +#> 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.7704 14.13707 0.316114 0.351470 +#> INC -0.9865 0.39508 0.008834 0.008834 +#> HOVAL -0.2904 0.09878 0.002209 0.002209 +#> lag.INC -0.8506 0.67466 0.015086 0.016187 +#> lambda 0.5688 0.17577 0.003930 0.009250 +#> sige 116.1069 27.61223 0.617428 0.789579 +#> +#> 2. Quantiles for each variable: +#> +#> 2.5% 25% 50% 75% 97.5% +#> (Intercept) 40.6365 65.4333 73.7510 81.2815 97.03388 +#> INC -1.7216 -1.2471 -1.0006 -0.7335 -0.15124 +#> HOVAL -0.4820 -0.3577 -0.2903 -0.2257 -0.09267 +#> lag.INC -2.0923 -1.2942 -0.8863 -0.4615 0.58595 +#> lambda 0.2251 0.4389 0.5731 0.6982 0.89484 +#> sige 74.0162 96.6611 112.0106 130.7735 181.83168 +#> print(summary(impacts(COL.err.Bayes))) +#> Impact measures (SDEM, MCMC, n): +#> Direct Indirect Total +#> INC -0.9865093 -0.8506053 -1.8371146 +#> HOVAL -0.2904128 NA -0.2904128 +#> ======================================================== +#> Standard errors: +#> Direct Indirect Total +#> INC 0.37860903 0.6465345 0.82479429 +#> HOVAL 0.09466321 NA 0.09466321 +#> ======================================================== +#> Z-values: +#> Direct Indirect Total +#> INC -2.605615 -1.315638 -2.227361 +#> HOVAL -3.067853 NA -3.067853 +#> +#> p-values: +#> Direct Indirect Total +#> INC 0.009171 0.1883 0.025923 +#> HOVAL 0.002156 NA 0.002156 +#> 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) 10 3436 937 3.670 +#> INC 2 930 937 0.993 +#> HOVAL 2 892 937 0.952 +#> lag.INC 3 1143 937 1.220 +#> lambda 13 3579 937 3.820 +#> sige 3 1096 937 1.170 +#> 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.1834 24.14830 0.763636 3.169805 +#> INC -0.9888 0.36234 0.011458 0.016240 +#> HOVAL -0.2861 0.09785 0.003094 0.003094 +#> lag.INC -0.8734 0.77381 0.024470 0.060097 +#> lag.HOVAL 0.1735 0.21943 0.006939 0.013468 +#> 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.280965 56.1010 76.4887 113.5316 +#> INC -1.6963 -1.228821 -0.9813 -0.7271 -0.3299 +#> HOVAL -0.4747 -0.350819 -0.2877 -0.2214 -0.1097 +#> lag.INC -2.5274 -1.342333 -0.8558 -0.3281 0.5407 +#> lag.HOVAL -0.2876 0.037091 0.1824 0.3249 0.5910 +#> 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.0341008 -1.2391410 -2.2732417 +#> HOVAL -0.2808381 0.1434285 -0.1374096 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> INC 0.3510018 0.9466895 1.0370366 +#> HOVAL 0.1010491 0.2773738 0.3171738 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> INC -2.966104 -1.4250264 -2.3048033 +#> HOVAL -2.775326 0.5576145 -0.3965541 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> INC 0.0030160 0.15415 0.021178 +#> HOVAL 0.0055146 0.57711 0.691696 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.3105 13.74662 0.307384 0.307384 +#> INC -0.9246 0.36837 0.008237 0.008237 +#> HOVAL -0.2959 0.09778 0.002187 0.002126 +#> lag.INC -0.5917 0.63940 0.014297 0.014297 +#> lag.HOVAL 0.2443 0.18959 0.004239 0.004239 +#> 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.18483 35.9928 45.0377 54.5246 72.8590 +#> INC -1.62602 -1.1728 -0.9198 -0.6633 -0.2117 +#> HOVAL -0.48546 -0.3597 -0.2944 -0.2310 -0.1012 +#> lag.INC -1.83696 -1.0334 -0.5786 -0.1739 0.6326 +#> lag.HOVAL -0.13032 0.1164 0.2417 0.3715 0.6133 +#> 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.0277896 -1.471320 -2.49911007 +#> HOVAL -0.2820493 0.197023 -0.08502627 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> INC 0.3831777 1.3179822 1.4600855 +#> HOVAL 0.1017455 0.3302926 0.3686927 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> INC -2.735660 -1.2618796 -1.8570000 +#> HOVAL -2.785094 0.5880284 -0.2417988 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> INC 0.0062255 0.20699 0.063311 +#> HOVAL 0.0053512 0.55651 0.808936 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.4193 13.40915 0.299838 0.299838 +#> DISCBD -4.7794 2.14113 0.047877 0.047877 +#> INC -0.9422 0.34677 0.007754 0.007754 +#> HOVAL -0.1806 0.09837 0.002200 0.002200 +#> lag.INC 0.4213 0.63243 0.014141 0.014141 +#> 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.6341 47.656256 56.4238 65.0877 84.808823 +#> DISCBD -9.2565 -6.160525 -4.7536 -3.3735 -0.714089 +#> INC -1.6336 -1.174123 -0.9427 -0.7090 -0.263713 +#> HOVAL -0.3719 -0.247402 -0.1791 -0.1148 0.006329 +#> lag.INC -0.8379 -0.005461 0.4129 0.8399 1.648619 +#> rho -0.2116 0.072536 0.1836 0.3067 0.535793 +#> sige 67.8946 87.054945 100.9815 116.4528 161.258258 +#> print(summary(impacts(COL.D1.Bayes, tr=trMatc), short=TRUE, zstats=TRUE)) +#> Impact measures (mixed, trace): +#> Direct Indirect Total +#> DISCBD -4.8194767 -1.06169757 -5.8811743 +#> INC -0.9312352 0.29032298 -0.6409122 +#> HOVAL -0.1821150 -0.04011868 -0.2222336 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> DISCBD 2.1896626 1.82105239 3.3751922 +#> INC 0.3494541 0.81612759 0.8914013 +#> HOVAL 0.1001927 0.06787487 0.1430090 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> DISCBD -2.224202 -0.7454394 -1.8451502 +#> INC -2.683306 0.3286082 -0.7510717 +#> HOVAL -1.835800 -0.7422472 -1.6384544 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> DISCBD 0.0261348 0.45601 0.065016 +#> INC 0.0072898 0.74245 0.452610 +#> HOVAL 0.0663872 0.45794 0.101327 #data(elect80, package="spData") #lw <- spdep::nb2listw(e80_queen, zero.policy=TRUE) #el_ml <- lagsarlm(log(pc_turnout) ~ log(pc_college) + log(pc_homeownership) @@ -403,7 +962,7 @@

    Examples

    #print(summary(el_B)) #print(el_ml$timings) #print(attr(el_B, "timings")) -} +# }     @@ -418,7 +977,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/SLX.html b/docs/reference/SLX.html index d6653d7..67f8dce 100644 --- a/docs/reference/SLX.html +++ b/docs/reference/SLX.html @@ -1,5 +1,5 @@ -Spatial Durbin linear (SLX, spatially lagged X) model — lmSLX • spatialregSpatial Durbin linear (SLX, spatially lagged X) model — lmSLX • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -314,7 +314,7 @@

    Examples

    #> Direct Indirect Total #> INC 8.2671e-06 0.024029 5.1387e-11 #> -if (FALSE) { +# \dontrun{ crds <- cbind(COL.OLD$X, COL.OLD$Y) mdist <- sqrt(sum(diff(apply(crds, 2, range))^2)) dnb <- spdep::dnearneigh(crds, 0, mdist) @@ -328,24 +328,70 @@

    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) +#> +#> Coefficients: +#> Estimate Std. Error t value Pr(>|t|) +#> (Intercept) 18.78498 11.23795 1.67157 0.10187 +#> INC -0.65428 0.29463 -2.22064 0.03169 +#> HOVAL -0.18244 0.08138 -2.24174 0.03019 +#> lag..Intercept. 6.27790 4.44284 1.41304 0.16484 +#> lag.INC -0.12958 0.28794 -0.45002 0.65496 +#> lag.HOVAL 0.02668 0.11273 0.23664 0.81406 +#> summary(impacts(SLX)) -} +#> Impact measures (SlX, glht, 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) +#> Error in predict.SlX(COL.SLX, newdata = COL.OLD, listw = lw): mismatch between newdata and spatial weights. newdata should have region.id as row.names all.equal(pslx0, pslx1) -#> [1] TRUE +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'current' in selecting a method for function 'all.equal': object 'pslx1' not found COL.OLD1 <- COL.OLD COL.OLD1$INC <- COL.OLD1$INC + 1 pslx2 <- predict(COL.SLX, newdata=COL.OLD1, listw=lw) +#> Error in predict.SlX(COL.SLX, newdata = COL.OLD1, listw = lw): mismatch between newdata and spatial weights. newdata should have region.id as row.names sum(coef(COL.SLX)[c(2,4)]) -#> [1] -2.479902 +#> [1] 5.623621 mean(pslx2-pslx1) -#> [1] -2.479902 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'mean': object 'pslx2' not found @@ -360,7 +406,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/SpatialFiltering.html b/docs/reference/SpatialFiltering.html index 8848dff..c55b176 100644 --- a/docs/reference/SpatialFiltering.html +++ b/docs/reference/SpatialFiltering.html @@ -1,5 +1,5 @@ -Semi-parametric spatial filtering — SpatialFiltering • spatialregSemi-parametric spatial filtering — SpatialFiltering • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -389,7 +389,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/aple.html b/docs/reference/aple.html index a441c45..4798cea 100644 --- a/docs/reference/aple.html +++ b/docs/reference/aple.html @@ -1,5 +1,5 @@ -Approximate profile-likelihood estimator (APLE) — aple • spatialregApproximate profile-likelihood estimator (APLE) — aple • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -147,9 +147,10 @@

    Examples

    #> aple(as.vector(scale(wheat$yield_detrend, scale=FALSE)), spdep::nb2listw(nbr12, style="W")) #> [1] 0.6601805 -if (FALSE) { +# \dontrun{ errorsarlm(yield_detrend ~ 1, wheat, spdep::nb2listw(nbr12, style="W")) -} +#> Error in errorsarlm(yield_detrend ~ 1, wheat, spdep::nb2listw(nbr12, style = "W")): argument "Durbin" is missing, with no default +# } @@ -164,7 +165,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/aple.mc-1.png b/docs/reference/aple.mc-1.png index 5a3a386..44ce7f6 100644 Binary files a/docs/reference/aple.mc-1.png and b/docs/reference/aple.mc-1.png differ diff --git a/docs/reference/aple.mc.html b/docs/reference/aple.mc.html index 2b8014c..532dbb3 100644 --- a/docs/reference/aple.mc.html +++ b/docs/reference/aple.mc.html @@ -1,5 +1,5 @@ -Approximate profile-likelihood estimator (APLE) permutation test — aple.mc • spatialregApproximate profile-likelihood estimator (APLE) permutation test — aple.mc • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -111,7 +111,7 @@

    See also

    Examples

    - 
    if (FALSE) {
+    
    # \dontrun{
 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,7 +127,20 @@ 
    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 <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L
@@ -148,9 +161,21 @@ 
    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])
-}
+# }
    @@ -165,7 +190,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/aple.plot-1.png b/docs/reference/aple.plot-1.png index 59d556c..00777d6 100644 Binary files a/docs/reference/aple.plot-1.png and b/docs/reference/aple.plot-1.png differ diff --git a/docs/reference/aple.plot-2.png b/docs/reference/aple.plot-2.png index 9917d5f..279cf41 100644 Binary files a/docs/reference/aple.plot-2.png and b/docs/reference/aple.plot-2.png differ diff --git a/docs/reference/aple.plot-3.png b/docs/reference/aple.plot-3.png index 592d998..29b206b 100644 Binary files a/docs/reference/aple.plot-3.png and b/docs/reference/aple.plot-3.png differ diff --git a/docs/reference/aple.plot.html b/docs/reference/aple.plot.html index 1dd8e1e..77012c9 100644 --- a/docs/reference/aple.plot.html +++ b/docs/reference/aple.plot.html @@ -1,5 +1,5 @@ -Approximate profile-likelihood estimator (APLE) scatterplot — aple.plot • spatialregApproximate profile-likelihood estimator (APLE) scatterplot — aple.plot • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -127,7 +127,7 @@

    See also

    Examples

    - 
    if (FALSE) {
+    
    # \dontrun{
 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,18 +141,72 @@ 
    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)
-}
+# }
    @@ -167,7 +221,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/do_ldet.html b/docs/reference/do_ldet.html index c150a12..b10a757 100644 --- a/docs/reference/do_ldet.html +++ b/docs/reference/do_ldet.html @@ -1,5 +1,5 @@ -Spatial regression model Jacobian computations — do_ldet • spatialregSpatial regression model Jacobian computations — do_ldet • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -127,10 +127,10 @@

    Arguments

    fill-in initiation coefficient value, default 0.1

    Imult
    -

    see Cholesky; numeric scalar which defaults to zero. The matrix that is decomposed is A+m*I where m is the value of Imult and I is the identity matrix of order ncol(A). Default in calling spdep functions is 2, here it cannot be missing and does not have a default, but is rescaled for binary weights matrices in proportion to the maximim row sum in those calling functions

    +

    see Cholesky; numeric scalar which defaults to zero. The matrix that is decomposed is A+m*I where m is the value of Imult and I is the identity matrix of order ncol(A). Default in calling spdep functions is 2, here it cannot be missing and does not have a default, but is rescaled for binary weights matrices in proportion to the maximim row sum in those calling functions

    super
    -

    see Cholesky; logical scalar indicating is a supernodal decomposition should be created. The alternative is a simplicial decomposition. Default in calling spdep functions is FALSE for “Matrix_J” and as.logical(NA) for “Matrix”. Setting it to NA leaves the choice to a CHOLMOD-internal heuristic

    +

    see Cholesky; logical scalar indicating is a supernodal decomposition should be created. The alternative is a simplicial decomposition. Default in calling spdep functions is FALSE for “Matrix_J” and as.logical(NA) for “Matrix”. Setting it to NA leaves the choice to a CHOLMOD-internal heuristic

    order

    default FALSE; used in LU_prepermutate, note warnings given for lu method

    @@ -214,10 +214,10 @@

    Details

    negative sparse spatial weights matrix

    pChol
    -

    a “CHMfactor” from factorising csrw with Cholesky

    +

    a “CHMfactor” from factorising csrw with Cholesky

    nChol
    -

    a “CHMfactor” from factorising nW with Cholesky

    +

    a “CHMfactor” from factorising nW with Cholesky

    method

    string: “Matrix”

    @@ -491,7 +491,7 @@

    Author

    @@ -609,7 +609,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/eigenw.html b/docs/reference/eigenw.html index 45d29c9..6e91b66 100644 --- a/docs/reference/eigenw.html +++ b/docs/reference/eigenw.html @@ -1,5 +1,5 @@ -Spatial weights matrix eigenvalues — griffith_sone • spatialregSpatial weights matrix eigenvalues — griffith_sone • spatialreg spatialreg - 1.3-2 + 1.3-3 @@ -38,14 +38,14 @@
  • Changelog @@ -221,7 +221,7 @@

    Examples

    rg_GS <- griffith_sone(P=7, Q=7, type="rook") all.equal(rg_eig, rg_GS) #> [1] TRUE -if (FALSE) { +# \dontrun{ run <- FALSE if (require("RSpectra", quietly=TRUE)) run <- TRUE if (run) { @@ -230,18 +230,22 @@

    Examples

    resn <- eigs(B, k=1, which="SR")$values print(Re(c(resn, res1))) } +#> [1] -3.695518 3.695518 if (run) { print(all.equal(range(Re(rg_eig)), c(resn, res1))) } +#> [1] TRUE if (run) { lw <- spdep::nb2listw(rg, style="W") rg_eig <- eigenw(similar.listw(lw)) print(range(Re(rg_eig))) } +#> [1] -1 1 if (run) { W <- as(lw, "CsparseMatrix") print(Re(c(eigs(W, k=1, which="SR")$values, eigs(W, k=1, which="LR")$values))) -}} +}# } +#> [1] -1 1     @@ -256,7 +260,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/gstsls-1.png b/docs/reference/gstsls-1.png index 2d3b30b..9df8046 100644 Binary files a/docs/reference/gstsls-1.png and b/docs/reference/gstsls-1.png differ diff --git a/docs/reference/gstsls.html b/docs/reference/gstsls.html index 54147f6..370bfd2 100644 --- a/docs/reference/gstsls.html +++ b/docs/reference/gstsls.html @@ -1,5 +1,5 @@ -Spatial simultaneous autoregressive SAC model estimation by GMM — gstsls • spatialregSpatial simultaneous autoregressive SAC model estimation by GMM — gstsls • spatialregImpacts in spatial lag models — impacts • spatialregImpacts in spatial lag models — impacts • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -489,9 +489,13 @@

    Examples

    #> Direct Indirect Total #> INC -1.0968247 -0.8939687 -1.9907934 #> HOVAL -0.2781941 -0.1363596 -0.4145537 -if (FALSE) { +# \dontrun{ try(impacts(mobj, evalues=ev), silent=TRUE) -} +#> Impact measures (mixed, evalues): +#> Direct Indirect Total +#> INC -1.0418080 -1.4804246 -2.52223255 +#> HOVAL -0.2836325 0.2302055 -0.05342697 +# } summary(impacts(mobj, tr=trMatc, R=200), short=TRUE, zstats=TRUE) #> Impact measures (mixed, trace): #> Direct Indirect Total @@ -580,21 +584,135 @@

    Examples

    #> INC 0.00098873 0.03543 0.00026248 #> HOVAL 0.00226152 0.46502 0.56533971 #> -if (FALSE) { +# \dontrun{ mobj1 <- lagsarlm(CRIME ~ INC + HOVAL, columbus, listw, type="mixed", method="Matrix", control=list(fdHess=TRUE)) summary(mobj1) +#> +#> Call:lagsarlm(formula = CRIME ~ INC + HOVAL, data = columbus, listw = listw, +#> type = "mixed", method = "Matrix", control = list(fdHess = TRUE)) +#> +#> Residuals: +#> Min 1Q Median 3Q Max +#> -37.15904 -6.62594 -0.39823 6.57561 23.62757 +#> +#> Type: mixed +#> Coefficients: (asymptotic standard errors) +#> Estimate Std. Error z value Pr(>|z|) +#> (Intercept) 45.592896 14.240041 3.2017 0.0013660 +#> INC -0.939088 0.343069 -2.7373 0.0061942 +#> HOVAL -0.299605 0.090393 -3.3145 0.0009182 +#> lag.INC -0.618375 0.607708 -1.0176 0.3088903 +#> lag.HOVAL 0.266615 0.181524 1.4688 0.1418998 +#> +#> Rho: 0.38251, LR test value: 4.1648, p-value: 0.041272 +#> Asymptotic standard error: 0.17375 +#> z-value: 2.2015, p-value: 0.027702 +#> Wald statistic: 4.8465, p-value: 0.027702 +#> +#> Log likelihood: -182.0161 for mixed model +#> ML residual variance (sigma squared): 95.051, (sigma: 9.7494) +#> Number of observations: 49 +#> Number of parameters estimated: 7 +#> AIC: 378.03, (AIC for lm: 380.2) +#> LM test for residual autocorrelation +#> test value: 0.101, p-value: 0.75063 +#> set.seed(1) summary(impacts(mobj1, tr=trMatc, R=1000), zstats=TRUE, short=TRUE) +#> Impact measures (mixed, trace): +#> Direct Indirect Total +#> INC -1.0418080 -1.4804246 -2.52223255 +#> HOVAL -0.2836325 0.2302055 -0.05342697 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> INC 0.34686331 0.8324743 0.8959688 +#> HOVAL 0.09085583 0.3131242 0.3483785 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> INC -3.023902 -1.8284275 -2.8695190 +#> HOVAL -3.120537 0.7910739 -0.1028036 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> INC 0.0024954 0.067485 0.004111 +#> HOVAL 0.0018052 0.428901 0.918119 summary(impacts(mobj, tr=trMatc, R=1000), zstats=TRUE, short=TRUE) +#> Impact measures (mixed, trace): +#> Direct Indirect Total +#> INC -1.0418080 -1.4804246 -2.52223255 +#> HOVAL -0.2836325 0.2302055 -0.05342697 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> INC 0.3123750 0.8291512 0.8821167 +#> HOVAL 0.0922758 0.2961092 0.3267773 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> INC -3.330304 -1.8060414 -2.8769263 +#> HOVAL -3.099811 0.8140556 -0.1376722 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> INC 0.00086751 0.070912 0.0040157 +#> HOVAL 0.00193644 0.415613 0.8904995 mobj2 <- lagsarlm(CRIME ~ INC + HOVAL, columbus, listw, type="mixed", method="Matrix", control=list(fdHess=TRUE, optimHess=TRUE)) summary(impacts(mobj2, tr=trMatc, R=1000), zstats=TRUE, short=TRUE) +#> Impact measures (mixed, trace): +#> Direct Indirect Total +#> INC -1.0418080 -1.4804246 -2.52223255 +#> HOVAL -0.2836325 0.2302055 -0.05342697 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> INC 0.32853332 0.8044059 0.8615603 +#> HOVAL 0.09291975 0.3169361 0.3571304 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> INC -3.231204 -1.8690203 -2.9771674 +#> HOVAL -3.026494 0.7323897 -0.1374858 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> INC 0.0012327 0.06162 0.0029093 +#> HOVAL 0.0024741 0.46393 0.8906468 mobj3 <- lagsarlm(CRIME ~ INC + HOVAL, columbus, listw, type="mixed", method="spam", control=list(fdHess=TRUE)) summary(impacts(mobj3, tr=trMatc, R=1000), zstats=TRUE, short=TRUE) -} -if (FALSE) { +#> Impact measures (mixed, trace): +#> Direct Indirect Total +#> INC -1.0418080 -1.4804246 -2.52223255 +#> HOVAL -0.2836325 0.2302055 -0.05342697 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> INC 0.32866009 0.8664620 0.9215792 +#> HOVAL 0.09837849 0.3110899 0.3492674 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> INC -3.190433 -1.7325351 -2.766712 +#> HOVAL -2.890221 0.7694814 -0.128720 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> INC 0.0014206 0.083178 0.0056625 +#> HOVAL 0.0038497 0.441608 0.8975792 +# } +# \dontrun{ data(boston, package="spData") Wb <- as(spdep::nb2listw(boston.soi), "CsparseMatrix") trMatb <- trW(Wb, type="mult") @@ -603,7 +721,126 @@

    Examples

    data=boston.c, spdep::nb2listw(boston.soi), type="mixed", method="Matrix", control=list(fdHess=TRUE), 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, control = list(fdHess = TRUE)) +#> +#> 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.89816225 0.24400203 7.7793 7.327e-15 +#> CRIM -0.00571021 0.00093857 -6.0839 1.173e-09 +#> ZN 0.00069091 0.00051874 1.3319 0.1828931 +#> INDUS -0.00111343 0.00307380 -0.3622 0.7171778 +#> CHAS1 -0.04163225 0.02738839 -1.5201 0.1284937 +#> I(NOX^2) -0.01034950 0.19358633 -0.0535 0.9573639 +#> I(RM^2) 0.00794979 0.00102109 7.7856 6.883e-15 +#> AGE -0.00128789 0.00048929 -2.6322 0.0084838 +#> log(DIS) -0.12404108 0.09510145 -1.3043 0.1921304 +#> log(RAD) 0.05863502 0.02257529 2.5973 0.0093957 +#> TAX -0.00049084 0.00012146 -4.0412 5.317e-05 +#> PTRATIO -0.01319853 0.00595331 -2.2170 0.0266227 +#> B 0.00056383 0.00011084 5.0867 3.643e-07 +#> log(LSTAT) -0.24724454 0.02265149 -10.9152 < 2.2e-16 +#> lag.CRIM -0.00464215 0.00173900 -2.6694 0.0075978 +#> lag.ZN -0.00037937 0.00070703 -0.5366 0.5915659 +#> lag.INDUS 0.00025064 0.00385911 0.0649 0.9482165 +#> lag.CHAS1 0.12518252 0.04083949 3.0652 0.0021750 +#> lag.I(NOX^2) -0.38640403 0.22253428 -1.7364 0.0824967 +#> lag.I(RM^2) -0.00451252 0.00148919 -3.0302 0.0024440 +#> lag.AGE 0.00149678 0.00068470 2.1860 0.0288128 +#> lag.log(DIS) -0.00453785 0.10056961 -0.0451 0.9640104 +#> lag.log(RAD) -0.00940702 0.03127787 -0.3008 0.7636002 +#> lag.TAX 0.00041083 0.00017859 2.3004 0.0214237 +#> lag.PTRATIO 0.00060355 0.00789994 0.0764 0.9391011 +#> lag.B -0.00050781 0.00014107 -3.5996 0.0003187 +#> lag.log(LSTAT) 0.09846781 0.03574183 2.7550 0.0058697 +#> +#> Rho: 0.59578, LR test value: 181.68, p-value: < 2.22e-16 +#> Asymptotic standard error: 0.037474 +#> z-value: 15.899, p-value: < 2.22e-16 +#> Wald statistic: 252.76, 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 +#> summary(impacts(gp2mMi, tr=trMatb, R=1000), zstats=TRUE, short=TRUE) +#> Impact measures (mixed, trace): +#> Direct Indirect Total +#> CRIM -0.0074555753 -0.0181548470 -0.0256104223 +#> ZN 0.0006979073 0.0000727985 0.0007707058 +#> INDUS -0.0012029822 -0.0009314672 -0.0021344494 +#> CHAS1 -0.0198526431 0.2265453879 0.2066927447 +#> I(NOX^2) -0.0955268252 -0.8859909283 -0.9815177535 +#> I(RM^2) 0.0079983430 0.0005050175 0.0085033605 +#> AGE -0.0011296134 0.0016463650 0.0005167515 +#> log(DIS) -0.1410601685 -0.1770277708 -0.3180879393 +#> log(RAD) 0.0641735546 0.0576102594 0.1217838140 +#> TAX -0.0004651543 0.0002672119 -0.0001979424 +#> PTRATIO -0.0147737151 -0.0163846522 -0.0311583673 +#> B 0.0005265343 -0.0003879424 0.0001385920 +#> log(LSTAT) -0.2578403220 -0.1102143882 -0.3680547101 +#> ======================================================== +#> Simulation results ( variance matrix): +#> ======================================================== +#> Simulated standard errors +#> Direct Indirect Total +#> CRIM 0.0010076566 0.0036285971 0.0040401747 +#> ZN 0.0004896474 0.0012385892 0.0012617113 +#> INDUS 0.0028820705 0.0060789543 0.0057788737 +#> CHAS1 0.0271174416 0.0816948950 0.0902306553 +#> I(NOX^2) 0.1798525319 0.2954129272 0.2426994629 +#> I(RM^2) 0.0010995238 0.0031443424 0.0036252679 +#> AGE 0.0004766510 0.0012848631 0.0013546129 +#> log(DIS) 0.0884514383 0.1106832550 0.0757291955 +#> log(RAD) 0.0197637995 0.0513249209 0.0502885515 +#> TAX 0.0001172963 0.0003281755 0.0003385766 +#> PTRATIO 0.0054150426 0.0127887825 0.0122251822 +#> B 0.0001051363 0.0002388721 0.0002375770 +#> log(LSTAT) 0.0224736420 0.0610706224 0.0660453166 +#> +#> Simulated z-values: +#> Direct Indirect Total +#> CRIM -7.4242633 -5.04782172 -6.3852730 +#> ZN 1.4302178 0.09961643 0.6528326 +#> INDUS -0.4184439 -0.10883366 -0.3231737 +#> CHAS1 -0.7260661 2.79074703 2.3085362 +#> I(NOX^2) -0.5230875 -3.05728030 -4.1089450 +#> I(RM^2) 7.3426845 0.17905211 2.3822949 +#> AGE -2.4149822 1.25401260 0.3396769 +#> log(DIS) -1.5778895 -1.66514935 -4.2766960 +#> log(RAD) 3.2376719 1.17611638 2.4727850 +#> TAX -3.9726814 0.78133925 -0.6189573 +#> PTRATIO -2.7323855 -1.29137856 -2.5612005 +#> B 5.0000793 -1.62377140 0.5800914 +#> log(LSTAT) -11.4429357 -1.77731839 -5.5372039 +#> +#> Simulated p-values: +#> Direct Indirect Total +#> CRIM 1.1346e-13 4.4688e-07 1.7109e-10 +#> ZN 0.1526545 0.9206489 0.513864 +#> INDUS 0.6756226 0.9133344 0.746564 +#> CHAS1 0.4677982 0.0052587 0.020969 +#> I(NOX^2) 0.6009133 0.0022336 3.9747e-05 +#> I(RM^2) 2.0939e-13 0.8578968 0.017205 +#> AGE 0.0157360 0.2098374 0.734100 +#> log(DIS) 0.1145910 0.0958830 1.8969e-05 +#> log(RAD) 0.0012051 0.2395484 0.013406 +#> TAX 7.1068e-05 0.4346030 0.535944 +#> PTRATIO 0.0062878 0.1965724 0.010431 +#> B 5.7307e-07 0.1044246 0.561853 +#> log(LSTAT) < 2.22e-16 0.0755159 3.0734e-08 #data(house, package="spData") #lw <- spdep::nb2listw(LO_nb) #form <- formula(log(price) ~ age + I(age^2) + I(age^3) + log(lotsize) + @@ -621,7 +858,7 @@

    Examples

    #summary(mobj) #moobj <- impacts(mobj, tr=trMat, R=1000) #summary(moobj, zstats=TRUE, short=TRUE) -} +# }     @@ -636,7 +873,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/index.html b/docs/reference/index.html index edcf325..2ff40b5 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -1,5 +1,5 @@ -Function reference • spatialregFunction reference • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -137,8 +137,7 @@

    All functions

    predict(<Sarlm>) print(<Sarlm.pred>) as.data.frame(<Sarlm.pred>)

    -

    Prediction for spatial simultaneous autoregressive linear -model objects

    +

    Prediction for spatial simultaneous autoregressive linear model objects

    LR.Sarlm() logLik(<Sarlm>) LR1.Sarlm() Wald1.Sarlm() Hausman.test(<Sarlm>) anova(<Sarlm>) bptest.Sarlm() impacts(<Sarlm>)

    @@ -184,7 +183,7 @@

    All functions
    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/invIrM-1.png b/docs/reference/invIrM-1.png index 708c8ff..9bd2900 100644 Binary files a/docs/reference/invIrM-1.png and b/docs/reference/invIrM-1.png differ diff --git a/docs/reference/invIrM.html b/docs/reference/invIrM.html index 9b85742..38ed769 100644 --- a/docs/reference/invIrM.html +++ b/docs/reference/invIrM.html @@ -1,5 +1,5 @@ -Compute SAR generating operator — invIrM • spatialregCompute SAR generating operator — invIrM • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@

  • Changelog @@ -174,13 +174,21 @@

    Examples

    col=c("black", "brown", "green", "orange", "pink", "red"), lty=1, bty="n") } -if (FALSE) { +# \dontrun{ x <- matrix(rnorm(length(nb7rt)), ncol=1) system.time(e <- invIrM(nb7rt, rho=0.9, method="chol", feasible=TRUE) %*% x) +#> user system elapsed +#> 0.002 0.000 0.003 system.time(e <- invIrM(nb7rt, rho=0.9, method="chol", feasible=NULL) %*% x) +#> user system elapsed +#> 0.003 0.000 0.003 system.time(e <- invIrM(nb7rt, rho=0.9, method="solve", feasible=TRUE) %*% x) +#> user system elapsed +#> 0.001 0.000 0.001 system.time(e <- invIrM(nb7rt, rho=0.9, method="solve", feasible=NULL) %*% x) -} +#> user system elapsed +#> 0.001 0.000 0.002 +# } @@ -195,7 +203,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/lagmess.html b/docs/reference/lagmess.html index a9a2f18..a6e501d 100644 --- a/docs/reference/lagmess.html +++ b/docs/reference/lagmess.html @@ -1,5 +1,5 @@ -Matrix exponential spatial lag model — lagmess • spatialregMatrix exponential spatial lag model — lagmess • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -250,7 +250,7 @@

    Examples

    #> system.time(obj2 <- lagmess(log(PRICE) ~ PATIO + log(AGE) + log(SQFT), data=baltimore, listw=lw)) #> user system elapsed -#> 0.034 0.000 0.034 +#> 0.03 0.00 0.03 (x <- summary(obj2)) #> Matrix exponential spatial lag model: #> @@ -457,7 +457,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/lextrB.html b/docs/reference/lextrB.html index 3fa582a..34ece15 100644 --- a/docs/reference/lextrB.html +++ b/docs/reference/lextrB.html @@ -1,5 +1,5 @@ -Find extreme eigenvalues of binary symmetric spatial weights — lextrB • spatialregFind extreme eigenvalues of binary symmetric spatial weights — lextrB • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -172,22 +172,25 @@

    Examples

    eigs(B, k=1, which="LR")$values } #> [1] 1 -if (FALSE) { +# \dontrun{ ab.listw <- spdep::nb2listw(boston.soi, style="S") er <- range(eigenw(similar.listw(ab.listw))) er +#> [1] -0.723495 1.110373 res_1 <- lextrS(ab.listw) c(res_1) -} +#> lambda_n lambda_1 +#> -0.7230376 1.1103694 +# } if (run) { B <- as(similar.listw(ab.listw), "CsparseMatrix") eigs(B, k=1, which="SR")$values } -#> [1] -0.9708644 +#> [1] -0.723495 if (run) { eigs(B, k=1, which="LR")$values } -#> [1] 1 +#> [1] 1.110373 @@ -202,7 +205,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/predict.sarlm.html b/docs/reference/predict.sarlm.html index 589cad6..b0e694b 100644 --- a/docs/reference/predict.sarlm.html +++ b/docs/reference/predict.sarlm.html @@ -1,7 +1,5 @@ -Prediction for spatial simultaneous autoregressive linear -model objects — predict.Sarlm • spatialregPrediction for spatial simultaneous autoregressive linear model objects — predict.Sarlm • spatialregLikelihood ratio test — LR.Sarlm • spatialregLikelihood ratio test — LR.Sarlm • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -253,7 +253,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/set.mcOption.html b/docs/reference/set.mcOption.html index 581b9ae..6470e5a 100644 --- a/docs/reference/set.mcOption.html +++ b/docs/reference/set.mcOption.html @@ -1,5 +1,5 @@ -Options for parallel support — set.mcOption • spatialregOptions for parallel support — set.mcOption • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -126,7 +126,7 @@

    Examples

    if(.Platform$OS.type == "windows") { # forking not permitted on Windows - start cluster # removed for Github actions 210502 -if (FALSE) { +# \dontrun{ print(get.mcOption()) cl <- makeCluster(get.coresOption()) print(clusterEvalQ(cl, exists("mom_calc"))) @@ -137,7 +137,7 @@

    Examples

    set.ClusterOption(NULL) print(clusterEvalQ(cl, exists("mom_calc"))) stopCluster(cl) -} +# } } else { mcOpt <- get.mcOption() print(mcOpt) @@ -185,7 +185,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/set.spChkOption.html b/docs/reference/set.spChkOption.html index a1d3df4..edeef26 100644 --- a/docs/reference/set.spChkOption.html +++ b/docs/reference/set.spChkOption.html @@ -1,5 +1,5 @@ -Control checking of spatial object IDs — set.ZeroPolicyOption • spatialregControl checking of spatial object IDs — set.ZeroPolicyOption • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -117,7 +117,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/similar.listw.html b/docs/reference/similar.listw.html index 420111a..853df67 100644 --- a/docs/reference/similar.listw.html +++ b/docs/reference/similar.listw.html @@ -1,5 +1,5 @@ -Create symmetric similar weights lists — similar.listw • spatialregCreate symmetric similar weights lists — similar.listw • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -132,7 +132,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/sparse_mat.html b/docs/reference/sparse_mat.html index 7789c86..427107a 100644 --- a/docs/reference/sparse_mat.html +++ b/docs/reference/sparse_mat.html @@ -1,5 +1,5 @@ -Spatial neighbour sparse representation — as.spam.listw • spatialregSpatial neighbour sparse representation — as.spam.listw • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -115,7 +115,7 @@

    See also

    Examples

    - 
    if (FALSE) {
+    
    # \dontrun{
 require(sf, quietly=TRUE)
 columbus <- st_read(system.file("shapes/columbus.shp", package="spData")[1], quiet=TRUE)
 #require(spdep, quietly=TRUE)
@@ -125,9 +125,28 @@ 
    Examples
    
   col.sp <- as.spam.listw(col.listw)
   str(col.sp)
 }
+#> Spam version 2.10-0 (2023-10-23) 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)
+#> Warning: style is M (missing); style should be set to a valid value
 listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B")
 W_C <- as(listw_NY, "CsparseMatrix")
 W_R <- as(listw_NY, "RsparseMatrix")
@@ -136,11 +155,14 @@ 
    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)
+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)
@@ -158,33 +180,46 @@ 
    Examples
    
 W <- as(lw, "CsparseMatrix")
 system.time(e <- invIrM(nb7rt, rho=0.98, method="solve", feasible=NULL) %*% x)
 #>    user  system elapsed 
-#>   0.003   0.000   0.002 
+#>   0.002   0.000   0.002 
 system.time(ee <- powerWeights(W, rho=0.98, X=x))
 #> Warning: not converged within order iterations
 #>    user  system elapsed 
-#>   0.182   0.004   0.186 
+#>   0.181   0.010   0.192 
 str(attr(ee, "internal"))
 #> List of 5
-#>  $ series: num [1:250] 0.286 0.233 0.199 0.175 0.157 ...
+#>  $ series: num [1:250] 0.287 0.234 0.201 0.178 0.16 ...
 #>  $ 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.00604604"
-if (FALSE) {
+#> [1] "Mean relative difference: 0.0060747"
+# \dontrun{
 system.time(ee <- powerWeights(W, rho=0.9, X=x))
+#>    user  system elapsed 
+#>   0.133   0.005   0.140 
 system.time(ee <- powerWeights(W, rho=0.98, order=1000, X=x))
+#>    user  system elapsed 
+#>   0.769   0.008   0.782 
 all.equal(e, as(ee, "matrix"), check.attributes=FALSE)
+#> [1] TRUE
 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.009   0.000   0.009 
 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")
+#> Error in errorsarlm(as(ee, "matrix")[, 1] ~ 1, listw = spdep::nb2listw(nb60rt),     method = "Matrix"): argument "Durbin" is missing, with no default
 coefficients(obj)
-}
+#> Error in eval(expr, envir, enclos): object 'obj' not found
+# }
    @@ -199,7 +234,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/spautolm-1.png b/docs/reference/spautolm-1.png deleted file mode 100644 index 629b736..0000000 Binary files a/docs/reference/spautolm-1.png and /dev/null differ diff --git a/docs/reference/spautolm.html b/docs/reference/spautolm.html index f61e39c..9f7edce 100644 --- a/docs/reference/spautolm.html +++ b/docs/reference/spautolm.html @@ -1,5 +1,5 @@ -Spatial conditional and simultaneous autoregression model estimation — spautolm • spatialregSpatial conditional and simultaneous autoregression model estimation — spautolm • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -329,12 +329,55 @@

    See also

    Examples

    require("sf", quietly=TRUE)
 nydata <- st_read(system.file("shapes/NY8_bna_utm18.gpkg", package="spData")[1], quiet=TRUE)
-if (FALSE) {
+# \dontrun{
 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(
@@ -342,44 +385,263 @@ 
    Examples
    
 identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10))
 #> [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")
+listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), style="B")
 eigs <- eigenw(listw_NY)
-if (FALSE) {
+# \dontrun{
 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.196   0.000   0.197 
 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.017201 
+#> 
+#> 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.126350e-02
+#> PCTAGE65P    3.75419996 0.62472153  6.009397 1.862142e-09
+#> PCTOWNHOME  -0.41988960 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.01720109 0.18504698 0.04378056 0.62990270 0.20343911
 system.time(esar1M <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
  data=nydata, listw=listw_NY, family="SAR", method="Matrix"))
+#>    user  system elapsed 
+#>   0.221   0.000   0.221 
 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.196   0.000   0.196 
 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.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: NA (not available for weighted model)
+#> 
 system.time(esar1wM <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
  data=nydata, listw=listw_NY, weights=POP8, family="SAR", method="Matrix"))
+#>    user  system elapsed 
+#>   0.212   0.000   0.214 
 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: NA (not available for weighted model)
+#> 
 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: NA (not available for weighted model)
+#> 
 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: NA (not available for weighted model)
+#> 
+# }
 ecar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata,
  listw=listw_NY, family="CAR", method="eigen",
  control=list(pre_eig=eigs))
@@ -409,11 +671,38 @@ 
    Examples
    
 #> Number of parameters estimated: 6 
 #> AIC: 563.66
 #> 
-if (FALSE) {
+# \dontrun{
 system.time(ecar1M <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
  data=nydata, listw=listw_NY, family="CAR", method="Matrix"))
+#>    user  system elapsed 
+#>   0.236   0.000   0.237 
 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))
@@ -444,12 +733,39 @@ 
    Examples
    
 #> Number of parameters estimated: 6 
 #> AIC: NA (not available for weighted model)
 #> 
-if (FALSE) {
+# \dontrun{
 system.time(ecar1wM <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
  data=nydata, listw=listw_NY, weights=POP8, family="CAR", method="Matrix"))
+#>    user  system elapsed 
+#>   0.243   0.000   0.244 
 summary(ecar1wM)
-}
-if (FALSE) {
+#> 
+#> 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: NA (not available for weighted model)
+#> 
+# }
+# \dontrun{
 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) +
@@ -466,63 +782,179 @@ 
    Examples
    
 el2 <- sqrt(n[sids.nhbr$to]/n[sids.nhbr$from])
 sids.nhbr$weights <- el1*el2
 sids.nhbr.listw <- spdep::sn2listw(sids.nhbr)
+#> Warning: style is M (missing); style should be set to a valid value
+#> Warning: 56, 87 are not origins
+#> Error in spdep::sn2listw(sids.nhbr): no-neighbour observations found, set zero.policy to TRUE
 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)
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 W.4 <- W[-outl, -outl]
+#> Error in eval(expr, envir, enclos): object 'W' not found
 sids.nhbr.listw.4 <- spdep::mat2listw(W.4)
+#> Error in eval(expr, envir, enclos): object 'W.4' not found
 esarI <- errorsarlm(ft.SID74 ~ 1, data=mdata, listw=sids.nhbr.listw,
  zero.policy=TRUE)
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarI)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarI' not found
 esarIa <- spautolm(ft.SID74 ~ 1, data=mdata, listw=sids.nhbr.listw,
  family="SAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarIa)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarIa' not found
 esarIV <- errorsarlm(ft.SID74 ~ ft.NWBIR74, data=mdata, listw=sids.nhbr.listw,
  zero.policy=TRUE)
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarIV)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarIV' not found
 esarIVa <- spautolm(ft.SID74 ~ ft.NWBIR74, data=mdata, listw=sids.nhbr.listw,
  family="SAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarIVa)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarIVa' not found
 esarIaw <- spautolm(ft.SID74 ~ 1, data=mdata, listw=sids.nhbr.listw,
  weights=BIR74, family="SAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarIaw)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarIaw' not found
 esarIIaw <- spautolm(ft.SID74 ~ both - 1, data=mdata, listw=sids.nhbr.listw,
  weights=BIR74, family="SAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarIIaw)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarIIaw' not found
 esarIVaw <- spautolm(ft.SID74 ~ ft.NWBIR74, data=mdata,
  listw=sids.nhbr.listw, weights=BIR74, family="SAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw' not found
 summary(esarIVaw)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'esarIVaw' not found
 ecarIaw <- spautolm(ft.SID74 ~ 1, data=mdata.4, listw=sids.nhbr.listw.4,
  weights=BIR74, family="CAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw.4' not found
 summary(ecarIaw)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'ecarIaw' not found
 ecarIIaw <- spautolm(ft.SID74 ~ both - 1, data=mdata.4,
  listw=sids.nhbr.listw.4, weights=BIR74, family="CAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw.4' not found
 summary(ecarIIaw)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'ecarIIaw' not found
 ecarIVaw <- spautolm(ft.SID74 ~ ft.NWBIR74, data=mdata.4,
  listw=sids.nhbr.listw.4, weights=BIR74, family="CAR")
+#> Error in eval(expr, envir, enclos): object 'sids.nhbr.listw.4' not found
 summary(ecarIVaw)
+#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'ecarIVaw' not found
 nc.sids$fitIV <- append(fitted.values(ecarIVaw), NA, outl-1)
+#> Error in eval(expr, envir, enclos): object 'ecarIVaw' not found
 plot(nc.sids[,"fitIV"], nbreaks=12) # Cressie 1993, p. 565
-}
-if (FALSE) {
+#> Error in `[.data.frame`(x, i, j, drop = drop): undefined columns selected
+# }
+# \dontrun{
 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.366162 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.15241 
+#> 
+#> 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.044296 
+#> 
+#> 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
+#> 
+# }
    @@ -537,7 +969,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/stsls.html b/docs/reference/stsls.html index eed1125..8cb5c1c 100644 --- a/docs/reference/stsls.html +++ b/docs/reference/stsls.html @@ -1,5 +1,5 @@ -Generalized spatial two stage least squares — stsls • spatialregGeneralized spatial two stage least squares — stsls • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -262,19 +262,19 @@

    Examples

    #> Simulation results ( variance matrix): #> ======================================================== #> Simulated standard errors -#> Direct Indirect Total -#> INC 0.4250637 0.7860257 0.9598058 -#> HOVAL 0.1068486 0.6601511 0.7218709 +#> Direct Indirect Total +#> INC 0.39668006 0.9695180 1.0529050 +#> HOVAL 0.09757127 0.3123932 0.3649017 #> #> Simulated z-values: -#> Direct Indirect Total -#> INC -2.474018 -1.0625140 -1.9657924 -#> HOVAL -2.826222 -0.4999244 -0.8755077 +#> Direct Indirect Total +#> INC -2.768354 -0.8583926 -1.833383 +#> HOVAL -2.954457 -0.8425926 -1.511339 #> #> Simulated p-values: #> Direct Indirect Total -#> INC 0.0133603 0.28800 0.049323 -#> HOVAL 0.0047101 0.61713 0.381298 +#> INC 0.0056340 0.39068 0.066746 +#> HOVAL 0.0031322 0.39946 0.130702 ev <- eigenw(lw) loobj2 <- impacts(COL.lag.stsls, R=200, evalues=ev) summary(loobj2, zstats=TRUE, short=TRUE) @@ -287,23 +287,23 @@

    Examples

    #> ======================================================== #> Simulated standard errors #> Direct Indirect Total -#> INC 0.3647280 1.2688013 1.3886153 -#> HOVAL 0.1175608 0.8263993 0.9011554 +#> INC 0.3841092 0.7160775 0.8241613 +#> HOVAL 0.1008357 0.4126274 0.4698167 #> #> Simulated z-values: -#> Direct Indirect Total -#> INC -2.955756 -0.7150236 -1.4296758 -#> HOVAL -2.404673 -0.4145933 -0.6939034 +#> Direct Indirect Total +#> INC -2.895776 -1.1602746 -2.357719 +#> HOVAL -2.743940 -0.6064747 -1.121576 #> #> Simulated p-values: -#> Direct Indirect Total -#> INC 0.003119 0.47459 0.15281 -#> HOVAL 0.016187 0.67844 0.48774 +#> Direct Indirect Total +#> INC 0.0037822 0.24594 0.018388 +#> HOVAL 0.0060707 0.54420 0.262043 require(coda) HPDinterval(loobj1) -#> lower upper -#> INC -1.8854657 -0.2997080 -#> HOVAL -0.4662763 -0.1121351 +#> lower upper +#> INC -1.7799556 -0.35225230 +#> HOVAL -0.4587508 -0.09515137 #> attr(,"Probability") #> [1] 0.95 COL.lag.stslsW <- stsls(CRIME ~ INC + HOVAL, data=COL.OLD, lw, W2X=FALSE) @@ -431,7 +431,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/docs/reference/trW-1.png b/docs/reference/trW-1.png index d441164..9871e3b 100644 Binary files a/docs/reference/trW-1.png and b/docs/reference/trW-1.png differ diff --git a/docs/reference/trW.html b/docs/reference/trW.html index 86baebb..77af24b 100644 --- a/docs/reference/trW.html +++ b/docs/reference/trW.html @@ -1,5 +1,5 @@ -Spatial weights matrix powers traces — trW • spatialregSpatial weights matrix powers traces — trW • spatialreg @@ -17,7 +17,7 @@ spatialreg - 1.3-2 + 1.3-3 @@ -32,14 +32,14 @@
  • Changelog @@ -140,17 +140,17 @@

    Examples

    W <- as(listw, "CsparseMatrix") system.time(trMat <- trW(W, type="mult")) #> user system elapsed -#> 0.002 0.000 0.003 +#> 0.003 0.000 0.003 str(trMat) #> num [1:30] 0 10.91 3.65 5.62 3.66 ... -#> - attr(*, "timings")= Named num [1:2] 0.002 0.003 +#> - attr(*, "timings")= Named num [1:2] 0.003 0.003 #> ..- attr(*, "names")= chr [1:2] "user.self" "elapsed" #> - attr(*, "type")= chr "mult" #> - attr(*, "n")= int 49 set.seed(1100) system.time(trMC <- trW(W, type="MC")) #> user system elapsed -#> 0.005 0.001 0.005 +#> 0.005 0.000 0.006 str(trMC) #> num [1:30] 0 10.91 3.69 5.36 3.64 ... #> - attr(*, "sd")= num [1:30] NA NA 0.598 0.495 0.489 ... @@ -166,13 +166,13 @@

    Examples

    } listwS <- similar.listw(listw) -W <- forceSymmetric(as(listwS, "CsparseMatrix")) +W <- forceSymmetric(as(listwS, "CsparseMatrix")) system.time(trmom <- trW(listw=listwS, m=24, type="moments")) #> user system elapsed -#> 0.001 0.000 0.002 +#> 0.002 0.000 0.001 str(trmom) #> num [1:24] 0 10.91 3.65 5.62 3.66 ... -#> - attr(*, "timings")= Named num [1:2] 0.001 0.002 +#> - attr(*, "timings")= Named num [1:2] 0.002 0.001 #> ..- attr(*, "names")= chr [1:2] "user.self" "elapsed" #> - attr(*, "type")= chr "moments" #> - attr(*, "n")= int 49 @@ -180,10 +180,10 @@

    Examples

    #> [1] TRUE system.time(trMat <- trW(W, m=24, type="mult")) #> user system elapsed -#> 0.004 0.000 0.004 +#> 0.003 0.000 0.003 str(trMat) #> num [1:24] 0 10.91 3.65 5.62 3.66 ... -#> - attr(*, "timings")= Named num [1:2] 0.003 0.004 +#> - attr(*, "timings")= Named num [1:2] 0.003 0.003 #> ..- attr(*, "names")= chr [1:2] "user.self" "elapsed" #> - attr(*, "type")= chr "mult" #> - attr(*, "n")= int 49 @@ -192,20 +192,27 @@

    Examples

    set.seed(1) system.time(trMC <- trW(W, m=24, type="MC")) #> user system elapsed -#> 0.006 0.000 0.006 +#> 0.006 0.000 0.005 str(trMC) #> num [1:24] 0 10.91 2.44 4.97 2.82 ... #> - attr(*, "sd")= num [1:24] NA NA 0.618 0.501 0.451 ... -#> - attr(*, "timings")= Named num [1:2] 0.006 0.006 +#> - attr(*, "timings")= Named num [1:2] 0.006 0.005 #> ..- attr(*, "names")= chr [1:2] "user.self" "elapsed" #> - attr(*, "type")= chr "MC" #> - attr(*, "n")= int 49 -if (FALSE) { +# \dontrun{ data(boston, package="spData") listw <- spdep::nb2listw(boston.soi) listwS <- similar.listw(listw) system.time(trmom <- trW(listw=listwS, m=24, type="moments")) +#> user system elapsed +#> 0.108 0.000 0.108 str(trmom) +#> num [1:24] 0 124.2 32.7 63.7 33.2 ... +#> - attr(*, "timings")= Named num [1:2] 0.108 0.108 +#> ..- attr(*, "names")= chr [1:2] "user.self" "elapsed" +#> - attr(*, "type")= chr "moments" +#> - attr(*, "n")= int 506 library(parallel) nc <- max(2L, detectCores(logical=FALSE), na.rm = TRUE)-1L # set nc to 1L here @@ -217,13 +224,16 @@

    Examples

    set.ClusterOption(cl) } system.time(trmomp <- trW(listw=listwS, m=24, type="moments")) +#> user system elapsed +#> 0.111 0.000 0.112 if(!get.mcOption()) { set.ClusterOption(NULL) stopCluster(cl) } all.equal(trmom, trmomp, check.attributes=FALSE) +#> [1] TRUE invisible(set.coresOption(coresOpt)) -} +# }     @@ -238,7 +248,7 @@

    Examples

    -

    Site built with pkgdown 2.0.7.

    +

    Site built with pkgdown 2.0.9.

    diff --git a/man/GMerrorsar.Rd b/man/GMerrorsar.Rd index 7fe716f..d3349b8 100644 --- a/man/GMerrorsar.Rd +++ b/man/GMerrorsar.Rd @@ -125,8 +125,7 @@ suppressMessages(nyadjmat <- as.matrix(foreign::read.dbf(system.file( suppressMessages(ID <- as.character(names(foreign::read.dbf(system.file( "misc/nyadjwts.dbf", package="spData")[1]))[-1])) identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10)) -nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY)) -listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B") +listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") summary(esar1f) diff --git a/man/MCMCsamp.Rd b/man/MCMCsamp.Rd index 4c392a9..69a2c8e 100644 --- a/man/MCMCsamp.Rd +++ b/man/MCMCsamp.Rd @@ -37,8 +37,7 @@ suppressMessages(ID <- as.character(names(foreign::read.dbf(system.file( "misc/nyadjwts.dbf", package="spData")[1]))[-1])) identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10)) #require("spdep", quietly=TRUE) -nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY)) -listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B") +listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), style="B") esar1f <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, listw=listw_NY, family="SAR", method="eigen") summary(esar1f) diff --git a/man/spautolm.Rd b/man/spautolm.Rd index 1430192..60f4b4f 100644 --- a/man/spautolm.Rd +++ b/man/spautolm.Rd @@ -143,8 +143,7 @@ suppressMessages(ID <- as.character(names(foreign::read.dbf(system.file( "misc/nyadjwts.dbf", package="spData")[1]))[-1])) identical(substring(ID, 2, 10), substring(as.character(nydata$AREAKEY), 2, 10)) #require("spdep", quietly=TRUE) -nyadjlw <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY)) -listw_NY <- spdep::nb2listw(nyadjlw$neighbours, style="B") +listw_NY <- spdep::mat2listw(nyadjmat, as.character(nydata$AREAKEY), style="B") eigs <- eigenw(listw_NY) \dontrun{ esar0 <- errorsarlm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data=nydata, diff --git a/vignettes/nb_igraph.Rmd b/vignettes/nb_igraph.Rmd index 27e6a86..bd680a5 100644 --- a/vignettes/nb_igraph.Rmd +++ b/vignettes/nb_igraph.Rmd @@ -102,7 +102,8 @@ rownames(B)[1: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: ```{r, echo=dothis, eval=dothis} -nb_B1 <- spdep::mat2listw(as(as(B, "TsparseMatrix"), "CsparseMatrix")) +nb_B1 <- spdep::mat2listw(as(as(B, "TsparseMatrix"), "CsparseMatrix"), + style="B", zero.policy=TRUE) nb_B1$style all.equal(nb_B1$neighbours, col2, check.attributes=FALSE) all.equal(attr(nb_B1$neighbours, "region.id"), attr(nb_B$neighbours, "region.id")) @@ -133,7 +134,7 @@ The computation of a sparse Cholesky decomposition for each value of the spatial ```{r, echo=dothis, eval=dothis} nW <- -B nChol <- Cholesky(nW, Imult=8) -n * log(rho) + (2 * c(determinant(update(nChol, nW, 1/rho))$modulus)) +n * log(rho) + (2 * c(determinant(update(nChol, nW, 1/rho), sqrt=TRUE)$modulus)) ``` ### Asymmetric sparse matrices @@ -231,7 +232,7 @@ First we'll see how to get from sparse matrices to graphs. The mode of a symmetr class(B) object.size(B) if (!require("igraph", quietly=FALSE)) dothis <- FALSE -g1 <- graph.adjacency(B, mode="undirected") +g1 <- graph_from_adjacency_matrix(B, mode="undirected") class(g1) object.size(g1) ``` @@ -243,7 +244,7 @@ We can also convert this graph pack to the same matrix, but note that `get.adjac # Matrix 1.4-2 vulnerability work-around ow <- options("warn")$warn options("warn"=2L) -B1 <- try(get.adjacency(g1), silent=TRUE) +B1 <- try(as_adjacency_matrix(g1), silent=TRUE) if (!inherits(B1, "try-error")) print(class(B1)) if (!inherits(B1, "try-error")) print(object.size(B1)) if (!inherits(B1, "try-error")) print(all.equal(B, as(B1, "CsparseMatrix"))) @@ -264,7 +265,7 @@ table(res$comp.id) The same result can be obtained using the `clusters` function in **igraph**: ```{r, echo=dothis, eval=dothis} -c1 <- clusters(g1) +c1 <- components(g1) c1$no == res$nc all.equal(c1$membership, res$comp.id) all.equal(c1$csize, c(table(res$comp.id)), check.attributes=FALSE) @@ -274,8 +275,8 @@ The same holds for the row-standardised variant: ```{r, echo=dothis, eval=dothis} W <- as(spdep::nb2listw(col2, style="W", zero.policy=TRUE), "CsparseMatrix") -g1W <- graph.adjacency(W, mode="directed", weighted="W") -c1W <- clusters(g1W) +g1W <- graph_from_adjacency_matrix(W, mode="directed", weighted="W") +c1W <- components(g1W) all.equal(c1W$membership, res$comp.id) ``` ### Shortest paths in weights matrices: **igraph** @@ -283,10 +284,10 @@ all.equal(c1W$membership, res$comp.id) 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 an `n` x `n` matrix: ```{r, echo=dothis, eval=dothis} -is.connected(g1) +is_connected(g1) dg1 <- diameter(g1) dg1 -sp_mat <- shortest.paths(g1) +sp_mat <- distances(g1) str(sp_mat) ```