第三章

全体を通じて使用するパッケージ

library(sp)

3.1 密度

library(spatstat)
# 正規分布の場合
px <- rnorm(500, mean=0.5, sd=0.15)
py <- rnorm(500, mean=0.5, sd=0.15)
# 一様分布の場合(参考)
# px <- runif(500)
# py <- runif(500)
pz <- as.data.frame(rep(1, 500))
colnames(pz) <- c("pz")
# quadratcount用
pnt <- ppp(px, py, c(0,1), c(0,1))
plot(pnt, type="p")
plot(density(pnt), 0.1)
contour(density(pnt), add=T)
plot(pnt, type="p", add=T)
# overlay用
pnt_xy <- cbind(px, py)
pnt_sp <- SpatialPoints(data.frame(px, py))
pnt_spdf <- SpatialPointsDataFrame(pnt_sp, pz)
# plot(pnt_spdf)
quadratcount(pnt, nx=2, ny=2)
quadratcount(pnt, nx=2, ny=2)/(0.5^2)
plot(pnt, type="p")
plot(quadratcount(pnt, nx=2, ny=2), add=T, col="red")
quadratcount(pnt, nx=4, ny=4)
quadratcount(pnt, nx=4, ny=4)/(0.25^2)
plot(pnt, type="p")
plot(quadratcount(pnt, nx=4, ny=4), add=T, col="red")
poly1 <- cbind(
c(0, 0.75, 0.75, 0.5, 0.5, 0.25, 0.25, 0, 0), 
c(0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 0))
poly2 <- cbind(
c(0.25, 0.75, 0.75, 0.25, 0.25), 
c(0.5, 0.5, 0.75, 0.75, 0.5))
poly3 <- cbind(
c(0.75, 1, 1, 0.75, 0.75, 0.5, 0.5, 0.75, 0.75), 
c(0, 0, 0.75, 0.75, 0.5, 0.5, 0.25, 0.25, 0))
poly4 <- cbind(
c(0, 1, 1, 0, 0), 
c(0.75, 0.75, 1, 1, 0.75))
poly1_pl <- Polygons(list(Polygon(poly1)), "poly1")
poly2_pl <- Polygons(list(Polygon(poly2)), "poly2")
poly3_pl <- Polygons(list(Polygon(poly3)), "poly3")
poly4_pl <- Polygons(list(Polygon(poly4)), "poly4")
poly_sp <- SpatialPolygons(list(poly1_pl, poly2_pl, poly3_pl, poly4_pl))
poly_spdf <- SpatialPolygonsDataFrame(poly_sp, data.frame(c(1:4), 
row.names=c("poly1", "poly2", "poly3", "poly4")))
plot(poly_spdf)
plot(pnt, type="p", add=T)
overlay(pnt_spdf, poly_spdf, fn=colSums)
area.poly_spdf <- c(0.25^2*6, 0.25^2*2, 0.25^2*4, 0.25^2*4)
overlay(pnt_spdf, poly_spdf, fn=colSums)/area.poly_spdf

3.2 属性値の基本統計量と標準化

library(spdep)
library(maptools)
library(classInt)
# データの読み込み
lph <- read.table("lph2010.csv", sep=",", header=TRUE)
x <- lph$Easting
y <- lph$Northing
z <- as.data.frame(cbind(lph$JCODE, lph$lph2010))
colnames(z) <- c("JCODE", "lph2010")
lph_sp <- SpatialPoints(data.frame(x, y))
lph_spdf <- SpatialPointsDataFrame(lph_sp, z)
# 図3.3
yoko_spdf <- lph_spdf[lph_spdf$JCODE>14000 & lph_spdf$JCODE<14200,]
pal0 <- c("grey80","black")
q_lph <- classIntervals(round(yoko_spdf$lph2010/10000,1), n=4, style="fisher")
# plot(q_lph, pal=pal0)
q_lph_Col <- findColours(q_lph,pal0)
plot(yoko_spdf,col=q_lph_Col, pch=20, cex=1.8, axes=TRUE, cex.axis=1.5)
legend("topright",fill=attr(q_lph_Col,"palette"),
legend=names(attr(q_lph_Col,"table")), cex=1.3, bty="n")
# plot(lph_spdf,col=q_lph_Col, pch=20, cex=3, axes=TRUE, cex.axis=1.5,
# xlim=c(502600,503100), ylim=c(128000,128500))
# legend("bottomright",fill=attr(q_lph_Col,"palette"),
# legend=names(attr(q_lph_Col,"table")), cex=1.3, bg="white")
# spplotを使う方法
# spplot(lph_spdf, c("lph2010"))
# library(RColorBrewer)
# spplot(lph_spdf, c("lph2010"), col.regions=brewer.pal(5, "Greys"))
# 平均・分散・標準編纂
# 標本平均
mean(lph$lph2010)
# 標本分散
var(lph$lph2010) * (length(lph$lph2010)-1)/length(lph$lph2010)
# 標本標準偏差
sqrt(var(lph$lph2010) * (length(lph$lph2010)-1)/length(lph$lph2010))
# 不偏分散
var(lph$lph2010)
# 不偏標準偏差
sd(lph$lph2010)
# 標準化
lph$lph2010.scale <- scale(lph$lph2010)
summary(lph$lph2010.scale)
# plot(density(lph$lph2010.scale))
hist(lph$lph2010.scale, col="grey", xlim=c(-2,4), ylim=c(0,4000), breaks=48,
main="", xlab="標準化後の地価", ylab="頻度", cex.axis=1.3, cex.lab=1.2)
# 図3.4
hist(scale(yoko_spdf$lph2010), col="grey", breaks=36,
xlim=c(-2,6), ylim=c(0,300), main="", xlab="標準化後の地価", ylab="頻度",
cex.lab=1.2, cex.axis=1.3)
# 尖度・歪度
library(e1071)
skewness(lph$lph2010.scale)
kurtosis(lph$lph2010.scale)
mean((lph$lph2010.scale-mean(lph$lph2010.scale))^3)/(sd(lph$lph2010.scale)^3)
mean((lph$lph2010.scale-mean(lph$lph2010.scale))^4)/(sd(lph$lph2010.scale)^4) 

