絶対に描いてはいけないグラフ入りスライド24枚

Japan.R 2014で絶対に描いてはいけないグラフ入りスライド12枚と題して発表してきました。発表時間は5分という制約があったので、あえなく当日のスライドは半分に。本来発表に使いたかったものを含めてスライド数は倍にしての公開です。

絶対に描いてはいけないグラフ入りスライド24枚 from itoyan110

Rでオブジェクトの存在を確認する方法

Rでは変数名の上書きをうっかりしないよう、プログラマは変数を意識してコードを書かなければならない。そこで、Rにおいてオブジェクトが既に存在しているか確認する方法についていくつか挙げてみる。

1. ls()関数
• オブジェクトの一覧を表示する関数である。search()関数を利用するとサーチパス上にある名前空間の一覧をすべて見ることができる。

  > search()
  [1] ".GlobalEnv"        "package:stats"     "package:graphics" 
  [4] "package:grDevices" "package:utils"     "package:datasets" 
  [7] "package:methods"   "Autoloads"         "package:base" 
  

  単にls()とすると、".GlobalEnv"の中にあるオブジェクト一覧が表示される。

  > a < -5
  > ls()
  [1] "a"
  

  他の名前空間にあるオブジェクト一覧を知りたいというときは、search()関数で分かった順番を基に、pos引数を使って指定すればよい。以下の例では、2番目の値にあったstatsパッケージのオブジェクト一覧を見ていることになる。

  > head(ls(pos=2))
  [1] "acf"        "acf2AR"     "add.scope"  "add1"       "addmargins"
  [6] "aggregate" 
  

2. exists()関数
• オブジェクトの存在をダイレクトに確認するのであれば、この関数を利用するとよい。

Rでコイン投げの分析を一捻り

Japan.R 2013でコイン投げの分析を一捻りと題して発表してきました。ノンパラメトリックな手法を学習しているうちに見つけたのでご紹介いたしました。

コイン投げの分析を一捻り (Japan.R 2013 LT) from itoyan110

Rでノンパラメトリック法 1

Nagoya.R #10にて、「Rでノンパラメトリック法 1」と題して発表してきました。1冊の本が終わるか私が納得するまでシリーズとして続ける予定なので、2回目以降もきっとあるはずということで1とスライドの題名につけました。今回はノンパラメトリックな手法として符号検定とウィルコクソンの符号付き順位和検定を取り上げています。手計算とRでの実行結果が一致することを確認しながら、理屈を知ろうという内容になっております。

Rでノンパラメトリック法 1 from itoyan110

Writing R Extensions (version 2.15.2) 日本語訳を公開しました!

R 2.15.2 用に書かれたWriting R Extensionsを日本語に翻訳したドキュメント「Rの拡張を書く」を公開しました。slideshareからダウンロードできます。

Rの拡張を書く （R 2.15.2） from itoyan110

Rを含め、統計学の更なる普及を期待して作成いたしました。是非ご活用ください。

Rによる統計解析ハンドブック 第8章演習問題解答

