Re: [R-sig-Geo] Learning Resources Spatial Regression Models from the ground up

[Dexter Locke]

in GeoDA  Luc Anselin and Sergio J. Rey. (2014). Modern Spatial
Econometrics in Practice: A Guide to GeoDa, GeoDaSpace and PySAL. [link to
book]


Exploring Spatial Data with GeoDa: A Workbook (2005; 244 pp.,5.1Mb)
https://geodacenter.github.io/docs/geodaworkbook.pdfhttps://geodacenter.github.io/docs/geodaworkbook.pdf
Chapters
17, 18, 21, 22 - 25.

More conceptually advanced, at least look at figure 1: Golgher, A. B., &
Voss, P. R. (2016). How to Interpret the Coefficients of Spatial Models:
Spillovers, Direct and Indirect Effects. Spatial Demography, 4(3), 175–205.
https://doi.org/10.1007/s40980-015-0016-y

-Dexter


On Wed, Apr 24, 2024 at 12:03 PM Josiah Parry 
wrote:

> Thank you, Chris! I can take a look at the first resource. At the moment my
> interest is specifically in spatial econometric models and less so about
> point patterns (for the time being).
>
> On Wed, Apr 24, 2024 at 11:51 AM Christopher W. Ryan  >
> wrote:
>
> > Josiah--
> >
> > I've found the following very helpful over the years:
> >
> > Geographic Information Analysis, by David O'Sullivan and David Unwin
> >
> > Spatial Point Patterns, by Adrian Baddeley, Ege Rubak, and Rolf Turner
> >
> > Applied Spatial Data Analysis with R, by Roger Bivand, Edzer Pebesma,
> > and Virgilio Gomez-Rubio
> >
> > Statistical Analysis of Spatial and Spatio-Temporal Point Patterns
> >
> > The last 3 are, as the titles imply, focused specifically on spatial
> > point patterns. The first is a bit more general, including methods for
> > areal data.
> >
> > I listed them in increasing order (in my opinion) of mathemtical
> > complexity.
> >
> > --Chris Ryan
> >
> > In
> > Josiah Parry wrote:
> > > Hey folks,
> > >
> > > I'm hoping to build up my knowledge around spatial regression
> techniques
> > > from the ground up—e.g. I'm not interested in R-INLA or other
> > exceptionally
> > > complex techniques.
> > >
> > > I'm hoping this listserv has some recommendations for what readings /
> > > models I should prioritize learning about in, possibly, an opinionated
> > > order.
> > >
> > > At the moment I've purchased "Modern Spatial Econometrics in Practice"
> by
> > > Luc Anselin and Sergio Rey and will try to work through that. But if
> > there
> > > are additional resources that folks recommend that are friendly for the
> > > not-so-math-inclined, I'd love to have a look at them!
> > >
> > > The Spatial Regression section of the R-spatial book (
> > > https://r-spatial.org/book/16-SpatialRegression.html) is good but with
> > less
> > > handholding than I might need.
> > >
> > >   [[alternative HTML version deleted]]
> > >
> > > ___
> > > R-sig-Geo mailing list
> > > R-sig-Geo@r-project.org
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
> >
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Recommendations for wildfire risk analysis (spatially explicit)

[Dexter Locke]

Possibly helpful:
https://www.publish.csiro.au/WF/WF21176
https://www.fs.usda.gov/research/treesearch/67085

-Dexter





On Sun, Jan 21, 2024 at 9:36 PM Ben Tupper  wrote:

> Hi,
>
> Perhaps this might be a good place to start?
>
> https://esajournals.onlinelibrary.wiley.com/doi/10.1002/eap.1898  The
> authors state that the work was done primarily in R.
>
> Cheers,
> Ben
>
>
> On Sat, Jan 20, 2024 at 9:43 AM Manuel Spínola 
> wrote:
> >
> > Dear list members,
> >
> > I am looking for suggestions for R packages or an approach using R for
> > wildfire risk analysis.
> >
> > Manuel
> >
> > --
> > *Manuel Spínola, Ph.D.*
> > Instituto Internacional en Conservación y Manejo de Vida Silvestre
> > Universidad Nacional
> > Apartado 1350-3000
> > Heredia
> > COSTA RICA
> > mspin...@una.cr 
> > mspinol...@gmail.com
> > Teléfono: (506) 8706 - 4662
> > Sitio web institucional: ICOMVIS
> > 
> > Sitio web personal: Sitio personal <
> https://mspinola-sitioweb.netlify.app>
> > Blog sobre Ciencia de Datos: Blog de Ciencia de Datos
> > 
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Counting points in cells and ploting only cells containing points

[Dexter Locke]

Here is one approach. Kind of a silly example with just one grid cell, but
would work with more polygons and more points.


library(sf)
library(tidyverse)

# make some points data
# modified example from
# https://r-spatial.github.io/sf/reference/geos_binary_pred.html
(pts <-
st_sfc(st_point(c(.5,.5)), st_point(c(1.5, 1.5)), st_point(c(2.5,
2.5))) |>
st_as_sf()|>
rowid_to_column(var = 'pts_id'))

(pol <-
st_polygon(list(rbind(c(0,0), c(2,0), c(2,2), c(0,2), c(0,0 |>
st_sfc() |>
st_as_sf()|>
rowid_to_column(var = 'pol_id')
  )

# look at the data, crude 'map'
plot(pts); plot(pol, add = TRUE) # take a look

# perform the intersection
pts_pol_int <-
  pts |>
  st_intersection(pol) # only going to retain the overlaps (drops
non-intersecting elements)

# count the overlaps
cont_pts_pol_int <-
  pts_pol_int |>
  st_drop_geometry() |> # pulls out the data frame (or tibble)
  group_by(pol_id) |>
  count()

cont_pts_pol_int # could be joined back to pol

HTH, DHL


On Wed, Nov 22, 2023 at 6:39 AM MIGUEL GOMEZ DE ANTONIO 
wrote:

> Dear community,
>
> I want to plot only the cells of a grid that contains points.
>
> I have a set of points (of class SpatialPointsDataFrame) and a grid (of
> class SpatialPolygonDataFrame).
>
> I am interested in ploting the set of cells of the grid that contains
> points:
>
> 1.   Count how many points fall in each cell of the grid.
>
> 2.   Plot only the cells that contains points.
>
> I tried:
>
> library(sf)
>
> Points.sf=st_as_sf(points)
>
> Grid.sf=st_as_sf(grid)
>
> A=intersects(grid.sf, points.sf)
>
> apply(A,1,any)
>
> There I get the cells that contains points but:
>
> 1. How can I count how many points fall in each cell
>
> 2. How could I plot the set of cells that contains points?
>
>
> Thank you very much for your help,
>
>
> Miguel Gómez de Antonio
> Profesor Titular de Universidad
> Dpto. Economía Aplicada, Pública y Política
> Universidad Complutense de Madrid
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Calculating median age for a group of US census blocks?

[Dexter Locke]

Hi Kevin and all,

Given the binned data, you could count the number of people per age class
for those 10 blocks. You can then express that in a number of
different ways, like percent under 25 years old, or by calculating the
dependency
ratio

.

I do think it is feasible to calculate an estimated mean from the counts
within groups representing ranges. See, for example, here:
https://stackoverflow.com/questions/18887382/how-to-calculate-the-median-on-grouped-dataset

Since you are working in Baltimore, you may consider looking at The
Baltimore Neighborhood Indicators Alliance https://bniajfi.org/vital_signs/.
They provide useful data on a range of issues (transportation, crime,
education, environment etc.) including summaries from Census-derived
demographics. What you are seeking may already exist. BNIA creates
neighborhoods or "community statistical areas" (n=55) based on aggregates
of Census data.

Although not pertaining to age, Baltimore City Planning has paid Census in
the past to aggregate from individual-level Census data to the more
colloquially-used definitions of Baltimore shown here (n = 273):
https://data.baltimorecity.gov/datasets/neighborhood-1/explore?location=39.284832%2C-76.620516%2C12.91

Best, Dexter
https://dexterlocke.com/





On Mon, Aug 7, 2023 at 3:02 PM Kevin Zembower via R-sig-Geo <
r-sig-geo@r-project.org> wrote:

> Josiah, thanks for your reply.
>
> Regarding my objective, I'm trying to compile census statistics for the
> blocks that make up the neighborhood where I live. It consists of ten
> census blocks, of which I selected three for simplicity in my example.
> The census block-group which contains these ten blocks also contains
> some blocks which are outside of my neighborhood and shouldn't be
> counted or included.
>
> Since I won't be able to calculate the median age from the age and count
> data, and since the individual data doesn't seem to be available, is it
> your thought that I can't produce a valid median age for a group of
> census blocks?
>
> Thanks so much for your advice.
>
> -Kevin
>
> On 8/7/23 14:38, Josiah Parry wrote:
> > Hey Kevin, I don't think you're going to be able to get individual level
> > data from the US Census Bureau. The closest you may be able to get is
> > the current population survey (CPS) which I believe is also available
> > via tidycensus. Regarding your first question, I'm not sure I follow
> > what your objective is with it. I would use a geography of census block
> > groups as the measure of median for census block groups. Otherwise it is
> > unclear how you are defining what a "group of blocks" is.
> >
> > On Mon, Aug 7, 2023 at 2:34 PM Kevin Zembower via R-sig-Geo
> > mailto:r-sig-geo@r-project.org>> wrote:
> >
> > Hello, all,
> >
> > I'd like to obtain the median age for a population in a specific
> group
> > of US Decennial census blocks. Here's an example of the problem:
> >
> > ## Example of calculating median age of population in census blocks.
> > library(tidyverse)
> > library(tidycensus)
> >
> > counts <- get_decennial(
> >   geography = "block",
> >   state = "MD",
> >   county = "Baltimore city",
> >   table = "P1",
> >   year = 2020,
> >   sumfile = "dhc") %>%
> >   mutate(NAME = NULL) %>%
> >   filter(substr(GEOID, 6, 11) == "271101" &
> >  substr(GEOID, 12, 15) %in% c(3000, 3001, 3002)
> >  )
> >
> > ages <- get_decennial(
> >   geography = "block",
> >   state = "MD",
> >   county = "Baltimore city",
> >   table = "P13",
> >   year = 2020,
> >   sumfile = "dhc") %>%
> >   mutate(NAME = NULL) %>%
> >   filter(substr(GEOID, 6, 11) == "271101" &
> >  substr(GEOID, 12, 15) %in% c(3000, 3001, 3002)
> >  )
> >
> > I have two questions:
> >
> > 1. Is it mathematically valid to multiply the population of a block
> by
> > the median age of that block (in other words, assign the median age
> to
> > each member of a block), then calculate the median of those numbers
> for
> > a group of blocks?
> >
> > 2. Is raw data on the ages of individuals available anywhere else in
> > the
> > census data? I can find tables such as P12, that breaks down the
> > population by age ranges or bins, but can't find specific data of
> > counts
> > per age in years.
> >
> > Thanks for your advice and help.
> >
> > -Kevin
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org 
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > 
> >
>
>
> ___
> 

Re: [R-sig-Geo] Calculate areas of multiple polygones within multiple polygones

[Dexter Locke]

Hi Basile,

I'm not sure I fully understand the question. How does a shapefile of the
buildings not have geometry?

Do you want the area of building per neighborhood, their counts, or both?

bat_neigh_int <- bat %>%
  mutate(bat_area = st_area(.)) %>% # adds area to input file, may want to
change units.
  st_join(. # add neighborhood identifier
"CODE_NEIGHBOURHOOD" via spatial join
  , neighborhoods
  , join = st_intersects # check the help, you may want st_within
or to try both ways
  , left = TRUE)

buil_area_per_neigh <- bat_neigh_int %>%
  st_drop_geometry() %>%
  group_by(CODE_NEIGHBOURHOOD) %>%
  summarise(neigh_build_area = sum(bat_area) # gets the area of buildings
in each neighborhood
, build_count = sum(n()))  # counts the
buildings in each neighborhood

Then join buil_area_per_neigh back to neighborhoods with left_join

Does this help? Again I am a little confused about the intended outcome.

Best, Dexter


On Mon, Dec 20, 2021 at 6:45 AM basile.lenormand via R-sig-Geo <
r-sig-geo@r-project.org> wrote:

> Hello every one!
>
> I am trying to compute the area of buildings inside neighbourhoods. I have
> a shapefile with the buildings and one with the neighbourhoods. I want to
> know the built-up area per neighbourhood but I do not understand how to
> procede.
> I already compute the areas of each neighbourhood and I reach to attribute
> at each neighbourhood the buildings which are within. I created a field
> "batcount" in the neighbourhoods shapefile.
> However I do not have the geometries of the buildings in this shapefile
> and I do not know how to procede in order to find the area of the buildings
> in each neighbourhood.
>
> Here is the code I used to join the two shapefiles in order to count the
> number of buildings in each neighbourhood.
>
> bat_count_iris <- st_join(bat, neighbourhoods) %>%
> st_set_geometry(NULL) %>%
> group_by(CODE_NEIGHBOURHOOD) %>%
> tally() %>%
> select(CODE_NEIGHBOURHOOD, batcount = n)
> bat_in_iris <- left_join(CODE_NEIGHBOURHOOD, bat_count_iris, by =
> 'CODE_NEIGHBOURHOOD') %>%
> mutate(batcount = ifelse(is.na(batcount), 0, batcount))
>
> Do you have any clue of how I can compute the areas of the buildings in
> each neighbourhood?
>
> Have a great day,
> Basile.
>
> Sent with [ProtonMail](https://protonmail.com) Secure Email.
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] self-intersecting polygons

[Dexter Locke]

Dear list,

I have a set of 15 and 30 minute drive-time polygons created
from mapboxapi::mb_isochrone. Ultimately I want to erase the smaller
polygons from the larger so that the polygons are non-overlapping,
concentric rings.

The whole dataset includes ~1600 polygons for each drive-time spread out
over the 50 US states plus Washington DC and Puerto Rico. 13 of them have
apparent geometry errors, but the example below includes just 4 polygons
with one self intersection in Hawaii.

script (also below) and data available here:
http://dexterlocke.com/wp-content/uploads/2021/11/reprex_files.zip (I
couldn't figure out how to get the script to load data into R from a URL)

library(sf)
library(tidyverse)

# read in data
(small <- st_read(paste0(getwd(), '/data/reprex_files/small_polys.shp')))
(large <- st_read(paste0(getwd(), '/data/reprex_files/large_polys.shp')))

# look on basemap (optional)
mapview::mapview(small, color = 'black') + mapview::mapview(large, color =
'red')

# create function to erase holes to make 'donuts'
# https://github.com/r-spatial/sf/issues/346
st_erase <- function(x, y) st_difference(x, st_union(st_combine(y)))

# apply the erasing function
donut <- large %>% st_erase(., small) # fails

# 1st option explored, turn off spherical geometry as suggested here:
# https://github.com/r-spatial/sf/issues/1762

# here
https://gis.stackexchange.com/questions/404385/r-sf-some-edges-are-crossing-in-a-multipolygon-how-to-make-it-valid-when-using/404454

# and here
#
https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data/68481205#68481205

sf_use_s2(FALSE)

donut <- large %>% st_erase(., small) # results in errors..

# are the polygons valid? Yes
small %>% st_is_valid()
large %>% st_is_valid()

# 2nd option explored, more intensive geometry cleaning, as suggested here:
#
https://stackoverflow.com/questions/68478179/how-to-resolve-spherical-geometry-failures-when-joining-spatial-data/68481205#68481205
and elsewhere
small$geometry <- small$geometry %>%
  s2::s2_rebuild() %>%
  sf::st_as_sfc()

large$geometry <- large$geometry %>%
  s2::s2_rebuild() %>%
  sf::st_as_sfc()

# re-attempt
donut <- large %>% st_erase(., small) # results in errors

# 3rd option explored, zero-distance buffer to help reset geometries
small <- small %>% st_buffer(., 0)
large <- large %>% st_buffer(., 0)

# re-attempt
donut <- large %>% st_erase(., small) # more self-intersection errors

# 4th option explored, reproject before buffering
# suggested here:
https://gis.stackexchange.com/questions/404385/r-sf-some-edges-are-crossing-in-a-multipolygon-how-to-make-it-valid-when-using/404388#404388
small <- small %>% st_transform(crs = 5070) %>% st_buffer(0) %>%
st_transform(crs = 4269)
large <- large %>% st_transform(crs = 5070) %>% st_buffer(0) %>%
st_transform(crs = 4269)

# re-attempt
donut <- large %>% st_erase(., small) # more self-intersection errors

# 5th option
small_s2 <- small %>%
  s2::s2_rebuild(.,
 options = s2::s2_options(edge_type = "undirected",
 duplicate_edges = FALSE,
 split_crossing_edges = TRUE,
 validate = TRUE)) %>%
  st_as_sf()

large_s2 <- large %>%
  s2::s2_rebuild(.,
 options = s2::s2_options(edge_type = "undirected",
 duplicate_edges = FALSE,
 split_crossing_edges = TRUE,
 validate = TRUE)) %>%
  st_as_sf()

# re-attempt
donut <- large_s2 %>% st_erase(., small_s2) # more self-intersection errors


# 6th option explored using rgeos and sp, as suggested here:
#
https://gis.stackexchange.com/questions/163445/getting-topologyexception-input-geom-1-is-invalid-which-is-due-to-self-intersec/163480#163480

# same problems..

# so where is the double-vertex?
large %>%
  filter(id == 1017) %>% # based on zooming and panning around via mapview
  st_cast('POINT') %>%
  mutate(geom_text = as.character(geometry)) %>%
  group_by(geom_text) %>%
  count() %>%
  arrange(desc(n))

# as indicated from some errors and from query directly above
bad_point <- tibble(row = 1, x = -155.056700377001,  y = 19.5963974859705)
%>%
  st_as_sf(., coords = c('x', 'y')) %>%
  st_set_crs(4269)

# take a look
large %>% mapview::mapview() + mapview::mapview(bad_point, color = 'red')

# Ok, so there are two vertices there.. how can I remove the duplicates and
get on with the rest of the analyses?

Thank you, Dexter

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] crop or mask large raster on inport

[Dexter Locke]

Thanks all.

Best, Dexter


On Wed, Nov 25, 2020 at 4:13 AM Roger Bivand  wrote:

> On Tue, 24 Nov 2020, Dexter Locke wrote:
>
> > Dear list,
> >
> > Is there a way to reduce the spatial extent of an impractically-large
> > raster (~50Gb) on import? Can a crop or mask be applied as the large
> raster
> > is read into R to create a more manageable subset?
>
> rgdal::GDAL.open() simply opens the data source but does not read it,
> returning an object with a handle to the open file - close with
> rgdal::GDAL.close(). rgdal::getRasterData() has a complete set of
> arguments for the kinds of things you mention, see:
> https://rgdal.r-forge.r-project.org/reference/GDALRasterBand-class.html
>
> The same arguments are used by rgdal::readGDAL(): offset, region.dim,
> output.dim, band. These all permit rectangular cropping, and output.dim=
> can decimate/resample/aggregate input cells to coarser output cells. These
> are the same mechanisms that GDAL utility programs use. They are also used
> by raster internally.
>
> The example in rgdal::readGDAL() includes:
>
> SP27GTIF <- readGDAL(system.file("pictures/SP27GTIF.TIF",
>   package = "rgdal")[1], output.dim=c(100,100))
>
> doing resampling on-the-fly. Lower down,
>
> fn <- system.file("pictures/erdas_spnad83.tif", package = "rgdal")[1]
> erdas_spnad83 <- readGDAL(fn, offset=c(50, 100), region.dim=c(400, 400),
>   output.dim=c(100,100))
>
> uses all of offset=, region.dim= and output.dim=. The example files are
> tiny, but the same GDAL code applies to larger files. Note that the axis
> ordering of the argument vectors is often counter-intuitive and needs some
> patience.
>
> The GDAL utilities can be run file-to-file using the gdalUtils package,
> in addition to being run from the command line.
>
> Otherwise, see the stars package and its proxy=TRUE read functionality,
> which leaves the data on disk until required in possibly subsetted and
> resampled form.
>
> Hope this helps,
>
> Roger
>
> >
> > I've looked at raster and rgdal helpfiles, but have not found this
> > functionality. This would be akin to the sf::st_read() 's new "query"
> > argument that takes SQL statements for filtering records while reading
> data
> > into R.
> >
> > The ultimate use case is to summarize a large categorical raster into
> many
> > thousands of polygons, and I thought reading small subsets in at a time
> > could be looped through or parallelized for summarizing smaller pieces of
> > the data. If there are other approaches to zonal functions (like
> > equivalents of ArcMap's Tabulate Area) that could also work.
> >
> > Thanks for your consideration, Dexter
> >
> >   [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] crop or mask large raster on inport

[Dexter Locke]

Dear list,

Is there a way to reduce the spatial extent of an impractically-large
raster (~50Gb) on import? Can a crop or mask be applied as the large raster
is read into R to create a more manageable subset?

I've looked at raster and rgdal helpfiles, but have not found this
functionality. This would be akin to the sf::st_read() 's new "query"
argument that takes SQL statements for filtering records while reading data
into R.

The ultimate use case is to summarize a large categorical raster into many
thousands of polygons, and I thought reading small subsets in at a time
could be looped through or parallelized for summarizing smaller pieces of
the data. If there are other approaches to zonal functions (like
equivalents of ArcMap's Tabulate Area) that could also work.

Thanks for your consideration, Dexter

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] credible interval for empirical Bayesian estimates of rates

[Dexter Locke]

Wonderful thank you.

GeoDA and the spdep::EBest are providing the same estimates for my
datasets. Thanks for the reprex example with the GUI via ">". That is a
helpful convention I'll use in the future..

Ok great, there is not yet a CI implementation. If I can manage one I'd be
delighted to share and contribute it.

Julie Silge: https://juliasilge.com/blog/bayesian-blues/

and David Robinson
http://varianceexplained.org/r/credible_intervals_baseball/

have some suggestions that are similar to each other's (no surprise there,
they work together a lot).

Thanks for the suggestions about alternative model specifications and
including covariates.

This has been very helpful, I appreciate this list immensely,
Dexter



On Fri, Apr 24, 2020 at 3:05 PM Roger Bivand  wrote:

> On Fri, 24 Apr 2020, Dexter Locke wrote:
>
> > Thanks Roger and list.
> >
> > I didn't think a repex was needed because a question was: why does
> > spdep::EBest(counts, population, family = 'binomial') give the same
> > results at GeoDa's, while EBest(.. binomial) is "binomial" while GeoDa
> > calls that "Poisson-Gamma". GeoDa can't give use a repex (GUI) and think
> > this is a question about terminology. The same results were achieved with
> > the packages while naming the model differently - why?
>
> Reproducible example:
>
> auckland <- st_read(system.file("shapes/auckland.shp",
>package="spData")[1], quiet=TRUE)
> res <- EBest(auckland$M77_85, 9*auckland$Und5_81)
> res0 <- EBest(auckland$M77_85, 9*auckland$Und5_81, family="binomial")
>
> > summary(res$estmm)
>  Min.  1st Qu.   Median Mean  3rd Qu. Max.
> 0.001487 0.002237 0.002549 0.002648 0.002968 0.004531
> > summary(res0$estmm)
>  Min.  1st Qu.   Median Mean  3rd Qu. Max.
> 0.001484 0.002235 0.002549 0.002648 0.002968 0.004536
>
> After calculating Und5_81 * 9 as a new column, running GeoDa -> Map ->
> Rates-Caluculated Map -> Empirical Bayes and exporting a shapefile, in R I
> see:
>
> > summary(auck$R_EBS)
>  Min.  1st Qu.   Median Mean  3rd Qu. Max.
> 0.001487 0.002237 0.002549 0.002648 0.002968 0.004531
>
> which is the same as the Poisson, and:
>
> > all.equal(res$estmm, auck$R_EBS)
> [1] TRUE
>
> GeoDa is providing the same EB Poisson as EBest, isn't it? If yours
> differ, are both implementations seeing the same data?
>
> >
> > Yes ?spdep::EBest directed me to the literature I'm struggling to access.
> > And Yes, I've been looking at the raw code and understand how the estmm
> is
> > generated.
> >
> > I've been using the epitools::pois.exact() and spdep::EBest. I can
> compare
> > the point estimates from pois.exact to those provided by EBest, but I'd
> > like to graph side by side their credible / confidence intervals.
> >
> > Its this last part on the credible intervals I'm interested in. How to
> get
> > credible intervals around estmm? This is my main question.
>
> EBest() was written to implement Bailey & Gatrell's textbook, which did
> not provide CI, and just used the Marshall Auckland dataset. If you'd like
> to implement them, I'd welcome a contribution.
>
> >
> > ASDAR is a reference I'm using all the time. Thanks for that gem.
> >
> > DCluster::empbaysmooth also does not provide a credible interval, either.
> >
>
> As you see from ch. 10, CI are described for the epitools implementation.
> My feeling is that the literature moved away from simple EB rates towards
> IID RE and spatially structured RE, with relevant covariates, say like the
> classic Scottish Lip Cancer dataset, and currently CARBayes is a solid
> package among many others. Simply using base population becomes too
> unsatisfactory. PHE uses funnel plots which do have CI of a kind, to draw
> attention from small base populations:
> https://nhsrcommunity.com/blog/introduction-to-funnel-plots/ which can be
> used to adjust class intervals for mapping.
>
> Roger
>
> > -Dexter
> > http://dexterlocke.com/
> >
> >
> >
> > On Fri, Apr 24, 2020 at 10:23 AM Roger Bivand 
> wrote:
> >
> >> On Fri, 24 Apr 2020, Dexter Locke wrote:
> >>
> >>> Dear esteemed list,
> >>>
> >>> I'm using spdep::EBest with family = 'binomial' for counts of events
> >> within
> >>> polygons that have an 'at risk' population. The resultant "estmm" is
> >>> 'shrunk' compared to the raw rate (both given by EBest and calculated
> "by
> >>> hand" rate. All good there.
> >>>
> >>> Using GeoDa version 1.14.0 24 Augu

Re: [R-sig-Geo] credible interval for empirical Bayesian estimates of rates

[Dexter Locke]

Thanks Roger and list.

I didn't think a repex was needed because a question was: why does
spdep::EBest(counts, population, family = 'binomial') give the same
results at GeoDa's, while EBest(.. binomial) is "binomial" while GeoDa
calls that "Poisson-Gamma". GeoDa can't give use a repex (GUI) and think
this is a question about terminology. The same results were achieved with
the packages while naming the model differently - why?

Yes ?spdep::EBest directed me to the literature I'm struggling to access.
And Yes, I've been looking at the raw code and understand how the estmm is
generated.

I've been using the epitools::pois.exact() and spdep::EBest. I can compare
the point estimates from pois.exact to those provided by EBest, but I'd
like to graph side by side their credible / confidence intervals.

Its this last part on the credible intervals I'm interested in. How to get
credible intervals around estmm? This is my main question.

ASDAR is a reference I'm using all the time. Thanks for that gem.

DCluster::empbaysmooth also does not provide a credible interval, either.

-Dexter
http://dexterlocke.com/



On Fri, Apr 24, 2020 at 10:23 AM Roger Bivand  wrote:

> On Fri, 24 Apr 2020, Dexter Locke wrote:
>
> > Dear esteemed list,
> >
> > I'm using spdep::EBest with family = 'binomial' for counts of events
> within
> > polygons that have an 'at risk' population. The resultant "estmm" is
> > 'shrunk' compared to the raw rate (both given by EBest and calculated "by
> > hand" rate. All good there.
> >
> > Using GeoDa version 1.14.0 24 August 2019 produces identical results for
> > its Empirical Bayesian rate. This was confirmed by plotting the EBest
> > output against GeoDa's rate and finding a perfect correlation along the 1
> > to 1 line. All good there.
>
> Please provide a reproducible example, as this may help with answers.
>
> >
> > Two questions:
> > 1. How can credible intervals around these smoothed rate estimates be
> > calculated?
> > 2. The spdep documentation calls this a binomial family, but the
> identical
> > results are obtained from GeoDa calls this "Poisson-Gamma" model here:
> > https://geodacenter.github.io/workbook/3b_rates/lab3b.html#fnref11 , so
> > what is actually being calculated? This question may help me answer the
> > first question..
>
> No, the default family is "poisson", with "binomial" available for
> non-rare conditions following Martuzzi, implemented by Olaf
> Berke, ?spdep::EBest.
>
> The code in spdep is easily accessible, so can be read directly. Please
> also compare with code for the EB Moran test, and with analogous code in
> the DCluster package, empbaysmooth(). Cf. ASDAR 2nd ed., ch. 10, section
> 10.2, pp. 322-328. The epitools::pois.exact() function is used for CIs.
> For code and data see https://asdar-book.org/bundles2ed/dismap_bundle.zip.
>
> >
> > Possibly the answers are addressed in the literature cited which I cannot
> > access right now at home without institutional library access.
> >
>
> Most institutions do have proxy or VPN access, but the code will be as
> useful. In PySAL, the code would also guide you, but even though GeoDa is
> open source, the C++ is fairly dense.
>
> Hope this helps,
>
> Roger
>
> > Thanks for your consideration,
> > Dexter
> > http://dexterlocke.com/
> >
> >   [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] credible interval for empirical Bayesian estimates of rates

[Dexter Locke]

Dear esteemed list,

I'm using spdep::EBest with family = 'binomial' for counts of events within
polygons that have an 'at risk' population. The resultant "estmm" is
'shrunk' compared to the raw rate (both given by EBest and calculated "by
hand" rate. All good there.

Using GeoDa version 1.14.0 24 August 2019 produces identical results for
its Empirical Bayesian rate. This was confirmed by plotting the EBest
output against GeoDa's rate and finding a perfect correlation along the 1
to 1 line. All good there.

Two questions:
1. How can credible intervals around these smoothed rate estimates be
calculated?
2. The spdep documentation calls this a binomial family, but the identical
results are obtained from GeoDa calls this "Poisson-Gamma" model here:
https://geodacenter.github.io/workbook/3b_rates/lab3b.html#fnref11 , so
what is actually being calculated? This question may help me answer the
first question..

Possibly the answers are addressed in the literature cited which I cannot
access right now at home without institutional library access.

Thanks for your consideration,
Dexter
http://dexterlocke.com/

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Parcel scale or aggregation

[Dexter Locke]

Yes, those functions should work provided your data are shaped correctly
(they sound like they are).

I imagine the results will be sensitive the chosen spatial weights matrix.
That choice indirectly gets to your point about aggregation. Some types of
weights may include spatial neighbors that might assimilate neighborhoods
or blocks.

Good luck!

-Dexter
http://dexterlocke.com/


On Tue, Apr 21, 2020 at 11:23 AM derek  wrote:

> Dexter.
>
> Your point that it is more of a research design question is well taken,
> but I figured if anyone would know it would be a group of folks using R for
> spatial analysis. The research question does pertain to the parcels and the
> explanatory variables will be computed at that scale. I just didn't know if
> that matter per say for functions like sarlm() or or gmerrorsar(). It
> sounds like not so much as they will be treat the same as rows in a data
> frame regardless of what scale they are at.
>
> Thanks,
>
> Derek
> On 4/21/2020 9:59 AM, Dexter Locke wrote:
>
> Hello.
>
> This is more of a research design question than a spatial
> analysis question.
>
> If you research question pertains to parcels and you have parcel data,
> then why aggregate?
>
> Neighborhood-level attributes can be important. They can be included by
> attributing parcels with their neighborhood characteristics, with random
> effects at the neighborhood scale, among other techniques.
>
> The chosen method depends on the research questions.
>
> -Dexter
>
>
>
> On Tue, Apr 21, 2020 at 9:34 AM John Morgan  wrote:
>
>> Hello, I am working on preparing data for to run in a spatial
>> autoregression model (probably SEM). You are one of the most well
>> grounded people I am aware of in this type of methodology. I am hoping
>> you can
>> help me with a simple question.
>>
>> My question is this: we have the data at the scale of the parcel (or
>> household) and there are a couple of hounded thousand records across the
>> size of a large metropolitan area. When feeding the data into the
>> model, is there
>> a reason/requirement to aggregate our data variable to some boundary scale
>> such as a city block? Or is it ok to keep it at the parcel scale? We are
>> interested in analyzing characteristics at e.g. household level similar to
>> a hedonic model.
>>
>> Thanks for any feedback.
>>
>> [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> --
> John Derek Morgan, Ph.D., GISP
> Assistant Professor of GIS
> Earth & Environmental Sciences
> University of West 
> Floridahttps://uwf.edu/go/gis/https://pages.uwf.edu/jmorgan3
> >>>Chat directly with me via Google Chat our Hangout<<<
>
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Parcel scale or aggregation

[Dexter Locke]

Hello.

This is more of a research design question than a spatial analysis question.

If you research question pertains to parcels and you have parcel data, then
why aggregate?

Neighborhood-level attributes can be important. They can be included by
attributing parcels with their neighborhood characteristics, with random
effects at the neighborhood scale, among other techniques.

The chosen method depends on the research questions.

-Dexter



On Tue, Apr 21, 2020 at 9:34 AM John Morgan  wrote:

> Hello, I am working on preparing data for to run in a spatial
> autoregression model (probably SEM). You are one of the most well
> grounded people I am aware of in this type of methodology. I am hoping you
> can
> help me with a simple question.
>
> My question is this: we have the data at the scale of the parcel (or
> household) and there are a couple of hounded thousand records across the
> size of a large metropolitan area. When feeding the data into the
> model, is there
> a reason/requirement to aggregate our data variable to some boundary scale
> such as a city block? Or is it ok to keep it at the parcel scale? We are
> interested in analyzing characteristics at e.g. household level similar to
> a hedonic model.
>
> Thanks for any feedback.
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Merging Shape files

[Dexter Locke]

Within st_join() see the options like "left = TRUE".

You may also want to look at st_intersects()

Here are other helpful reference websites:
https://r-spatial.github.io/sf/reference/geos_binary_pred.html

https://cran.r-project.org/web/packages/sf/vignettes/sf3.html

-Dexter


On Tue, Mar 17, 2020 at 10:29 AM El Mechry El Koudouss <
eelkoudo...@fordham.edu> wrote:

> Dear readers,
> I downloaded a shapefile with data on crime incidents in Washington DC in
> 2016 (available here:
> https://opendata.dc.gov/datasets/crime-incidents-in-2016) and another
> shapefile with data on parks and recreation areas, also in Washington DC
> (available here:
> https://opendata.dc.gov/datasets/parks-and-recreation-areas).
> I was able to read both shapefiles using st_read() from the sf package. My
> question is, how can I merge the two data sets? I tried st_join() but it
> simply resulted in empty columns from the parks dataset getting added to
> the crime data. Any guidance would be much appreciated.
> library(sf)
> parks <- st_read("Parks_and_Recreation_Areas.shp")
> crime <- st_read("Crime_Incidents_in_2016.shp")
> df <- st_join(crime, parks)
>
> --
> El Mechry, El Koudouss (Meshry)
> Graduate Research Assistant
> Center for International Policy Studies
> Fordham University
> Department of Economics
> Website: www.meshry.com
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] poly2nb neighbour itself should be considered a neighbour

[Dexter Locke]

Dear Robert,

It sounds like what you are looking for is typically called a second order 
neighbor. Higher order neighbors can also included in a weights matrix such as 
your neighbors’, neighbors’, neighbor which is a third-order neighbor. I think 
you are seeking second-order neighbors. 

See the spdep vignettes, and the section 5 Higher-Order Neighbors subsection 
here: https://cran.r-project.org/web/packages/spdep/vignettes/nb.pdf in 
particular. The spdep::nblag might be what you need, but without additional 
information it is hard to know. 

Good luck,
Dexter
dexterlocke.com



> On Nov 3, 2019, at 2:12 PM, Robert R  wrote:
> 
> Dear All,
> 
> I would like to know if the function "poly2nb" ("spdep" package.) let me 
> create a neighborhood of itself, i.e., not only its queen neighbors 
> (queen=TRUE), but a neighbour itself should also be considered a neighbour.
> 
> I am looking to create a queen weight neighborhood matrix afterwards using 
> "nb2listw".
> 
> Any help would help me a lot.
> 
> Many thanks
> 
>[[alternative HTML version deleted]]
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] CSV with Geometry Column to SF object

[Dexter Locke]

For what it's worth, Edzer's reproducible example worked just fine for me
and did not throw the error described.

I'm using R version 3.6.0 (2019-04-26) -- "Planting of a Tree" on a Mac.

-Dexter



On Tue, Oct 15, 2019 at 3:30 PM Edzer Pebesma 
wrote:

> Without sharing a (minimal) reproducible example, it is unlikely that
> someone else can help you find out whether this points to a problem in
> your data, or in the software.
>
> On 10/15/19 2:31 PM, argunaw . wrote:
> >   When I run this, I get the following error:
> >
> > *Error in CPL_sfc_from_wkt(x) : OGR error*
> >
> > When I run it on my data, I get the same error.
> >
> > On Tue, Oct 15, 2019 at 2:16 PM Edzer Pebesma
> > mailto:edzer.pebe...@uni-muenster.de>>
> > wrote:
> >
> > You may try something along these lines:
> >
> > # read data.frame with read.csv; here, we create an example by hand:
> > df = data.frame(
> > a = 1:3, b = 3:1, geom = c("LINESTRING(0 0, 1 1)",
> > "LINESTRING(1 1,2
> > 2)", "LINESTRING(5 5,6 6)")
> > )
> >
> > library(sf)
> > sf = st_sf(df, geom = st_as_sfc(df$geom))
> > sf
> >
> >
> > On 10/15/19 1:58 PM, argunaw . wrote:
> > > I'm not sure how it was exported from postgis- the person who gave
> > me the
> > > file wasn't the one who downloaded it unfortunately.
> > >
> > > The file is a line file of roads. The files main columns are road
> ID
> > > numbers (type integer) and the geometry column (type geometry,
> > long strong
> > > of letters and numbers). Only the IT admins where I am have the
> > postgis
> > > load/import tools in the pgadmin/sql interface. The rest of us can
> > download
> > > from sql and create new tables from other sql databases, but not
> > create a
> > > new table from a csv file.
> > >
> > > On Tue, Oct 15, 2019 at 1:11 PM Alex Mandel
> > mailto:tech_...@wildintellect.com>>
> > > wrote:
> > >
> > >> On 10/15/19 8:43 AM, argunaw . wrote:
> > >>> Hello Everyone,
> > >>>
> > >>> I have a csv file with a postgis "geometry" column. I've been
> > trying to
> > >>> import it in to R as a SF file, with the goal of exporting it to
> a
> > >> postgis
> > >>> database, but to no avail. I've used the following methods:
> > >>>
> > >>> 1. file <- st_read("name.csv", stringsAsFactors=F,
> > geometry_column=geom)
> > >>>
> > >>> 2. file <- fread("name.csv", headers=True)
> > >>> file <- st_as_sf(file)
> > >>>
> > >>> How can I import a csv with a postgis "geometry" column in to R
> as a
> > >>> spatial/SF object?
> > >>>
> > >>>   [[alternative HTML version deleted]]
> > >>>
> > >>> ___
> > >>> R-sig-Geo mailing list
> > >>> R-sig-Geo@r-project.org 
> > >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >>>
> > >>
> > >> Can you paste an example somewhere, is it binary data or some
> kind of
> > >> plain text column? Do you know how it was exported from Postgis?
> > >>
> > >> If it's a dump from a postgis database did you try loading the
> table
> > >> directly to postgis with it's own load/import, or sql tools?
> > >>
> > >> Thanks,
> > >> Alex
> > >>
> > >
> > >   [[alternative HTML version deleted]]
> > >
> > > ___
> > > R-sig-Geo mailing list
> > > R-sig-Geo@r-project.org 
> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> > >
> >
> > --
> > Edzer Pebesma
> > Institute for Geoinformatics
> > Heisenbergstrasse 2, 48151 Muenster, Germany
> > Phone: +49 251 8333081
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org 
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] [FORGED] Re: Aggregating points based on distance

[Dexter Locke]

FWIW: I agree with Rolfs nomination.

+1

-Dexter
http://dexterlocke.com



On Thu, Mar 14, 2019 at 5:30 PM Rolf Turner  wrote:

>
> On 14/03/19 7:33 AM, Barry Rowlingson wrote:
>
> 
>
> > The problem with "modern" syntax is that it's subject to rapid change
> > and often slower than using base R, which has had years to stabilise
> > and optimise.
>
> 
>
> Fortune nomination.
>
> cheers,
>
> Rolf
>
> --
> Honorary Research Fellow
> Department of Statistics
> University of Auckland
> Phone: +64-9-373-7599 ext. 88276
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Fitting GWR for Panel Data in R

[Dexter Locke]

See

Fotheringham, A. S., Crespo, R., & Yao, J. (2015). Geographical and
Temporal Weighted Regression (GTWR). *Geographical Analysis*, 1–22.
https://doi.org/10./gean.12071

Gollini, I., Lu, B., Charlton, M., Brunsdon, C., & Harris, P. (2015).
GWmodel: An R Package for Exploring Spatial Heterogeneity Using
Geographically Weighted Models. *Journal of Statistical Software*, *63*(17),
1–50. https://doi.org/10.1080/10095020.2014.917453
HTH, Dexter


On Wed, Jun 27, 2018 at 7:28 AM Arnold Salvacion via R-sig-Geo <
r-sig-geo@r-project.org> wrote:

> Good day!
>
> Does anyone here have experience or code for fitting geographically
> weighted regression (GWR) for panel data. For example, I need to determine
> the effect of independent variable X1, X2, and X3 to Y which are all
> collected for T1 to Tn.  I have experience in fitting GWR for
> cross-sectional data but have no experience on doing such on panel data.
>
> Any suggestion is highly appreciated.
>
> Thanks in advance!
>
> Arnold Salvacion
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] nb object in spdep from categorical attribute

[Dexter Locke]

Thanks so much.

The nb2blocknb(), aggregate.nb(), and union.nb() are fantastic functions!
Between these three functions I that I had not seen before I can accomplish
what I am seeking. It looks like the new spdep vignette has been greatly
expanded, which is also fantastic.

Thanks again,
Dexter


On Wed, Jun 20, 2018 at 8:24 AM Roger Bivand  wrote:

> On Tue, 19 Jun 2018, Dexter Locke wrote:
>
> > Dear list,
> >
> > Does anyone know how to make nb objects (to eventually make spatial
> > weights) with spdep using an attribute in a shapefile's data frame?
> >
> > Consider, for example, all property parcels on the same street segment
> > should be defined as neighbors. In the parcel polygons I have the
> > associated street segment name.
> >
> > KNN forces sparsely settled areas to become neighbors, which is
> undesirable
> > for the subsequent analyses.
> >
> > Distances between parcels vary, so distance-based neighbors are also
> > undesirable in this case.
> >
> > Contiguity-based methods do not connect parcels on opposite sites of the
> > same street segment.
> >
> > Is there a way to use an attribute to define what constitutes a neighbor?
>
> Yes, spdep::nb2blocknb() is a specific case where say observations without
> positional data are known to belong to an aggregate entity for which we
> have positional data. The output is an nb object with diagonal blocks (if
> the observations are sorted by block ID, not in the case below), as this
> example shows. You'll need the development version to run it, as released
> spdep::aggregate.nb() failed with empty aggregate sets - the nc median
> poish grid has such empty sets:
>
> # devtools::install.github("r-spatial/spdep")
> library(spdep)
> data(nc.sids)
> ncCC89.nb[sapply(ncCC89.nb, length) == 0L] <- 0L
> image(as(nb2listw(ncCC89.nb, zero.policy=TRUE), "CsparseMatrix"))
> anb <- aggregate(ncCC89.nb, interaction(nc.sids$M.id, nc.sids$L.id))
> image(as(nb2listw(anb, style="B"), "CsparseMatrix"))
> # when the neighbour object is not empty
> bnb <- nb2blocknb(anb, as.character(interaction(nc.sids$M.id,
>nc.sids$L.id)))
> image(as(nb2listw(bnb, style="B"), "CsparseMatrix"))
> # and when the neighbour object is empty
> anb1 <- lapply(anb, function(x) x <- 0L)
> attributes(anb1) <- attributes(anb)
> bnb1 <- nb2blocknb(anb1, as.character(interaction(nc.sids$M.id,
>nc.sids$L.id)))
> image(as(nb2listw(bnb1, style="B"), "CsparseMatrix"))
>
> This isn't your case, but crafting an nb object by hand (adding
> attributes after creating the list) then using spdep::union.nb or similar
> to combine with a positional nb (such as poly2nb) may work. For parcels
> across a stream or street, maybe look at the snap= argument if the streets
> are narrower than frontages. Could you provide a small reproducible
> example?
>
> Hope this clarifies,
>
> Roger
>
> >
> > All ideas welcome,
> > Dexter
> >
> >   [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
> http://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] nb object in spdep from categorical attribute

[Dexter Locke]

Dear list,

Does anyone know how to make nb objects (to eventually make spatial
weights) with spdep using an attribute in a shapefile's data frame?

Consider, for example, all property parcels on the same street segment
should be defined as neighbors. In the parcel polygons I have the
associated street segment name.

KNN forces sparsely settled areas to become neighbors, which is undesirable
for the subsequent analyses.

Distances between parcels vary, so distance-based neighbors are also
undesirable in this case.

Contiguity-based methods do not connect parcels on opposite sites of the
same street segment.

Is there a way to use an attribute to define what constitutes a neighbor?

All ideas welcome,
Dexter

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Generate the local coefficients in spgwr

[Dexter Locke]

The R script here might be helpful: DOI: 

https://doi.org/10.5061/dryad.3vh79

-Dexter 

> On May 4, 2018, at 8:09 PM, Danlin Yu  wrote:
> 
> Alysson:
> 
> I believe if you calibrate your GWR model using the gwr() routine, the 
> resulted gwr object shall have all the elements you need in the SDF element, 
> as detailed in the package's manual:
> 
> The SDF is a SpatialPointsDataFrame (may be gridded) or 
> SpatialPolygonsDataFrame object (see package "sp") with fit.points, weights, 
> GWR coefficient estimates, R-squared, and coefficient standard errors in its 
> "data" slot.
> 
> Once you extract this information, you can output them using generic 
> functions like write.table to output to txt files, csv files and the like.
> 
> Hope this helps.
> 
> Cheers,
> 
> Dr. Danlin Yu
> 
> 
>> On 5/4/2018 5:17 PM, Alysson Luiz Stege wrote:
>> Hello, my name is Alysson.
>> 
>> I am estimating a model by the geographically weighted regression method
>> and using the spgwr package. I would like to know how to generate the local
>> coefficients, the value of the statistic t or the level of significance,
>> the standard deviation and then export them into a txt file, for example.
>> 
>> Thank's
>> 
>> Alysson
>> 
>>[[alternative HTML version deleted]]
>> 
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 
> -- 
> ___
> Danlin Yu, Ph.D.
> Professor of GIS and Urban Geography
> Department of Earth & Environmental Studies
> Montclair State University
> Montclair, NJ, 07043
> Tel: 973-655-4313
> Fax: 973-655-4072
> Office: CELS 314
> Email: y...@mail.montclair.edu
> webpage: csam.montclair.edu/~yu
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] lm.morantest

[Dexter Locke]

Hi Ana,

Maybe you need to use the zero.policy = T option when creating your weights
matrix? Alternatively, consider a different style of weights that will
surely create connections.

Since your data are multi-level, consider also the ICC (Intraclass
Correlation Coefficient).

The Moran's I test of the residuals and ICC will help determine the degree
of spatial autocorrelation, and hierarchical correlation, respectively.

HTH,
Dexter



On Mon, Aug 22, 2016 at 11:13 AM, Anita Morales 
wrote:

> Dear all,
>
> I'm relative new to the list, being this my first question. I have a
> problem when trying to estimate a Moran test for residuals in the second
> level when there are missing areas. My dataset contains information on
> individual nested in cities, the shapefile contains information of the 342
> cities, while my dataset contains individuals that belong to 288 cities.
>
> I've used multilevel models to estimate area level variation and now I want
> to test the autocorrelation for the residuals in the second level (cities).
> I'm using "lm.morantest"; however, I'm getting an error about "areas with
> no neighbours" (I already checked that and I'm quite sure that all my areas
> are connected). I think the problem may be caused by the missing data
> (residuals) for some areas. Since the listw object created, based on the
> shape file, contains neighbourhood weights for 342 cities, while from the
> estimated multilevel model, I've got residuals data for only 288 cities.
>
> Could you please advise me as to how to test for autocorrelation in the
> residuals when there are missing data?
> --
> Ana
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] spatial anova

[Dexter Locke]

Normally I would agree with Steve, but this is a topic I'm interested in
and have searched extensively. I have not yet found a useable
implementation. If there is some R code that can do this, I would love to
hear about it.

In April Babak Naimi  posted on this list a very similar
question, but I did not see any replies.

For the underlying maths, see these papers:

Griffith, D. A. (1978). A spatially adjusted ANOVA model. *Geographical
Analysis*, *X*(3), 298–301. doi:10.1016/0166-0462(92)90034-X

Griffith, D. A. (1992). A spatially adjusted N-way ANOVA model. *Regional
Science and Urban Economics*, *22*(3), 347–369.
doi:10.1016/0166-0462(92)90034-X

These are also related:

Long, J. a., Nelson, T. a., & Wulder, M. a. (2010). Local indicators for
categorical data: Impacts of scaling decisions. *Canadian Geographer*, *54*(1),
15–28. doi:10./j.1541-0064.2009.00300.x

Pietrzak, M. B., Wilk, J., Kossowski, T., & Bivand, R. (2014). The
application of local indicators for categorical data (LICD) in the spatial
analysis of economic development, (14). Retrieved from
http://ideas.repec.org/p/pes/wpaper/2014no14.html

The old SpaceStat used to have a spatial anova option (white paper by
Anselin):

http://siteresources.worldbank.org/INTPGI/Resources/342674-1092157888460/Anselin.spacestatTutorial.pdf

SpaceStat 4 appears to be able to perform a spatial anova:

https://www.biomedware.com/files/SpaceStat_4.0_Documentation.pdf


Just to be sure I hadn't missed anything I double checked and found this:

library(spdep); help(joincount.multi) and joincount.test(). So it looks
like Roger Bivand provided the needed code once again! Thanks so much.

HTH,
Dexter




On Tue, Jun 28, 2016 at 12:20 PM, Steve Friedman 
wrote:

> I would start by reading google and R posts regarding spatial anova. There
> is plenty of help among that literature.  Formulate your model, share it,
> and then you might have better questions.
>
> Steve
>
> On Tuesday, June 28, 2016, Csany22 csany  wrote:
>
> > Hi R-sig-geo people,
> >
> >
> > what would be the best way to run spatial anova in R? I have a dependent
> > (continuous), independent (categorical) variables and lat and long fo
> reach
> > of the cell(pixel) based observations. I need to run anova. as would like
> > to compare the means (Tukey) of the independent variable.
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org 
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
> --
> Steve Friedman
> Cell: 517-648-6290
> Web:  www.envisioningecosystemsdecisions.com
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] test spatial autocorrelation of level 2 residuals in lme()

[Dexter Locke]

Apologies for cross-posting, its unclear if this is a spatial or a mixed
models question. Maybe its both?

I am interested in extracting the residuals from level 2 of a mixed model
(created using lme()), building a spatial weights matrix, and then testing
for spatial autocorrelation using Moran's I.

While I have found methods of incorporating spatial effects into a mixed
model using corStruct, I am interested in first evaluating *if* that is an
appropriate model for a given dataset by examining the level 2 residuals'
spatial patterning - or lack thereof.

Which slot contains the residuals for level 2 and how are they ordered?

Best,
Dexter

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Recovering GWR Predicted Values from GGWR function

[Dexter Locke]

Hi list,

Is there an in print citation to attribute this idea to The aim of spgwr
has always been to explore the weaknesses of GWR as a technique, not to
help people use it?

Or should I use something like (R. Bivand, personal communication)?

Thanks and happy R coding,
Dexter

On Wed, May 20, 2015 at 2:27 PM, Roger Bivand roger.biv...@nhh.no wrote:

 On Wed, 20 May 2015, Brown, Jason wrote:

  Hello,

 I'm estimating a Poisson GWR using the ggwr function.  I would like to
 recover the predicted values implied by the gwr model.  However, thus far
 I've only been able to recover the fitted values from the global model
 using $lm$fitted.values.  These obviously are the same values I get from a
 standard Poisson model. Is it possible to get the gwr predicted values
 instead?


 No, they are not provided. They would only be available if the data points
 and the fit points were identical, or if newdata was available (and was
 admitted by the function, which it is not). The aim of spgwr has always
 been to explore the weaknesses of GWR as a technique, not to help people
 use it, so no provision for GLM-GWR predictions is likely to be forthcoming
 here (there are no coefficient standard errors either).

 Hope this clarifies,

 Roger


 Thank you,
 Jason



 [[alternative HTML version deleted]]

 ___
 R-sig-Geo mailing list
 R-sig-Geo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-geo


 --
 Roger Bivand
 Department of Economics, Norwegian School of Economics,
 Helleveien 30, N-5045 Bergen, Norway.
 voice: +47 55 95 93 55; fax +47 55 95 91 00
 e-mail: roger.biv...@nhh.no


 ___
 R-sig-Geo mailing list
 R-sig-Geo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-geo


[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo