$%2017/06/12
\newcommand\compl[1]{{#1^\mathtt{c}}}
\newcommand\pare[1]{{(#1)}}
\newcommand\Pare[1]{\left(#1\right)}
\newcommand\curl[1]{\{#1\}}
\newcommand\Curl[1]{\left\{#1\right\}}
\newcommand\squa[1]{[#1]}
\newcommand\Squa[1]{\left[#1\right]}
\newcommand\abs[1]{\lvert#1\rvert}
\newcommand\Abs[1]{\left\lvert#1\right\rvert}
\newcommand\floor[1]{\lfloor#1\rfloor}
\newcommand\Floor[1]{\left\lfloor#1\right\rfloor}
\newcommand\ceil[1]{\lceil#1\rceil}
\newcommand\Ceil[1]{\left\lceil#1\right\rceil}
\newcommand\setN[0]{\mathbb{N}}
\newcommand\setZ[0]{\mathbb{Z}}
\newcommand\setQ[0]{\mathbb{Q}}
\newcommand\setR[0]{\mathbb{R}}
\newcommand\setC[0]{\mathbb{C}}
\newcommand\sfrac[2]{#1/#2}
\newcommand\od[2]{\frac{d#1}{d#2}}
\newcommand\pd[2]{\frac{\partial#1}{\partial#2}}
\newcommand\sod[2]{\sfrac{d#1}{d#2}}
\newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}}$
$\newcommand\L{\mathcal{L}}
\newcommand\R{\mathcal{R}}
\newcommand\A{\mathscr{A}}
\newcommand\B{\mathscr{B}}
\newcommand\G{\mathscr{G}}
\newcommand\M{\mathbb{M}}
\newcommand\I{\mathscr{I}}
\newcommand\N{\mathscr{N}}
\newcommand\setOR[0]{\overline{\mathbb{R}}}
\DeclareMathOperator{\vol}{vol}
\DeclareMathOperator\supp{supp}$

# R言語で統計処理をお手軽にやってみよう

## 教科書とサポートサイト
データ解析のための統計モデリング入門をもとに進める．
以下はサポートサイト
https://kuboweb.github.io/-kubo/ce/IwanamiBook.html

Rの文法事項は以下を見れば大概のっている．
http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html

## インストール
Windowsの場合，以下サイトに従ってインストール
https://qiita.com/FukuharaYohei/items/8e0ddd0af11132031355

## 基礎的なコマンド
カレントディレクトリを変更：
setwd("C:/Users/honyarara/Desktop")

txt,csv,xlsxファイルを読み込み：
https://htsuda.net/stats/dataset.html

直接入力によりデータを作成：
https://www.cis.doshisha.ac.jp/mjin/R/02.html

例えば，
c(20,20,30,10,30,10)
は20,20,30,10,30,10の要素からなるデータを表す．

data<-c(20,20,30,10,30,10)
にて，dataという変数を用意し，保存しておくことになる．

http://cse.naro.affrc.go.jp/takezawa/r-tips/r/12.html
にはそのほか典型的なデータの作り方が書いてある．
http://takenaka-akio.org/doc/r_auto/series.html
にも書いてある．

ベクトルの要素を直接取り出したり，変更するときには[]を使う．

data*2
で各要素が２倍になったデータをかえす．

data+3
で各要素に3を加えたデータをかえす．

f <- function(m) m*m
でmを受取m*mをかえす関数f(m)を定義する．

sapply(data,f)
で各要素にf(m)を作用させたデータをかえす．

idata <- numeric(30)+1
で30個の１が並んだデータを作成する

以下のようにfor文を使うこともできる．
appendは末尾に指定データを加えたデータベクトルを作ることを意味する．
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/30.html

数字を丸めるには，例えば以下の命令を使う．
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/37.html

data
と打つと中身を直接表示．

length(data)
でデータの個数を表示

summary(data)
で標本平均，最小，最大値，四分数位などを表示．

table(data)
で度数分布を表示

hist(data, breaks = seq(-5, 35,10))
でヒストグラム表示

mean(data),var(data), sd(data)は標本平均，標本分散，標本標準偏差をそれぞれ出力する．

## 疑似乱数
ポアソン分布に従う平均3.5の乱数を２００個集めたデータを作成：
data2 <- rpois(200, 3.5)

ヒストグラムで表示する場合：
hist(data2, breaks = seq(-0.5,13,1))

たとえば，
f <- function(m) mean(rpois(50, m))
で実行のたび，「平均mのポアソン分布に従う50要素からなるデータ」を作成し，その平均をかえす関数fを作成する．
実行のたびにfの値は異なることに注意する．

idata <- numeric(3000)+3.5
f <- function(m) mean(rpois(50, m))
result = sapply(idata, f)
hist(result, breaks = seq(2.5,4.5,0.1))

でfの実行結果を3000個並べたデータが作れる.supply関数がmap関数的機能を実現していることに注意．

同じことは
result = rep(mean(rpois(50, 3.5)), times = 3000)
hist(result, breaks = seq(2.5,4.5,0.1))
でいけそうな気もするが，この場合は一回作った乱数を使いまわされてしまう(repは値をコピーして並べるという意味である)．

この問題を回避するには，例えば，以下のコードを使うことが考えられる．
f <- function() mean(rpois(50, 3.5))
result <- numeric() #空ベクトルを作成
for(i in 1:3000) result <- append(result,f())
hist(result, breaks = seq(1.5,5.5,0.1))
難度もresultのオブジェクト自体を上書きしているので美しくはないが．．．

## 3.2節csvファイルの処理

d <- read.csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv")
のように，read.csvファイルでcsvファイルを取り込んだデータベクトルを作成．dに代入．
教科書と同じ結果が得られるようにするには，
d <- read.csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv", stringsAsFactors=TRUE)
とするように細工が必要．

d$x
でラベルがxの列の中身を表示．
d$y
ならラベルがyの列の中身を表示．

read.csv()関数は文字を含む列があった場合，その列のデータを因子(factor)というクラスのデータとして取り込む．
サポートサイトのdata3a.csvでいえば，
d$f
でfactor型がCとTの2水準(level)からなる値で構成されていることが確認できる．

class(d)
でdのクラスを表示．
class(d&y)
でラベルがyの列には何のクラスのデータが代入されているか表示．
#### 手動でデータを作る場合
手動で打ち込む場合は以下の手順で．

#### 手動で表データを作る場合
行列にnames属性を与えることができる．

x <- matrix(1:6, nrow=2, ncol=3)
colnames(x) <- c("left","center","right")
rownames(x) <- c(1,2)

スマートにしたいなら，
行の数を返す関数nrow()を使って，
x <- matrix(1:24, nrow=4, ncol=3)
colnames(x) <- c("left","center","right")
rownames(x) <- c(1:nrow(x))
と書けば，列の名前をc(1,2,3,4)と直打ちしなくて済むようになる．

x[2,3]などでデータの各要素に直接アクセスできる．
x[2,3]<- 30などで値を直接書き換えてしまうこともできる．

x[2,3] <- "hoge"により，データのひとつだけを文字列にしようとすると，全部のデータが文字列にされてしまう．

行列のクラスはあくまで"matrix" "array" であって，データフレームとは異なることに注意すること．

d2 <- data.frame(x)
と，行列をデータフレーム型に変換することで，data.frameのオブジェクトを作成できる．

文字列と数字が混在している列のクラスは
文字列と判定されるようである．
一つのセルを"character"型にすると，列丸ごと"character"型に変換されているっぽい．

このことは，Rの型に関するルールに由来している．
as.integer()関数などで列まるごとキャストできる
このコードでは，正確にいうと，キャストして新しく作った列データを，データフレームの該当列に上書きコピーしている．

## 3.4節
### 箱ひげ図の処理
d <- read.csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv", stringsAsFactors=TRUE)
plot(d$f,d$y)
で箱ひげ図を作図可能
### 3.4.1線形予測子
ある個体$i$のサイズが$x_i$でその種子数が$y_i$というデータが得られているとする．
説明変数$x_i$に対して，種子数の平均$\lambda_i$が
\begin{align}
\lambda_i = \exp(\beta_1+\beta_2x_i)
\end{align}
であると仮定する．ただし，$\beta_1,\beta_2$はパラメータ．

これを変形すると以下が得られる．

\begin{align}
\log\lambda_i = \beta_1+\beta_2x_i
\end{align}

上記右辺を線形予測子と呼ぶ．
上記左辺をリンク関数と呼ぶ．
$\log\lambda_i$の形になっているリンク関数を対数リンク関数とよぶ．
### 3.4.2あてはめ

与えられたパラメータ$\beta_1,\beta_2$に対数尤度$\log L$は以下で得られる．
\begin{align}
\log L(\beta_1,\beta_2)=
\sum_i
\log
\frac{\lambda_i(\beta_1,\beta_2)^{y_i}\exp(-\lambda_i(\beta_1,\beta_2))}{y_i!}
\end{align}

各パーツの意味を以下に図示する．

与えられたデータ$x_i,y_i$をもとに，
対数尤度$\log L$を最大化する$\beta_1,\beta_2$を求めるには以下のコマンドを利用する．

fit <- glm(y ~ x, data = d, family = poisson)

summary(fit)で詳細な結果が表示される．

標準誤差(Standard error, SE)：推定値$\hat{\beta}_1$,$\hat{\beta}_2$ のばらつき
z value: 最尤推定値をSEで除した値
→尤度が単峰で鋭ければ鋭いほど大きな値
→ただし，$(\beta_1,\beta_2)$空間における演算をしているのではなく，
$\beta_1$に関するz valueを求めたいときには，$\beta_2$を最尤推定値に固定して，$\beta_1$をずらしたときの尤度を考えているっぽい？
→p52の図3.6は２つの分布を同じグラフに無理やりのせているっぽい？

Pr(>|z|)は定義が正確に書いていない．
glm()に限定した場合の定義しか書かれていないので，
「統計学的な検定に使えそうな概念」程度でスキップして良いのでは？

最大対数尤度を，「当てはまりの良さ」とよぶ．
最大対数尤度は，対数尤度$\log L(\beta_1,\beta_2)$が最大になっている箇所における尤度のことである．
最大対数尤度を評価するには以下のコマンドを用いる．

d <- read.csv("https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/poisson/data3a.csv", stringsAsFactors=TRUE)
fit <- glm(y ~ x, data = d, family = poisson)
logLik(fit)
'log Lik.' -235.3863 (df=2)

### 3.4.3 ポアソン回帰モデルによる予測

plot(d$x, d$y, pch = c(21,19)[d$f])
xx <- seq(min(d$x), max(d$x), length = 100)
lines(xx, exp(1.29 + 0.0757 * xx), lwd = 2)

seq命令でｘ軸方向に量子化をかけているベクトルをつくり，
ベクトルの四則演算により，y軸方向のベクトルを作り，それをlinesでplotしていると思われる
（散布図的なものの点列を線でむすぶ命令っぽい）．
実際にデータの中身を見てみると，
であり，確かにそのようになっている（アフィンなグラフを描く場合には２点でも結局同じ表示になるので，でてくるもの自体はかわらない）．

#### todo