空間統計学(統計ライブラリ)　自然科学から人文・社会科学まで

#空間統計学(統計ライブラリ)　自然科学から人文・社会科学まで　p.112

#空間ダービンモデル(Durbin):SDM

library(dplyr)

#各都道府県の距離データ(https://www.gsi.go.jp/KOKUJYOHO/kenchokan.html)
d=read.csv("C:/Users/kozakai ryouta/Documents/都道府県の距離データ.csv")

name=as.character(d[,1]);d=d[,-1]

d=as.matrix(d)

rownames(d)=name

d[d=="" | is.na(d)>0]=0

d=matrix(as.numeric(d),ncol=ncol(d))

rownames(d)=c(name);colnames(d)=c(name)

d=d+t(d);d=d/10

d=d[-nrow(d),-ncol(d)]

W=1/d;diag(W)=0

#y:交通事故死傷者数、x1:人口密度、x2:道路投資額、x3:道路改良率、x4:トラック・バス・乗用車登録台数

data=data.frame(prefecture=c("Hokkaido","Aomori","Iwate","Miyagi","Akita","Yamagata","Fukushima","Ibaragi","Totigi","Gumma","Saitama","Tiba","Tokyo","Kanagawa","Niigata","Toyama","Ishikawa","Fukui","Yamanashi","Nagano","Gifu","Shizuoka","Aichi","Mie","Shiga","Kyoto","Oosaka","Hyogo","Nara","Wakayama","Tottori","Shimane","Okayama","Hiroshima","Yamaguchi","Tokushima","Kagawa","Ehime","Kouchi","Hukuoka","Saga","Nagasaki","Kumamoto","Ooita","Miyazaki","Kagoshima"),x1=c(6364,7135,7266,7186,6274,6497,6740,6687,6917,6725,7928,7403,13500,8537,7315,6730,8875,7589,7249,6319,7673,7270,8149,6385,7049,10861,11548,9986,6728,7596,7324,6166,5889,8170,5147,6997,6460,6657,7979,7301,6564,8866,6928,6983,6016,7688),x2=c(90547,12634,20101,15486,12471,12446,27026,19486,19919,16611,42768,37630,147481,55648,27569,11709,13660,11161,10092,22957,20135,36377,48255,19562,9889,16925,146758,75350,16405,13007,8005,10334,22146,23029,15488,9905,7407,14795,12387,31172,8151,11117,14634,11504,11415,12546),x3=c(19.7,20,16.3,20.6,18.1,21.2,14.2,9.6,28.8,16.3,15.5,28.4,49.1,31.8,14.6,35.3,30.8,32.9,24.2,11.4,18.3,19.6,21.4,14.6,21.4,23.6,40.9,24.8,10.7,16.8,26.6,10.2,11,16.3,23.3,10.8,21.9,16.2,11.8,17.7,24.7,17.5,18.5,20.9,25.2,32.5),x4=c(640.9,128.2,113.8,171.7,97.4,134.4,166.5,225.1,187.5,212,337.3,319.4,1584.3,559.2,234.3,121.4,120,93.5,93.2,238.1,255.1,419.1,825.1,160.5,91.8,233.6,807.5,404.6,85.8,111,44.1,52.9,156.6,256.1,126.8,71.4,81.5,108.1,72.7,367,70.6,89.6,136.6,97.5,96.3,119.4),y=c(44523,10395,8669,12808,7640,7094,17556,20320,18715,16755,34601,27811,88406,48030,18213,8839,12086,8871,8961,15280,17837,37421,50998,14809,12610,35478,75497,55464,7916,14259,6244,5400,19015,36771,13827,9372,11088,9960,8516,51175,10384,10652,15987,10216,8460,12932))

Y=data$y;X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2","x3","x4")]))

X_star=cbind(X,W%*%X)

n=nrow(data)

beta0=solve(t(X_star)%*%X_star)%*%t(X_star)%*%Y

betad=solve(t(X_star)%*%X_star)%*%t(X_star)%*%(W%*%Y)

epsilon0=Y-X_star%*%beta0

epsilond=W%*%Y-X_star%*%betad

loglik=function(p){

return(-n*log(sum((epsilon0-epsilond*p)^2)/n)/2+log(det(diag(rep(1,nrow(W)))-p*W)))  

}

P=0.3

ite=1000

eta=0.001

h=0.01

for(l in 1:ite){

P=P+eta*(loglik(P+h)-loglik(P))/h  

print(loglik(P))

}

beta=beta0-P*betad

sigma=sum((epsilon0-epsilond*P)^2)/n

pre=P*W%*%Y+X_star%*%beta

plot(c(1:length(Y)),Y,xlab="index",ylab="predict",xlim=c(1,length(Y)),ylim=c(min(c(Y,pre)),max(c(Y,pre))),col=2,type="o")

par(new=T)

plot(c(1:length(Y)),pre,xlab="index",ylab="predict",xlim=c(1,length(Y)),ylim=c(min(c(Y,pre)),max(c(Y,pre))),col=3,type="o")