3.3 地域属性の差の比較

# SPMデータ
library(RColorBrewer)
spm.shp <- readShapePoints("tma_spm.shp")
# 図3.5
# カラーパレットを指定する方法
pal0 <- c("grey","grey2")
q_spm <- classIntervals(spm.shp$SPM07, n=5, style="quantile")
q_spm_Col <- findColours(q_spm,pal0)
plot(spm.shp, col=q_spm_Col, pch=20, cex=1.5)
legend("topright",fill=attr(q_spm_Col,"palette"),
legend=names(attr(q_spm_Col,"table")),cex=1.2, bty="n")
# spplotを使う方法
spplot(spm.shp, c("SPM07"), cex=1.5, col.regions=brewer.pal(5, "Greys"))
tma_spm <- read.table("tma_spm.csv", sep=",", header=TRUE)
hist(tma_spm$SPM07)
plot(density(tma_spm$SPM07))
# 地区別集計
tma_spm$KCODE <- floor(tma_spm$ID/1000000)
tapply(tma_spm$SPM07, factor(tma_spm$KCODE), length)
tapply(tma_spm$SPM07, factor(tma_spm$KCODE), mean)
tapply(tma_spm$SPM07, factor(tma_spm$KCODE), var)
tapply(tma_spm$SPM07, factor(tma_spm$KCODE), sd)
# 都県別に公示地価地点データを抽出
spm_11 <- tma_spm[tma_spm$KCODE==11,]
spm_12 <- tma_spm[tma_spm$KCODE==12,]
spm_13 <- tma_spm[tma_spm$KCODE==13,]
spm_14 <- tma_spm[tma_spm$KCODE==14,]
# コルモゴロフ・スミルノフ検定
ks.test(spm_11$SPM07, "pnorm", mean(spm_11$SPM07), sd(spm_11$SPM07))
ks.test(spm_12$SPM07, "pnorm", mean(spm_12$SPM07), sd(spm_12$SPM07))
ks.test(spm_13$SPM07, "pnorm", mean(spm_13$SPM07), sd(spm_13$SPM07))
ks.test(spm_14$SPM07, "pnorm", mean(spm_14$SPM07), sd(spm_14$SPM07))
# 図3.6(a)
hist(spm_13$SPM07, col="grey", xlim=c(0.04, 0.10), ylim=c(0,25),
xlab="SPM(mg/m3)", ylab="頻度", main="", cex.lab=1.2, cex.axis=1.3) 
# 図3.6(b)
plot(ecdf(spm_13$SPM07), do.point=FALSE, verticals=TRUE, main="",
lwd=2, cex.axis=1.3, cex.lab=1.2)
z <- seq(0.04, 0.09, by=0.001)
lines(z, pnorm(z, mean=mean(spm_13$SPM07), sd=sd(spm_13$SPM07)), lty=2, lwd=2)
# 分散の差の検定
var.test(spm_13$SPM07, spm_14$SPM07)
# 平均値の差の検定(等分散)
t.test(spm_13$SPM07, spm_14$SPM07, var.equal=T)
# 平均値の差の検定(不等分散)
t.test(spm_13$SPM07, spm_14$SPM07, var.equal=F)
# Wilcoxonの順位和検定
wilcox.test(spm_13$SPM07, spm_14$SPM07)
# 地価データ
# 地区別集計
lph$KCODE <- floor(lph$JCODE/1000)
tapply(lph$lph2010, factor(lph$KCODE), length)
tapply(lph$lph2010, factor(lph$KCODE), mean)
tapply(lph$lph2010, factor(lph$KCODE), var)
tapply(lph$lph2010, factor(lph$KCODE), sd)
# 都県別に公示地価地点データを抽出
lph_08 <- lph[lph$KCODE==8,]
lph_09 <- lph[lph$KCODE==9,]
lph_10 <- lph[lph$KCODE==10,]
lph_11 <- lph[lph$KCODE==11,]
lph_12 <- lph[lph$KCODE==12,]
lph_13 <- lph[lph$KCODE==13,]
lph_14 <- lph[lph$KCODE==14,]
# コルモゴロフ・スミルノフ検定
ks.test(lph_09$lph2010, "pnorm", mean(lph_09$lph2010), sd(lph_09$lph2010))
ks.test(lph_10$lph2010, "pnorm", mean(lph_09$lph2010), sd(lph_10$lph2010)) 
# 分散の差の検定
var.test(lph_09$lph2010, lph_10$lph2010) 
# 平均値の差の検定(等分散)
t.test(lph_09$lph2010, lph_10$lph2010, var.equal=T)
# 平均値の差の検定(不等分散)
t.test(lph_09$lph2010, lph_10$lph2010, var.equal=F)
# Wilcoxonの順位和検定
wilcox.test(lph_09$lph2010, lph_10$lph2010)
# ベイズ的方法
# 横須賀市
lph_14201 <- lph[lph$JCODE==14201,]
# 鎌倉市
lph_14204 <- lph[lph$JCODE==14204,]
# 逗子市
lph_14208 <- lph[lph$JCODE==14208,]
# 三浦市
lph_14210 <- lph[lph$JCODE==14210,]
# 葉山町
lph_14301 <- lph[lph$JCODE==14301,]
length(lph_14201$lph2010)
length(lph_14204$lph2010)
length(lph_14208$lph2010)
length(lph_14210$lph2010)
length(lph_14301$lph2010)
# 確率密度関数
plot(density(lph_14201$lph2010), xlim=c(0, 5e+5), ylim=c(0,3.0e-5),
main="", xlab="地価(円/)", ylab="密度", lty=1, lwd=2,
cex.axis=1.3, cex.lab=1.2)
lines(density(lph_14204$lph2010), lty=2, lwd=2)
lines(density(lph_14208$lph2010), lty=3, lwd=2)
lines(density(lph_14210$lph2010), lty=4, lwd=2)
lines(density(lph_14301$lph2010), lty=5, lwd=2)
rect(320000, 2e-5, 500000, 3e-5,col="white")
segments(350000, 2.9e-5, 400000, 2.9e-5, lty=1, lwd=2)
text(450000,2.9e-5, "横須賀市", cex=1.3)
segments(350000, 2.7e-5, 400000, 2.7e-5, lty=2, lwd=2)
text(450000,2.7e-5, "鎌倉市", cex=1.3)
segments(350000, 2.5e-5, 400000, 2.5e-5, lty=3, lwd=2)
text(450000,2.5e-5, "逗子市", cex=1.3)
segments(350000, 2.3e-5, 400000, 2.3e-5, lty=4, lwd=2)
text(450000,2.3e-5, "三浦市", cex=1.3)
segments(350000, 2.1e-5, 400000, 2.1e-5, lty=5, lwd=2)
text(450000,2.1e-5, "葉山町", cex=1.3)
# 箱ひげ図
boxplot(lph_14201$lph2010, lph_14204$lph2010, lph_14208$lph2010,
lph_14210$lph2010, lph_14301$lph2010,
names=c("横須賀市","鎌倉市","逗子市","三浦市","葉山町"),
cex.axis=1.3, lwd=2)

