####################################
# CSVファイルの読み込み
####################################
#Mac読み書きディレクトリを定義
#ソースを使い回す場合は書き換えてください
basedata<-"/Users/saity/Dropbox/R/basedata/"
output<-"/Users/saity/Dropbox/R/output/"

#CSVファイルの読み込み
crime<-read.csv(paste(basedata,"crime_prefecture.csv",sep=""), row.names=1)


####################################
# データサマリー分析
####################################

#ヘッダー情報
head(crime)

#データサマリー出力
summary(crime)

#ヒストグラム（棒グラフで犯罪数の分布を知る）
#第2引数で区切り幅を指定出来る
hist(crime$arrest_number)
hist(crime$arrest_number, 20)

#散布図を描画(認知・検挙・人口密度・登録外人・外人・離婚・自殺…)
#こんな風にいっぺんにやると、散布図がゴチャっとしてわかりにくい
pairs(crime[c(6:19)])

#工夫しながら二つに分割しましょう
pairs(crime[c(6:12)])
pairs(crime[c(6,13:19)])

#相関マトリクス
#1は正の相関 -1は負の相関
#1や-1に近いデータ程散布図の散らばりが少ない(0.8以上は相関がかなり高い)
cor(crime[3:19])

####################################
# χ二乗検定
####################################
#検挙率と離婚カイ二乗
table1<-xtabs(~arrest_number+divorces,data=crime)
write.csv(table1,paste(output,"crime&divorces_hist.csv",sep=""))
chisq.test(table1)

#検挙率と自殺カイ二乗
table2<-xtabs(~arrest_number+alcohol,data=crime)
write.csv(table2,paste(output,"crime&suicide_hist.csv",sep=""))
chisq.test(table2)

#消費アルコール量と失業率
table3<-xtabs(~crime_rate+crime_cognition_number,data=crime)
write.csv(table3,paste(output,"crime&suicide_hist.csv",sep=""))
chisq.test(table3)


####################################
# クラスター分析
####################################
#クラスタリングに必要なrpartライブラリの読み込み
library(rpart)

#クラスタリング：距離を求める[9～19列目]
#データに不公平が無いように[data-平均/標準偏差＝標準得点]として出力している(プラス/マイナスで表される)
(sdata1<-scale(crime[6:19]))
(dist1<-dist(sdata1))


####################################
# 重回帰分析
####################################

#基準・説明変数間の相関をチェック
cor(crime[6:19])

#重回帰分析
#lm(基準変数 ~ 説明変数1 + 説明変数2 + 説明変数3 … )
# ***のように、*がついている数が多い程高優位
multiple_regression1<-lm(arrest_number ~ area + population_density + foreigners + marriages + divorces + suicide + police + unemployment + sunsine_of_h + rain_of_d + income + alcohol + sleep ,crime)
summary(multiple_regression1)

#関係のなさそうなデータを抜くと、新たに優位となるデータが出てくる
multiple_regression1<-lm(arrest_number ~ population_density + marriages + divorces + suicide + police + sunsine_of_h + rain_of_d + income + alcohol ,crime)
summary(multiple_regression1)


#多重共線性チェック（10を超えていると危ない)
diag(solve(cor(crime[6:19])))
#（10は超えないが）高い数値である検挙率・婚姻・離婚を省く
diag(solve(cor(crime[c(9,12:19)])))

#標準化回帰係数
sy<-scale(crime$arrest_number)
sx<-scale(crime[c(9,12:19)])
rg2<-lm(sy~sx)
summary(rg2)

#ステップワイズ法
library(MASS)
rg0<-lm(arrest_number ~ crime_rate + foreigners + marriages + divorces + suicide + police + unemployment + sunsine_of_h + rain_of_d + income + alcohol + sleep ,crime)

rg0<-lm(arrest_number ~ crime_rate + foreigners + marriages + divorces + suicide + police + unemployment + sunsine_of_h + rain_of_d + income + alcohol + sleep ,crime)

rg1<-stepAIC(rg0)


####################################
# ロジスティクス分析
####################################
#ロジスティック回帰
attach(crime)
Aclass<-arrest_number>=350

lrg1<-glm(Aclass~ population_density + foreigners + marriages + divorces + suicide + crime_cognition_number + police + unemployment + sunsine_of_h + rain_of_d + income + alcohol + sleep, binomial)

lrg1<-glm(Aclass~ foreigners + marriages + divorces + suicide + police + unemployment + sunsine_of_h + rain_of_d + income + alcohol + sleep, binomial)

#確率予測
prob<-fitted(lrg1)
write.csv(prob,paste(output,"prov.csv"))

####################################
# @TODO:警官数と人口密度の相関を調べてみる
####################################
# @TODO:検挙数と検挙率で、高い所・低い所の特徴を調べて再分析してみる
####################################