epiDisplayの使い方
多変量ロジスティック回帰分析結果を出力できるepiDisplayの使い方。
インストールはいつもの感じでパッケージをインストールする。
install.packages("epiDisplay")ライブラリを呼び出す
library(
epiDisplay)オッズ比と、調整済みオッズ比と95%信頼区間を出力
logistic.display(res_forward)
調整済みオッズ比の結果だけ出力
logistic.display(res_forward, simplified=TRUE)
R Commanderのインストール
Rには、コマンド入力ではなくメニュー等からの選択で実行する「R Commander」とゆうGUIがある。
インストール方法は
install.packages("Rcmdr")
でパッケージをインストールする。
ミラーサイトを選択して実行。
実行後
library(Rcmdr)
で実行すると、不足しているパッケージが表示され、インストールするかどうかを聞かれるのでインストールする。
インストールが終了すると、Rcmdrパネルが表示されます。
ちなみに一度Rcmdrパネルを閉じると、もう一度library(Rcmdr)を入力しても表示されず、Ｒを再起動したら表示される。library()だから２回目以降は効いていないのか？たぶんなにか方法があるんだろうと思うので、判ったら追記しよう。
RとQGISとPostgreSQLと
FOSS4G Advent Calendar 2017 18日目の記事です。
毎年皆さま新しい技術記事のなか、今さらながらＲを始めたのでQGISのプラグインから使ってみた記事です。
R：オープンソースの統計解析向けプログラミング言語。初心者です。
QGIS：オープンソースのGISソフト。お世話になってます。
PostgreSQL：オープンソースのデータベースソフト。お世話になってます。
Python：いろいろ使えるスクリプト言語。チューニング次第で早くなったり遅くなったり。
で、この方々の関係がどうなっているかとゆうと
R - QGIS
プロセッシングでＲスクリプトを利用できる。
https://docs.qgis.org/2.18/ja/docs/training_manual/processing/r_intro.html
R- PostgreSQL
RPostgreSQLパッケージで利用できる。
R- Python
今はRを使わないでPythonで統計をやるのが流行りらしい。
Python からRを利用できるライブラリのPypeRとゆうものがあるらしい。
レイヤとしてテーブル内容を読み込めるし、プラグインからSQLを実行できる。
ドライバモジュール使用で接続できる。
らしい。
とゆうことでQGISのプラグインをPythonで作成して、PyperでＲを実行。R内でRPostgreSQLを使ってデータベースを更新してQGISに更新結果を反映してみる試み。
QGISとＲとPostgreSQLはインストール済み。
データは国土数値情報の市区町村データと統計局estatの統計データから出生数・人口総数・婚姻数等をPostgreSQLに突っ込む。
DB 名：rsample
市区町村名テーブル：shikuchouson。これのbornに出生数、born2に推定出生数、born_diffに差異を入れる。
統計データテーブル：tbl2017_a。
a1101 ：総人口、a4101：出生数、a9101：婚姻件数
これをRで出生数を目的変数、人口総数、婚姻数を説明変数とし、重回帰分析して出生数の実数、推計数、推計差分をQGISプラグインから実行してテーブルを更新して表示してみる。
まずQGISプラグインを作成する。作り方は矢崎さんの記事などを参考。
QGISプラグインの作り方（パッケージ生成編） - Qiita
私の環境(QGIS2.18)だとresources_rc.pyをresources.pyにリネームしないとプラグインが壊れていますとか言われた。GUIはこんな感じ。この推計出生数ボタンを押された時にＲスクリプトを実行してＤＢを更新させる。
PYTHONHOME： QGISで使用するpythonのインストール先
PYTHONPATH：pythonのsite-packagesのインストール先
R_HOME : Rのインストール先
PATH：R_HOMEを追加
そしてPYTHONHOMEのpythonにRのライブラリPypeRをpipでインストールする。対応するpythonのsite-packagesにインストールされるはず。
pip install pyper
インストールが正常にいき、パスがちゃんと通っていればQGISのプラグイン>pythonコンソール上でimportできるはず。エラーがでたら失敗です。pythonのバージョンいろいろ入れている人はパス関係を確認してください。私はここで、はまりました。
最後にRにPostgreSQLに接続するためのRPostgreSQLパッケージをインストール。
Rを起動して
install.packages("DBI")
install.packages("RPostgreSQL")
これで使用準備は終わり。後はＲのスクリプトを作成する。今回のスクリプトは以下のようなものを born.Rとして保存。
#RPostgreSQLライブラリ呼び出し
require("DBI")
require("RPostgreSQL")
# DB連結
con <- dbConnect(PostgreSQL(), host="localhost", port=5434, dbname="rsample", user="postgres", password="postgres")
# データセット作成
dataset <- dbGetQuery(con,"SELECT * FROM tbl2017_a")
# lmで重回帰分析
ans <- lm(a4101~a1101+a9101, dataset)
# パラメータ取得
intercept <- as.character(ans$coefficients["(Intercept)"])
a1101 <- as.character(ans$coefficients["a1101"])
a9101 <- as.character(ans$coefficients["a9101"])
# 推計式作成
formura = paste(intercept,"+",a1101,"*a1101 +",a9101,"*a9101")
# 推計値を入力するSQL
sql = paste("UPDATE shikuchouson SET born2=a.born2 FROM",
"(SELECT scode,round(", formura, ",0) AS born2 FROM tbl2017_a) a",
"WHERE a.scode=TO_NUMBER(n03_007, '99999');")
# 実際値と推計値の差異を入力するSQL
sql2 = "UPDATE shikuchouson SET born_diff=abs(born-born2);"
# SQL実行
dbGetQuery(con,sql)
dbGetQuery(con,sql2)
#---------------------------------------------------------------------------------------
# Ｒ推計値セット
#---------------------------------------------------------------------------------------
import pyper as py
def setR(self):
success = True
try:
r = py.R()
r("source(file='C:/tmp/R/QGIS/born.R')")
QtGui.QMessageBox.information(self.dlg, u'確認', u'Ｒ推測値セットしました', QtGui.QMessageBox.Close)
except:
QtGui.QMessageBox.critical(self.dlg, u'エラー', u'Ｒ推測値セットに失敗しました', QtGui.QMessageBox.Close)
success = False
import pyper as py でpyperをインポート。pyを別名にする。
r = py.R() でRのインスタンスを作成
r("source(file={ファイル名})") でRのスクリプトを実行する。
これはR側で全部処理した場合ですが、python側で処理したものをRとやりとりする場合は
r.assign('data', test) # rにtestをdataとして渡す
res = r.get("result") # rのresultをresとして受け取る
のような形で受け渡しできます。
推計値で色分け。
推計値差異で色分け。濃いとこが差異多い。
ちなみに北海道では札幌市を各区別としたら、旭川市が一番出生数が多いらしい。
Ｒは地図用のパッケージいろいろあるみたいですね。来年勉強します。
ロジスティック回帰分析（その１）
ロジスティック回帰分析とは発生確率を予測する手法。
基本的な考え方は線形回帰分析と同じだが、予測結果が０～１の値を取る。
０%（０．００）～100%（１．００）のような確率で予測を行う場合に有意な手法。
Ｒではglm()で以下のように実行できる
ans = glm(Y~X1+X2+・・・Xn, data=DATA, family=binomial(link="logit"))
summary(ans)
線形解析分析と同様に解析結果は β0+β1*X1+β2*X2・・・βn*Xnの形になりますが、推測値として０～１の範囲内に収まる確率値を得るために対数をとり、以下の式になります。
log(Y/(1-Y)) = β0+β1*X1+β2*X2・・・βn*Xn
ここでYは目的変数が発生する確率であり１－Ｙは発生しない確率を意味します。
分解すると
ｔ＝β0+β1*X1+β2*X2・・・βn*Xn
log(Y/(1-Y)) = t
Y/(Y-1) = exp(t)
et = exp(t)
Y = et /(1+et)
となります。
重回帰分析（その２）
さて、桐生６日目最終日。本当は５日目までのデータをとりたかったが、竹井選手が優勝戦に乗れず、やる気がなくなったので４日目までのデータで重回帰分析してみた。
目的変数 result （着順）
説明変数１ shinnyu （進入コース）
説明変数２ motor_2win （モータ２連帯率）
説明変数３ bort_2win （ボート２連帯率）
ans = lm(result~shinnyu+motor_2win+bort_2win, dataset)
ans
推定着順 = 4.234239 + 0.380960 * shinnyu -0.056836 * motor_2win -0.001379 * bort_2win
これを元に最終日優勝戦の推定着順を求めてみると
resultが小さいほど１着に近いと考える。
まず１号艇頭は固い＆６号艇は要らない。２，３，４，５は団子とゆう、まぁあまり面白くない推定になった。まぁ２連帯率だけだからね。
ちなみに展示走行で５号艇がスタート気張りそう＆モータ勝率で２着固定にして取れました。
１－＞５－＞２ ￥2420。大瀧選手おめでと～！
