0. はじめに
もっとも基本的な統計量である平均値に加え、標準偏差、標準誤差の計算をRで行っていきましょう。
だんだん統計らしくなり、「実際に使っている感」が増すと思いますので、ぜひ一緒にコードを打ってもらえたらと思います！
今回も前回までにやったことがいっぱい登場しますが、各所に前回までのリンクを載せてありますので、わからない箇所は必要に応じて振り返ってもらえたらと思います。
【前回の記事はこちら】
0-1. この記事の到達目標
- 前回までの列の取り出しや抽出の復習をする。
- 平均値、標準偏差、標準誤差を計算できる。
- どんな処理過程を踏んでいるのか理解を進める。
0-2. 初めて登場する関数・演算子
mean(ベクトル)- ベクトルの平均値を計算します。
sd(ベクトル)- ベクトルの標準偏差を計算します。
sqrt(数値)- 数値の平方根を返します。
function(引数){処理内容}- 自作関数を定義します。
std.error(ベクトル)- ベクトルの標準誤差を計算します。
plotrixパッケージ内に含まれます。
1. 前準備
前回、前々回も用いたトウモロコシのデータを今回も練習として用います。
以下のコードをコピペ実行で、このトウモロコシのデータdfが取得できますので是非いっしょにやってみましょう。
ID<-c(1:5,1:5)
Treatment<-c(rep("Fertilizer",5),rep("Control",5))
hight<-c(180,178,192,182,185,165,171,173,168,169)
yield<-c(410,375,368,391,389,352,345,364,378,342)
df<-data.frame(ID,Treatment,hight,yield)
dfのなかみを確認してみると、以下のような感じになっていると思います。
> df
ID Treatment hight yield
1 1 Fertilizer 180 410
2 2 Fertilizer 178 375
3 3 Fertilizer 192 368
4 4 Fertilizer 182 391
5 5 Fertilizer 185 389
6 1 Control 165 352
7 2 Control 171 345
8 3 Control 173 364
9 4 Control 168 378
10 5 Control 169 342
もし、前回までの復習としてExcelからデータを読み込んでみたい方は、
ワーキングディレクトリを設定し、ファイルを読み込み、要因の列をfactorに変換しおわったところから始めましょう。
※トウモロコシのデータの読み込みについては第4回へ
※要因の列のfactorに変換することについては第5回へ
一例としてこんな感じになるんじゃないかと思います。
#ワーキングディレクトリの設定
> setwd("C:/Users/Surku/Desktop/R")
#データの読み込み
> df<-read.table("トウモロコシ.csv",header = TRUE,sep=",")
> df[,2]<-as.factor(df[,2])
2. 平均値,標準偏差,標準誤差の計算
それでは、早速平均値、標準偏差、標準誤差の計算をやっていきましょう。
今回は、Fertilizer区とControl区のグループで計算していきます。
前回のsubset()関数がたくさん出てきますので、忘れてしまった方は必要に応じて前回をご覧ください。
2-1. mean()で平均値の計算
平均値を計算するにはmean()関数を用います。引数に平均したい数値が入ったベクトルを代入します。もちろん文字列や要因といったデータ型を指定するとエラーになります。
試しに、テキトーにベクトルを作って計算してみます。
> number<-c(1,2,3,4,5)
> mean(number)
[1] 3
> character<-c("いち","に","さん","よん","ご")
> mean(character)
[1] NA
警告メッセージ:
mean.default(character) で:
引数は数値でも論理値でもありません。NA 値を返します
ではトウモロコシのデータを使って、Fertilizer区の背丈の平均値を計算してみましょう。
Fertilizer区の平均値を計算するのに、Control区の値は不要なのでsubset()でFertilizer区のみ抽出を行い、新たにdf_Fertilizerというデータフレームを作ります。
次に、データフレームから背丈の列（3列目）だけ取り出します。
mean()の引数に入れることができるのはベクトルであり、データフレームを入れることはできません。取り出した列はベクトルのデータ型を持っているので、mean()に入れることができ、計算させることができます。
#TreatmentがFertilizerのみ抽出
> df_Fertilizer<-subset(df,df[,2]=="Fertilizer")
> df_Fertilizer
ID Treatment hight yield
1 1 Fertilizer 180 410
2 2 Fertilizer 178 375
3 3 Fertilizer 192 368
4 4 Fertilizer 182 391
5 5 Fertilizer 185 389
#データフレームのままmean()に入れるとエラー
> mean(df_Fertilizer)
[1] NA
警告メッセージ:
mean.default(df_Fertilizer) で:
引数は数値でも論理値でもありません。NA 値を返します
#背丈の列だけ取り出したものはベクトル
> df_Fertilizer[,3]
[1] 180 178 192 182 185
#ベクトルだと計算できる
> mean(df_Fertilizer[,3])
[1] 183.4
同じ要領で、Control区の収量(yield)の平均値も計算してみましょう。
> df_Control<-subset(df,df[,2]=="Control")
> mean(df_Control[,4])
[1] 356.2
2-2. sd()で標準偏差の計算
標準偏差と標準誤差の違いについては次回以降の記事にまとめるつもりで、ここではとりあえずコードの書き方だけおさえていただければと思います。
標準偏差の計算の仕方はとっても簡単で、前項でmean()となっていたところをsd()に置き換えるだけでOKです。sdとはstandard deviationの略です。
せっかくなので、Fertilizer区とControl区両方で、背丈と収量の標準偏差を求めてみましょう。
> df_Fertilizer<-subset(df,df[,2]=="Fertilizer")
#Fertilizer区の背丈
> sd(df_Fertilizer[,3])
[1] 5.458938
#Fertilizer区の収量
> sd(df_Fertilizer[,4])
[1] 16.22652
> df_Control<-subset(df,df[,2]=="Control")
#Control区の背丈
> sd(df_Control[,3])
[1] 3.03315
#Control区の収量
> sd(df_Control[,4])
[1] 14.83914
2-3. sd()から標準誤差を計算
標準誤差の計算はちょっと厄介です。
なぜなら標準誤差を一発で計算する関数がデフォルトで入っていないためです。
ここでは、デフォルトの関数を組みわせて標準誤差を計算してみましょう。
標準誤差の公式は次の通りです。
標準誤差= \frac{標準偏差}{\sqrt{サンプルサイズ}}
標準偏差の計算には先ほどのsd()関数が使えます
ここでいうサンプルサイズはあまり聞きなれないかもしれませんが、ここではデータの個数と考えてください。
データの個数を調べるためにはlength()関数を使うことができます。
この関数は第2回で登場しました。
平方根の計算にはsqrt()関数が使えます。
引数に数値を入れると平方根を返してくれます。試しにいくつか数字を入れてみましょう。
> sqrt(2)
[1] 1.414214
> sqrt(169)
[1] 13
これら組みわせて、標準誤差の公式を書き換えてみると
標準誤差＝\frac{sd()}{sqrt( length() )}
となります。
早速、sd()とlength()の引数に、Fertilizer区の背丈データを入れてみましょう。
#TreatmentがFertilizerのみ抽出
> df_Fertilizer<-subset(df,df[,2]=="Fertilizer")
#背丈の列だけ取り出した
> df_Fertilizer_height<-df_Fertilizer[,3]
#一応、ベクトルであることを確認
> df_Fertilizer_height
[1] 180 178 192 182 185
#標準誤差の公式に代入
> sd(df_Fertilizer_height)/sqrt(length(df_Fertilizer_height))
[1] 2.441311
長々とコードを書きましたが、もしこれを1行で書くならこんな感じです。
ごちゃごちゃしていますがどんな手順を踏んだのか復習になれば幸いです。
> sd(subset(df,df[,2]=="Fertilizer")[,3])/
+ sqrt(length(subset(df,df[,2]=="Fertilizer")[,3]))
[1] 2.441311
2-4. 簡単に標準誤差計算できないかな?
先ほどはsd()から標準誤差を計算しましたが、ちょっとややこしい部分がありましたよね。もっと簡単に標準誤差をできないか、その方法を@StrawBerryMoonさんより提示いただきましたので、こちらでも紹介します。
1つ目は、標準誤差を計算する関数を自分で定義する方法です。
標準誤差の公式からse()という関数を自分で定義し、引数にベクトルを指定しています。
関数の定義方法は第9回で扱いましたので、必要に応じてご確認ください。
> se<-function(x){
+ sd(x)/sqrt(length(x))
+ }
> df_Fertilizer<-subset(df,df[,2]=="Fertilizer")
> se(df_Fertilizer[,3])
[1] 2.441311
2つ目はplotrixパッケージ内のstd.error()関数を用いる方法です。
パッケージのインストール方法を読み込みについては次回で取り上げますのでご覧ください！
ちなみにplotrixパッケージは図表を描画するためによく使用されるパッケージになっています。
> install.packages("plotrix")
> library("plotrix")
> std.error(df_Fertilizer[,3])
[1] 2.441311
3．次回