This post is solutions of the chapter 8 of `A Handbook of Statistical Analyses Using R, Second Edition'. Chapter 6 and 7 will be posted later.

このポストはRによる統計解析ハンドブック(第2版)の第8章の解答になります。第6, 7章については後日掲載する予定です。

1

Although it is difficult to discern from the histogram, result of which kernel band is 3000 indicates that there are 6 peaks in density function. As kernel band gets lerger, we can't discern all of 6 peaks because they are mixed.

ヒストグラムからは判別しにくいが、カーネル幅が3000のときの結果から、密度関数からは峰が6つあることが示唆される。カーネル幅が大きいと峰が混ざるため、すべて判別できなくなる。

data("galaxies", package="MASS")

oldpar <- par(mfrow=c(2, 2))
do <- function(dat, ker, wid){
  h.info <- hist(dat, plot=FALSE)
  d.info <- density(dat, kernel=ker, width=wid)
  yr <- range(c(0, h.info$density, d.info$y))
  tit <- paste(ker, "kernal width = ", wid, collapse=" ")
  hist(dat, probability=TRUE, ylim=yr, main=tit)
  lines(d.info, lwd=2)
}

do(galaxies, "gaussian", 3000)
do(galaxies, "gaussian", 5000)
do(galaxies, "triangular", 3000)
do(galaxies, "triangular", 5000)
par(oldpar)

2

Data can be seen around (20, 10), and we can see that death rate is constant with birth one.

データは(20, 10)近辺に多く見られ、0.001の等高線からは出生率に関係なく死亡率は一定に近い傾向があることが見て取れる。

library(KernSmooth)
data("birthdeathrates", package="HSAUR2")
bdr.d <- bkde2D(birthdeathrates, bandwidth=sapply(birthdeathrates, dpik))
contour(bdr.d$x1, bdr.d$x2, bdr.d$fhat,
        xlab=names(birthdeathrates)[1], ylab=names(birthdeathrates)[2], 
        main="estimated density of birthdeathrates")

3

Note that plot.Mclust doesn't reflect option such as col and ylim. Graphs show that there is modality for female, and tendency that male develops schizophrenia earlier than female (Note that vertical axis is different between two graphs).

plot.Mclustはcol、ylimなどのオプションを指定しても反映されないことに注意。グラフを見ると、女性は二峰性があり、男性は女性よりも若いころに発症する傾向がある(縦軸の値が異なることに注意)。

library(mclust)
data("schizophrenia", package="HSAUR2")
# previous consideration (事前の考察)
boxplot(age~gender, data=schizophrenia, main="age distribution")
# split data(データの分割)
female <- subset(schizophrenia, gender=="female", select=age)
male <- subset(schizophrenia, gender=="male", select=age)
mc.fem <- Mclust(female)
mc.mal <- Mclust(male)

oldpar <- par(mfrow=c(1, 2))
plot(mc.fem, what="density", xlim=c(0, 70), ylim=c(0, 0.06))
plot(mc.mal, what="density", xlim=c(0, 70), ylim=c(0, 0.06))
par(oldpar)

Rによる統計解析ハンドブック 第5章演習問題解答

This post is solutions of the chapter 5 of `A Handbook of Statistical Analyses Using R, Second Edition'.

このポストはRによる統計解析ハンドブック(第2版)の第5章の解答になります。

1

Estimated values are larger than actual ones for pair of beef and low or that of cereal and high. Otherwise these are lower than that. B, C, L, H in the below graph stand for Beef, Cereal, Low and High.

牛肉で低たんぱく質、穀物で高たんぱく質は推定値が大きめになり、それ以外は推定値が低めに出ている。グラフ中のB、C、LとHはそれぞれBeef、Cereal、LowとHighの頭文字をとったものであることに注意。

data("weightgain", package="HSAUR2")
resid <- weightgain$weightgain - aov(weightgain ~ source + type, data=weightgain)$fitted.values
plot(resid, type="n", main="residuals")
lab <- paste(abbreviate(weightgain[,1],minlength=1),
             abbreviate(weightgain[,2],minlength=1), sep="")
text(resid, label=lab)

2

The result of plot.design shows no effect on gender and learner. So I conducted anova only with race and school. Race is only affects absent days in conclusion. Note that if anova is conducted with all vatiable, f-value of learner is smaller than 1 and it should be pooled.

plot.designからgender、learnerは関係なさそうとなる。そこで、人種、学年のみを使って分散分析を行った。結論としては欠席日数に効いているのは人種であると言える。なお、全変数で分散分析をかけると、learnerのF値が1より小さくプーリングすべきと考えられる。

data("schooldays", package="HSAUR2")
plot.design(schooldays)
aov.res <- aov(absent ~ race * school,  data=schooldays)
anova(aov.res)
interaction.plot(schooldays$race, schooldays$school, schooldays$absent)

> anova(aov.res)
Analysis of Variance Table

Response: absent
             Df  Sum Sq Mean Sq F value    Pr(>F)    
race          1  2645.7 2645.65 12.4615 0.0005561 ***
school        3  1325.5  441.85  2.0812 0.1052455    
race:school   3  3738.2 1246.08  5.8693 0.0008221 ***
Residuals   146 30996.7  212.31                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

3

By executing below code, we find that the hypothesis of the uniformity of mean values are rejected between any two groups.

以下のコードを実行すると、任意の2グループ間で、平均値の同等性の仮定は棄却されることが分かる。

data("students", package="HSAUR2")
students.manova <- manova(cbind(low, high) ~ treatment, data=students)
summary(students.manova, test="Hotelling-Lawley") # p-value is very small!
summary.aov(students.manova)

summary(manova(cbind(low, high) ~ treatment, data=students, 
        subset = treatment %in% c("AA", "C")))
summary(manova(cbind(low, high) ~ treatment, data=students, 
        subset = treatment %in% c("C", "NC")))