# 2群間の平均
# 逗子市と葉山町
y1 <- lph_14208$lph2010
y2 <- lph_14301$lph2010
n1<-length(y1) ; n2<-length(y2)
# prior parameters
mu0<-mean(c(y1,y2)) ; g02<-var(y1)
del0<-0 ; t02<-var(y1)
s20<-15; nu0<-1
# starting values
mu<- ( mean(y1) + mean(y2) )/2
del<- ( mean(y1) - mean(y2) )/2
# Gibbs sampler
MU<-DEL<-S2<-NULL
Y12<-NULL
set.seed(1)
for(s in 1:10000) 
{
 ##update s2
 s2<-1/rgamma(1,(nu0+n1+n2)/2, 
       (nu0*s20+sum((y1-mu-del)^2)+sum((y2-mu+del)^2) )/2)
 ##update mu
 var.mu<-  1/(1/g02+ (n1+n2)/s2 )
 mean.mu<- var.mu*( mu0/g02 + sum(y1-del)/s2 + sum(y2+del)/s2 )
 mu<-rnorm(1,mean.mu,sqrt(var.mu))
 ##update del
 var.del<-  1/(1/t02+ (n1+n2)/s2 )
 mean.del<- var.del*( del0/t02 + sum(y1-mu)/s2 - sum(y2-mu)/s2 )
 del<-rnorm(1,mean.del,sqrt(var.del))
 ##save parameter values
 MU<-c(MU,mu) ; DEL<-c(DEL,del) ; S2<-c(S2,s2) 
 Y12<-rbind(Y12,c(rnorm(2,mu+c(1,-1)*del,sqrt(s2) ) ) )
}         
# 図3.9
plot(density(Y12[,1]), col="black", lwd=3, cex.axis=1.3, cex.lab=1.2,
main="", xlab="地価(円/屐", ylab="密度")
lines(density(Y12[,2]), col="grey", lwd=3)
text(3e+5,1e-5,"逗子市", cex=1.3)
text(3e+5,9e-6,"葉山町", cex=1.3)
lines(x=c(2.5e+5,2.75e+5),y=c(1e-5,1e-5), lwd=3, col="black")
lines(x=c(2.5e+5,2.75e+5),y=c(9e-6,9e-6), lwd=3, col="grey")
####
# 以下、作業メモ
#
plot(MU+DEL,MU-DEL)
cor(MU+DEL,MU-DEL)
# par(mfrow=c(1,2),mar=c(3,3,1,1),mgp=c(1.75,.75,0))
plot( density(MU,adj=2),main="",xlab=expression(mu),ylab="density",lwd=2 )
ds<-seq(mu0-sqrt(g02),mu0+sqrt(g02),length=100)
lines(ds,dnorm(ds,mu0,sqrt(g02)),lwd=2,col="gray" )
legend(172280,0.015,legend=c("prior","posterior"),lwd=c(2,2),
col=c("black","gray"),  bty="n")
#
plot( density(DEL,adj=2), main="",xlab=expression(delta),ylab="density",lwd=2 )
ds<-seq(-sqrt(t02),sqrt(t02),length=100)
lines(ds,dnorm(ds,0,sqrt(t02)),lwd=2,col="gray" )
legend(50, 0.015,legend=c("prior","posterior"),lwd=c(2,2),
col=c("black","gray"),  bty="n")
#
quantile(DEL,c(.025,.5,.975))
quantile(DEL*2,c(.025,.5,.975))
mean(DEL>0)
mean(Y12[,1]>Y12[,2])
#
# 多群間の平均
# 省略

