CSRの帰無仮設を棄却する検定手順(1)背景と分析地区(データクリーニング後)
-
背景
点過程データを分析する際、その最初のステップとして、CSR(Complete Spatial Randomness)という帰無仮設を置き、最近隣距離法若しくは法格法を用いて、その仮設を棄却する確率を算出する。*1*2ここでは、Rのspatstatというパッケージ*3を用いて、佐賀県と福岡県の意匠登録者の空間点過程データをBivand, Pebesma, Gómez-Rubio (2008)が説明するDiggle(2003) の例に従って、分析する。
-
分析地区
spパッケージのSpatialPointsオブジェクトに対して、spatstatがpppオブジェクトを利用するため、spatstatの機能を利用するようにspSagaFukuokaの地区を抽出した後、pppオブジェクトに変換しなければならない。
-
佐賀県と福岡県の全地区(pppSagaFukuoka)
> ## Bivand, Pebesma, Gómez-Rubio (2008)pp156-158 > library(maptools) > summary(spSagaFukuoka) Object of class SpatialPoints Coordinates: min max x 129.81595 131.11875 y 33.00453 33.95653 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] Number of points: 1255 > ## convert to scaled rectangle > spSagaFukuokaScaled <- elide(spSagaFukuoka, scale=TRUE, unitsq = FALSE) > summary(spSagaFukuokaScaled) Object of class SpatialPoints Coordinates: min max x 0 1.0000000 y 0 0.7307327 Is projected: NA proj4string : [NA] Number of points: 1255 > pppSagaFukuoka <- as( spSagaFukuokaScaled, "ppp") > summary(pppSagaFukuoka) Planar point pattern: 1255 points Average intensity 1255 points per square unit Coordinates are given to 8 decimal places Window: rectangle = [0, 1] x [0, 1] units Window area = 1 square unit
-
福岡市―筑紫野市(pppFukuokaChikushino)
> ## done in earlier session > llCRS <- CRS("+proj=longlat +ellps=WGS84") > ## coordinates > Sr1 = Polygon(cbind( + c(130.191354,130.626400,130.626400,130.191354, 130.191354), + c(33.718652,33.718652,33.423697,33.423697,33.718652) + )) > Srs1 = Polygons(list(Sr1), "bbFukuokaChikushino") > spbbFukuokaChikushino = SpatialPolygons(list(Srs1), 1:1, proj4string = llCRS) > summary(spbbFukuokaChikushino) Object of class SpatialPolygons Coordinates: min max x 130.1914 130.62640 y 33.4237 33.71865 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] > spFukuokaChikushino <- spSagaFukuoka[spbbFukuokaChikushino,] > summary(spFukuokaChikushino) Object of class SpatialPoints Coordinates: min max x 130.19690 130.58931 y 33.42709 33.71665 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] Number of points: 601 >
> ## convert to scaled rectangle > plot(spFukuokaChikushino, pch=1) > spFukuokaChikushinoScaled <- elide(spFukuokaChikushino, scale=TRUE, unitsq = FALSE) > summary(spFukuokaChikushinoScaled) Object of class SpatialPoints Coordinates: min max x 0 1.0000000 y 0 0.7379017 Is projected: NA proj4string : [NA] Number of points: 601 > pppFukuokaChikushino <- as(spFukuokaChikushinoScaled, "ppp") > summary( pppFukuokaChikushino) Planar point pattern: 601 points Average intensity 814.4716 points per square unit Coordinates are given to 8 decimal places Window: rectangle = [0, 1] x [0, 0.7379017] units Window area = 0.737902 square units > png(file="pppFukuokaChikushino.png", width=600, height=480) > plot( pppFukuokaChikushino, pch=1) > dev.off()
-
> ########## ChuoHakata > Sr2 = Polygon(cbind( + c(130.354992,130.486107,130.486107,130.354992, 130.354992), + c(33.613987,33.613987,33.535556,33.535556,33.613987) + )) > Srs2 = Polygons(list(Sr2), "bbChuoHakata") > spbbChuoHakata = SpatialPolygons(list(Srs2), 1:1, proj4string = llCRS) > summary(spbbChuoHakata) Object of class SpatialPolygons Coordinates: min max x 130.35499 130.48611 y 33.53556 33.61399 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] > spChuoHakata <- spSagaFukuoka[spbbChuoHakata,] > summary(spChuoHakata) Object of class SpatialPoints Coordinates: min max x 130.35518 130.48578 y 33.53611 33.61334 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] Number of points: 348 >
> ## convert to scaled rectangle > spChuoHakataScaled <- elide(spChuoHakata, scale=TRUE, unitsq = FALSE) > summary(spChuoHakataScaled) Object of class SpatialPoints Coordinates: min max x 0 1.0000000 y 0 0.5913476 Is projected: NA proj4string : [NA] Number of points: 348 > pppChuoHakata <- as(spChuoHakataScaled, "ppp") > summary(pppChuoHakata) Planar point pattern: 348 points Average intensity 588.4863 points per square unit Coordinates are given to 8 decimal places Window: rectangle = [0, 1] x [0, 0.5913476] units Window area = 0.591348 square units
-
春日市(pppKasuga)
> ## Earlier Saga-Fukuoka, Fukuoka-Chikushino Cities, Chuo-Hakata Wards > ## subsets show high degrees of aggregation. But this is to be expected > ## from terrain and population. This subset is from an area of Kasuga City > ## Chikushino cities that has had a lot of planned development with > ## that are newer developments with large areas of single family dwellings. > ## 33.556279, 130.481984 > ## 33.516162, 130.409705 > latNE = 33.556279 > longNE = 130.481984 > latSW = 33.516162 > longSW = 130.409705 > Sr1 = Polygon(cbind( + c(longNE ,longNE ,longSW ,longSW , longNE ), + c(latNE,latSW ,latSW ,latNE,latNE) + )) > Srs1 = Polygons(list(Sr1), "bbKasuga") > spbbKasuga = SpatialPolygons(list(Srs1), 1:1, proj4string = llCRS) > summary(spbbKasuga) Object of class SpatialPolygons Coordinates: min max x 130.40971 130.48198 y 33.51616 33.55628 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] > spKasuga <- spSagaFukuoka[spbbKasuga,] > summary(spKasuga) Object of class SpatialPoints Coordinates: min max x 130.40993 130.4798 y 33.51632 33.5561 Is projected: FALSE proj4string : [+proj=longlat +ellps=WGS84] Number of points: 50 >
> ## convert to scaled rectangle > spKasugaScaled <- elide(spKasuga, scale=TRUE, unitsq = FALSE) > summary(spKasugaScaled) Object of class SpatialPoints Coordinates: min max x 0 1.0000000 y 0 0.5692616 Is projected: NA proj4string : [NA] Number of points: 50 > > pppKasuga <- as(spKasugaScaled, "ppp") > summary(pppKasuga) Planar point pattern: 50 points Average intensity 87.83308 points per square unit Coordinates are given to 8 decimal places Window: rectangle = [0, 1] x [0, 0.5692616] units Window area = 0.569262 square units
-