summary(manova(cbind(low, high) ~ treatment, data=students, 
        subset = treatment %in% c("NC", "AA")))

Rによる統計解析ハンドブック 第4章演習問題解答

This post is solutions of chapter 4 of `A Handbook of Statistical Analyses Using R, Second Edition'.

このポストはRによる統計解析ハンドブック(第2版)の第4章の解答になります。

1

Test based on chi-square can't reject independency of region and site of lesion, but Fisher's exact test rejects that.

カイ2乗値に基づく検定では地域と病変部位は独立でないと言えないが、Fisherの正確検定では棄却される。

data("orallesions", package="HSAUR2")
chisq.test(orallesions) # p-value = 0.1096, accept
fisher.test(orallesions) # p-value = 0.01904, reject

2

labels in graph warp each other, so should I change row names of the dataframe to avoid wapping in advance?

グラフのラベルが重なるのですが、予めrownames関数でデータフレームの行名を変えなければならないのでしょうか。

library("vcd")
data("orallesions", package="HSAUR2")
assoc(orallesions)
mosaic(orallesions)

3

Difference of p-value is derived from the that of mean. To check this, boxplot is helpful. By executing sum(p.diff>0) in R, it turns out that elements of p.diff is always larger than 0.

違いは平均値の差の大きさに由来する。これを確認するためには箱ひげ図が便利である。Rでsum(p.diff>0)とすることで、すべて0より大きな値をとっていることが分かる。

library("coin")

n.test <- 100 # number of test
p.diff <- matrix(nrow=n.test, ncol=4) # matrix to store the difference of p-value
colnames(p.diff) <- c("diff=0.5", "diff=1", "diff=1.5", "diff=2")
sz <- 12 # 15 is upper limit on my computer. Be sure to make sz less than 15!
for( i in 1:n.test ){
  for( j in 1:4 ){
    g1 <- rnorm(sz, mean=0, sd=1)
    g2 <- rnorm(sz, mean=j/2, sd=1)
    lev <- c(rep(1, sz), rep(2, sz))
    ddd <- c(g1, g2)
    gs <- as.data.frame(cbind(ddd, lev))
    colnames(gs) <c("value", "gr")
    gs$gr <- factor(gs$gr)
    p.diff[i, j] <pvalue(wilcox_test(value ~ gr, data=gs, distribution=exact()))
               - pvalue(independence_test(value ~ gr, data=gs, distribution=exact()))
  }
}

boxplot(p.diff)

Rのpairs関数の出力をカスタマイズする

pairs関数(対散布図とも言う)はデフォルトの設定では黒い点でプロットするだけなので、グラフにより情報を持たせるため、グラフをカスタマイズする。対散布図は行方向、列方向に同じ数の図が描かれるが、i行j列の位置にある図とj行i列の位置にあるグラフは実質同じなので、工夫したグラフを試してみようということである。

ここではおなじみのirisデータを基に、右上のグラフには相関係数の絶対値、左下のグラフにはデータを色で分けたグラフと、直線性が強い組についてのみ回帰直線を描いてみることにする。

右上、左下のグラフを制御するのはpairs関数の引数のうち、panel、upper.panel、lower.panelである。デフォルトではpanelは散布図を描くようになっており、upper/lower.panelはそれを使うようになっているので、点による散布図が書かれることになる。upper/lower.panelの引数に与える関数を書き換えることで、グラフに手を加えることができる。

upper <- function(x, y, ...){ # 右上の散布図を設定
  oldpar <- par(usr = c(0, 1, 0, 1)) # 3行下でテキストの位置を(0.5, 0.5)にするための調整
  v <- abs(cor(x, y))
  sz <- 1 + v * 2 # テキストのサイズを直線性の強さで調整
  text(0.5, 0.5, sprintf("%.3f", v), cex = sz)
  par(oldpar)
}

lower <- function(x, y, ...){ # 左下の散布図を設定
  points(x, y, pch = 21, bg = c("red", "green3", "blue")[unclass(iris$Species)])
  if( abs(cor(x, y)) > 0.8 ){ # 直線性の強いデータの組だけ回帰式を作成
    lmobj <- lm(y~x)
    abline(lmobj)
  }
}

pairs(iris[,1:4], upper.panel=upper, lower.panel=lower)

気をつける必要があるのは、lower関数の1行目で、ここでplotを用いて書こうとするとエラーになる点である。低水準作図関数を使うようにすること。また、lower関数内で用いたpoints関数のbgパラメータのように、データをそのまますべて渡せば、すべての散布図で自動的にうまく描いてくれる。1つ1つの散布図を気にしてパラメータを変更する必要はない。

上のコードを実行すると、次のような画像が得られる。

Rによる統計解析ハンドブック 第3章演習問題解答

This post is solutions of chapter 3 of `A Handbook of Statistical Analyses Using R, Second Edition'.