3.4 地域格差

# 3.4.1 ジニ係数
library(reldist)
gini(lph$lph2010)
inc <- read.table("inc2006.csv", sep=",", header=TRUE)
library(reldist)
gini(inc$INC2006)
#
library(ineq)
inc.gini <- Lc(inc$INC2006)
inc.gini
# 図3.10
plot(inc.gini, cex.axis=1.3, cex.lab=1.2, cex.main=1.3, lwd=3)
# 3.4.2 地域変動係数
Cv <- sd(lph$lph2010)/mean(lph$lph2010)
Cv
# 3.4.3 地域特化係数
CLpop05 <- read.table("CLpop05.csv", sep=",", header=TRUE)
CLpop05$pop05CL3 <- (CLpop05$POP3 / CLpop05$POP05) / (sum(CLpop05$POP3) /  sum(CLpop05$POP05))
jpn_pref <- readShapePoly("jpn_pref.shp",IDvar="PREF_CODE")
plot(jpn_pref, col="grey")
ID.match <- match(jpn_pref$PREF_CODE, CLpop05$PREFCODE)
ID.match
jpn_CL <-CLpop05[ID.match,]
summary(jpn_CL)
jpn_pref_CL<- spCbind(jpn_pref,jpn_CL)
names(jpn_pref_CL)
summary(jpn_pref_CL)
# カラーパレットの作成
library(maptools)
library(classInt)
pal0 <- c("white","grey","grey2")
# 等量分類
q_pref <- classIntervals(round(jpn_pref_CL$pop05CL3,2), n=5, style="quantile")
plot(q_pref, pal=pal0)
q_pref_Col <- findColours(q_pref,pal0)
plot(jpn_pref_CL,col=q_pref_Col)
#title("地域特化係数(65歳以上人口)")
legend("topleft", cex=1.5, fill=attr(q_pref_Col,"palette"),
legend=names(attr(q_pref_Col,"table")),bty="n")

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-10-27 (水) 14:42:28 (3279d)