このポストはRによる統計解析ハンドブック(第2版)の第3章の解答になります。

1

To decide which is superior, I chose wilcox.test as relative merit is decided by the absolute difference from 43 meters. t.test doesn't tell us which index is better estimator. conf.int in wilcox.test is needed to get estimated value which indicates the difference of rank. In below code execution printed value(=-3.1999...) shows that feet is superior than meter though null hypothesis couldn't rejected as p-value is slightly greater than 0.05.

どちらが優れているかを決めるため、wilcox.testを選んだ。というのも、優劣は43からの絶対誤差で決まるからである。t.testはどちらの指標が優れた推定量であるかを示してはくれない。wilcox.testのconf.intは順位の違いを示す推定値を得るために必要である。p値が0.05よりわずかに大きいため帰無仮説(2群は順位差がない)は棄却できないものの、実行時に表示される値はfeetがmeterよりも優れていることを示している。

data("roomwidth", package="HSAUR2")
convert <- ifelse(roomwidth$unit == "feet", 1, 3.28)
wt <- wilcox.test(I(abs(width*convert-43))~unit, data=roomwidth, conf.int=TRUE)
print(wt$estimate) # =-3.200

2

data("water", package="HSAUR2")
by(water, water$location, function(x){ return(cor(x$mortality, x$hardness))})

3

There is a mistake in the textbook written in Japanese. The numerator of formula in the question should be (njk-Ejk)/sqrt(Ejk). Otherwise it can't take a nagative value.

日本語版テキストの公式は間違いで、分子の方は(njk-Ejk)/sqrt(Ejk)が正しい。そうしないと(N(0, 1)の分布なのに)負の値がとれないよね～

data("pistonrings", package="HSAUR2")

calc.res <- function(x){
  sta.res <- chisq.test(x)$residuals
  x.exp <- outer(rowSums(x), colSums(x))/sum(x) # creates Ejk
  rr <- rowSums(x)
  cc <- colSums(x)
  nn <- sum(x)
  adj.res <- matrix(nrow=nrow(x), ncol=ncol(x))
  for( j in 1:nrow(x) ){
    for( k in 1:ncol(x) ){
      adj.res[j,k] <- (x[j,k]-x.exp[j,k])/sqrt(x.exp[j,k])/sqrt((1-rr[j]/nn)*(1-cc[k]/nn))
    }
  }
  adj.res <- as.data.frame(adj.res)
  rownames(adj.res) <- rownames(x)
  colnames(adj.res) <- colnames(x)
  return(list("sta.res"=sta.res, "adj.res"=adj.res))
}
pis.res <- calc.res(pistonrings)

4

The result of prop.test indicates that null hypothesis couldn't be rejected because p-value is greater than 0.05.

prop.testの結果はp値が0.05より大きいことから、(比率が一定という)帰無仮説が棄却できないことを示している。

data("rearrests", package="HSAUR2")
prop.test(rearrests, correct=FALSE)

Rによる統計解析ハンドブック 第2章演習問題解答

This post is solutions of chapter 2 of `A Handbook of Statistical Analyses Using R, Second Edition'. Note that my interpretations of data may differ from yours, but they are not the only solutions.

このポストはRによる統計解析ハンドブック(第2版)の第2章の解答になります。私のデータの解釈はあなたのそれらとは異なるかも知れませんが、唯一の解答ではないことに注意してください。

1

I aggregated expenditure for each gender. Compared to men, women tend to spend much money on housing and less on food.

性別ごとに支出金額を集計した。男性と比較して、女性は住まいへの費用が多く、食費への支出が少ない傾向にある。

lst <- split(household, household$gender)
expenses <- sapply(lst, function(x){colSums(x[,1:4])})
normalize <- matrix(rep(colSums(expenses),4), ncol=2, byrow=T)
ratios <- expenses/normalize
xpos <- barplot(ratios)
ypos <- diff(rbind( c(0, 0), apply(ratios, 2, cumsum)))/2 + 
        rbind( c(0, 0), apply(ratios, 2, cumsum))[1:4,] # calculate y coord to insert strings
text(rep(xpos, each=4), ypos, rep(rownames(ratios), 2))

2

As men get older, suicide rate increases. All of outliers belong to Hungary. Assume that library(HSAUR2) is already done. If not, execute data("suicides2", package="HSUAR2").

年齢が上がるに連れて自殺率が上がる。外れ値はいずれもハンガリーになっている。library(HSAUR2)は実行済みと仮定してください。もししていなければ、data("suicides2", package="HSUAR2")を実行するように。

boxplot(suicides2)

3

This answer shows how to change points into strings when using function pairs. By controlling parameter `panel', we can control both lower.panel and upper.panel. In case of changing counterpart, we should define functions for each panel. q2 required me to plot Life.Expectancy - Homicide relationship conditioned by Income, but no indication about it in text book. So I divided data whether it is higher or lower than the median income.

対散布図で描画するときに、点を文字に変えるにはどうしたら良いかという解答がこれになります。panelを制御することで、lower.panelとupper.panelを両方とも制御することができます。一方だけ変えたい場合は、それぞれに対して関数を定義してやる必要があります。q2では、平均収入を条件とした平均余命と殺人率のプロットを描くのですが、平均収入を条件とするということに関する指示がありません。ここでは、収入が中央値より高いか低いかでデータを分けました。

# states.abb <- abbreviate(rownames(USstates), minlength=2)
# q1
states.abb <- c("AL", "CA", "IA", "MS", "NH", "OH", "OR", "PA", "SD", "VT")
pairs(USstates, panel=function(x, y, ...){ text(x, y, states.abb) })
# q2
cond <- ifelse(USstates$Income < median(USstates$Income), "Lower", "Upper")
xyplot(Life.Expectancy ~ Homicide | cond, data=USstates)

4

When I used pairs, I found counterfeits were slightly short in diagonal length. To confirm that, I used boxplot.

対散布図を使ったところ、偽物がわずかに対角線の長さが短いことに気づきました。それを確認するために、箱ひげ図を利用してみました。

data("banknote", package="alr3")
pairs(banknote, col=banknote$Y+1) # black=genine, red=fake
boxplot(banknote$Diagonal~banknote$Y)

Rによる統計解析ハンドブック 第1章演習問題解答

This post is solutions of chapter 1 of `A Handbook of Statistical Analyses Using R, Second Edition'. As far as I know, there's no solutions on Web, so I decided to post my answers and hope to sophiscate them through discussing if needed.

このポストはRによる統計解析ハンドブック(第2版)の第1章の解答になります。私の知る限り、解法はWeb上にありませんので、私の解答を掲載することを決心しました。必要があれば議論して解答をより良くしたいと思っています。

Because almost no contents in the text book is written in this post, so please get one previously.

テキストに書かれている内容はほぼこのポストにありませんので、本は予め入手しておいてください。

1

library(HSAUR2)
prof.meds <- tapply(Forbes2000$profits, Forbes2000$country, median, na.rm=T)
print(prof.meds[c("United States", "United Kingdom", "France", "Germany")])

In the following answers, assume that `library(HSAUR2)' is already executed.

以後library(HSAUR2)は既に実行されていることを仮定せよ。

2

ger.deficit <- Forbes2000$country == "Germany" & Forbes2000$profits < 0
print(Forbes2000[ger.deficit, ])

3

Insurance Company is majority of this data (bermuda).

bermuda <- Forbes2000$country == "Bermuda"
print(table(Forbes2000[bermuda, "category"]))

4

big.profit <- Forbes2000[Forbes2000$profits >= sort(Forbes2000$profits, dec=T)[50],]
big.profit <- big.profit[!is.na(big.profit$profits),]
plot(big.profit$assets, big.profit$sales, log="xy", type="n")
text(big.profit$assets, big.profit$sales, abbreviate(big.profit$country))

5

print(tapply(Forbes2000$profits, Forbes2000$country, mean, na.rm=T))
larger <- function(x, th){
  return(length(which(x > th)))
}
print(tapply(Forbes2000$profits, Forbes2000$country, larger, 5))

Rで線形判別分析

フィッシャーの線形判別分析(LDA:Linear Discriminant Analysis)と呼ばれる手法についての簡単な紹介。

2群(実際には3つ以上のときもあるが)の線形判別分析というのは、あるデータに対し、それが2群のどちらに属するかを正負により判定する線形式で与えるというものである。正なら群1、負なら群2というように分けられる。今回は、Rのirisデータにある「がく幅、がく長」の2つが与えられたときに、それがSetosa、Versicolorのどちらに属するか判定するという、一番基本的な線形判別のコードを試してみた。

アルゴリズムに関しては、Wikipedia(信頼しすぎてはいけないが)の判別分析の項目に見ることができる。さて、RのMASSパッケージには線形判別のための関数が用意されているので、これを利用することにする。

library(MASS)

set.seed(1)
iris2 <- iris[1:100,c(1,2,5)] # がく長、がく幅列とsetosaとversicolorが含まれる行のみを使う
# 学習用データと検証用データにデータを分ける
mdl.idx <- sample(1:100, 50, rep=F)
# iris2の種類の水準にvirginicaが残っているため警告がでる
mdl <- lda(Species~., data=iris2[mdl.idx,]) # 線形判別モデル
# モデルに検証用データを突っ込む
result <- predict(mdl, iris2[-mdl.idx,])

xr <- range(iris2[,1]) # 描画範囲の設定
yr <- range(iris2[,2])
# 学習用データの描画
plot(iris2[mdl.idx,1], iris2[mdl.idx,2], col=iris2[mdl.idx,3], 
     xlim=xr, ylim=yr, xlab=names(iris)[1], ylab=names(iris)[2])
par(new=T)
# 検証用データの描画
plot(iris2[-mdl.idx,1], iris2[-mdl.idx,2], col=iris2[-mdl.idx,3],
     xlim=xr, ylim=yr, xlab="", ylab="", pch=2)
title("irisデータの線形判別")
# グラフに凡例を追加
leg <- c("train:setosa", "train:versicolor", "test:setosa", "test:versicolor")
legend(6.4, 4.3, leg, pch=c(1,1,2,2), col=c(1,2,1,2))
# 線形判別直線の描画範囲設定
xc <- seq(min(xr), 6.3, length=100)
# 判別直線の定数項は計算が必要(applyの結果が定数に相当する)
yc <- (apply(mdl$means %*% mdl$scaling, 2, mean)/mdl$scaling[2] 
       - mdl$scaling[1]/mdl$scaling[2] * xc )
points(x=xc, y=yc, type="l", col="green", lwd=2) # 判別の境界線描画
par(new=F)

注意しなければいけないのは、線形判別直線を引くところである。単純に実行すると、

y = -apply(mdl$means %*% mdl$scaling, 2, mean) + mdl$scaling[1]*がく長 + mdl$scaling[2]*がく幅

という判別式が作成されるので、がく幅=の形に変形してxc(がく長)からyc(がく幅)を計算している。変形する際は、境界を考えているので、y=0になることにも注意。

上のコードを実行すると次のような結果が得られる。ちょっと際どいデータもあるが、概ね完全な分離ができているといっても良いであろう。

R言語のsubset関数の調査

とあるR使いとコードについて話をしたとき、subset関数について話が出た。データ操作でよく使われそうなものなので、どういう挙動か見てみることにする。

subset関数はベクトル、マトリックス、データフレームから条件を満たす行のみを抽出する。

使用方法は次のようになっている。

• subset(x, ...)
• subset(x, subset, ...)
• subset(x, subset, select, drop = FALSE, ...)(マトリックスやデータフレームに対して)

引数は次のような意味をもつ。

• xは部分集合を求めたいオブジェクト。ここではデータフレームなどのデータと思って差支えない。
• subsetは残しておきたい要素、行の論理値表現。欠損値はfalseとみなされる。
• selectはデータフレーム、マトリックスから抜き出したい列を選ぶ表現を書く。
• dropは添え字オペレータ[で扱えるものにする。例えば1列のマトリックスはtrueに設定すればベクトル化されることになる。
• ...にはさらに渡す引数や他のメソッドからの引数を与える。

subsetは総称的関数であり、マトリックス、データフレーム、ベクトルに対するメソッドが与えられている(リストも含む)。パッケージやユーザはさらにメソッドを追加することもできる。通常のベクトルであれば単にx[subset & !is.na(subset)]とすればよい。データフレームに対しては、引数subsetは行に対して働く。subsetはデータフレームで評価され、列は変数として(名前で)言及されることもある。引数selectはデータフレームとマトリックスに対するメソッドにおいてのみ存在する。最初に列名を対応する列番号に置き換え、得られた整数ベクトルを索引づけるために用いる。これにより、標準的な索引づけの慣習は、例えば列の範囲を簡単に特定するあるいは1列をドロップするといったことができるようにしている。引数dropはマトリックスとデータフレームに対して索引づけのメソッドに移される: マトリックスに対する既定の動作はそのような添字をつけるのとは異なる。因子は部分集合化した後には空の水準があるかもしれない: 使われていない水準は自動的に削除される。

返り値は、ベクトルであれば選ばれた要素を、データフレームやマトリックスでは選ばれた行と列を引数のxに似たオブジェクトで返す。

この関数はインタラクティブな使用を意図した便利な関数である。プログラミングでは"["のような標準的な部分集合化する関数を利用するのがよい。特に普通ではない評価をするような部分集合を求める引数では、思わぬ結果を引き起こすことになるかもしれない。

Rで実験計画法 後編

4ヶ月ほど前に開かれたTokyo.R #22にて発表に利用した、実験計画法の話題に関する後編のスライドです。ソフトウェアのテストと絡めて、直交表の使い方などについて話をしてきました。前編についてはこちらからご覧ください。

Rで実験計画法 後編 from itoyan110

一様分布から標準正規乱数を発生させる

標準正規分布に従う乱数を発生させる簡単な方法として、一様乱数をいくつか発生させてその和をとり、期待値で引き算するという方法がある。そこで、得られた値が本当に標準正規分布といって差し支えないか調べてみることにした。以下のコードでは、6つの一様乱数(取りうる値の範囲は[0,1])とし、それらの和の期待値3を引いたものを大量に発生させ、値の分布が一様乱数であるか検定したものである。6としているのは、とりうる値を[-3, 3]とするためである。標準正規分布であれば、多くの場合、そこから発生する乱数は[-3, 3]の範囲に収まるためである。

# 一様乱数から正規乱数を生成するプログラム

set.seed(1)
pts <- 6 # 一つの標準正規乱数を発生させるのに用いる一様乱数の数
N <- 5000 # 発生させる正規乱数
num.data <- N*pts
rn <- runif(num.data)
rnd.mat <- matrix(rn,nrow=N) # 一様乱数の行列
num <- rowSums(rnd.mat)-pts/2

hist(num,breaks=20,xlim=c(-4,4),xlab="",ylab="freq",main="")
oldpar <- par(new=T)
curve(dnorm,-4,4,xlab="",ylab="",axes=FALSE)
par <- oldpar

# コルモゴロフ・スミルノフ(KS)検定
# 標準正規分布 N(0, 1)とデータを比較している
ks.test(num, "pnorm", mean=0, sd=1)

コルモゴロフ・スミルノフ検定は、データがある分布に従っているかを検定するというものである。上のコードは1標本KS検定に相当する。帰無仮説はデータは分布から発生したものではない、対立仮説はデータは分布から発生したものとなる。統計量は従うと想定される分布の累積分布関数と、データのそれを比較した結果から得られる差の値に基づく。統計量が小さくなると、帰無仮説が棄却されるということになる。結論は、KS検定ではp値はほぼ0となり、帰無仮説は標準正規分布ではないということになった。ただし、中心極限定理からデータ数を6から増やせば、標準正規分布に従うことが想定される。そこで、pts<-12としてみると標準正規乱数といえるようなp値となった（有意水準5%で帰無仮説を棄却できる）。必要なptsのサイズは標準化する分布(今回は一様分布)の対称の程度に依存していることに注意。今回利用した一様分布U(0,1)は0.5で対称なので、比較的少ないデータ数で標準正規分布が生成できるのである。

RのapplyとrowSumsの比較

各所で行方向の和を求める際にはapplyよりもrowSumsがオススメであるということを聞いた。本当にそうなのか、確認するためにコードを書いてみた。下記の2パターンでapplyとrowSumsの実行時間がどのように変化するかということを調べている。サイズによる処理時間の変化を観察することにする。

• 列数10で固定で、行数のみ増やす
  
• 行数10で固定で、列数のみ増やす
  

# 行方向のみ数が変化する(列数は10で固定)
for( sz in (1:10)*50000 ){
 data <- matrix(runif(sz*10),nrow=sz)
 print(system.time(rowSums(data)))
 # print(system.time(apply(data,1,sum))) # applyの実行時間を測る場合はコメントを解除
}

# 列方向のみ数が変化する(行数は10で固定)
for( sz in (1:10)*50000 ){
 data <- matrix(runif(sz*10),ncol=sz)
 print(system.time(rowSums(data)))
 # print(system.time(apply(data,1,sum)))
}

上のような時間を出力するコードを私のマシンで実行し、出力された経過時間(system.time関数の返り値の第3引数)をプロットしたのが下の二つの図である。上図は列数固定の結果である。横軸は行数で縦軸は経過時間(単位は秒)を表している。下図は行数固定とした結果である。横軸は列数で、縦軸は経過時間(単位は秒)を表している。経過時間は、表示など処理のすべてにかかった時間である。下図は時間が0になっているが、表示の都合上小数点2桁までしか表されず、それ以下の細かい値まで取れていないためである。各計測点で2回経過時間を計測した平均をとっているが、2回の測定とも、概ね同じような実行時間になっていた。この結果より次のことが言える。

• applyは行数に比例して処理時間がかかるのに対し、rowSumsは数十万行の規模でもほぼ一定の短い時間で処理が可能
  
• 行数が一定の場合でもrowSumsの性能がapplyを上回るものの、その差は列数固定の場合よりは小さい
  

なお、rowSumsの仲間で、列方向の和を求める関数colSumsが存在する。また、マトリックスの行/列方向の平均を求める関数として、rowMeans/colMeansというものが存在する。今回のような単純な集計に対しては、applyよりも更に速い処理ができる関数があるということは覚えておいて損はないと思う。

Rで実験計画法 前編

Tokyo.R#16にて、実験計画法の話題について話をしてきました。ソフトウェアテストに使われる直交表について学ぶついでに、本や資料をいくつか読んだので、それらをまとめました。

Rで実験計画法 前編

View more presentations from itoyan110

Rのmissing関数に関する調査

あるパッケージの関数を眺めていたらmissingという関数が出てきて、気になったので調べてみた。missing関数の呼び方は極めて単純である。

missing(x)

missingは関数の引数として欠けているか否かを評価する。missing関数は関数の中で評価される前に利用しないと意味がない。つまり、missing(x)とする前に、引数xに一切手を加えてはならないのである。また、ネストするとmissing関数は使えなくなる。同じ階層のみで利用可能である(将来的には変更されるらしいが)。

plot関数はx,y座標のyを省くと、勝手にxがyになり1,2,3,...をx座標に設定するが(plot(c(1,3,5))を試してみよ)、それはmissing関数を用いて実装することができる。Exampleの関数myplotがその例である。

myplot <- function(x,y) {
 if(missing(y)) { # yがなかったら
  y <- x # xをyに置き換えて
  x <- 1:length(y) # xを再設定 
 }
 plot(x,y)
}

Rのtable関数に関する調査

R言語のtable関数はクロス集計のために用いられる関数である。Excelではピボットテーブルに相当する機能である。呼び出しの形式は以下の通りである。

table(..., exclude = if (useNA == "no") c(NA, NaN), useNA = c("no","ifany", "always"), 
       dnn = list.names(...), deparse.level = 1) 

多くの場合、オプションは気にしなくても良いようである。単純な使い方から順に説明する。

1つのカテゴリ変数をtable関数に与えると、カテゴリごとの集計が行われる。次の例を見て欲しい。

> table(state.region) # アメリカ50州の国勢調査対象地域の分類
# 地域ごとの集計結果が出力される
state.region
    Northeast         South North Central          West 
            9            16            12            13 
> table(rbinom(10,1,0.5)) # 0,1が0.5の割合で現れる2項分布から10個データをランダムに生成
# 0,1になった個数が、それぞれ3、7であったことを示している(ランダムなので結果は毎回異なる)
0 1 
3 7 

2変数を与えるとクロス集計もできる。

# state.x77はアメリカ50州のデータ。平均所得の範囲を確認
> range(state.x77[,"Income"]) 
[1] 3098 6315
# 1000刻みで数値をカテゴリに変換する
> state.income <- cut(state.x77[,2],breaks=seq(from=3000,to=7000,by=1000))
> table(state.income,state.region)
               state.region
state.income    Northeast South North Central West
  (3e+03,4e+03]         2    10             0    1
  (4e+03,5e+03]         5     5            10    9
  (5e+03,6e+03]         2     1             2    2
  (6e+03,7e+03]         0     0             0    1

この表から、南部の収入がやや低い傾向を窺うことができる。余談ではあるが、1人当たりの平均所得が一番高い州(上の表では右下の1のこと)はアラスカであった。資源が豊富で、油田による収益が大きいのがその理由である。