R関数について

現在R、というより、統計の基礎を1から勉強中です。Rを使いこなすためにはある程度の統計的知識が必要であるのですが、考え方を変えて、Rを使いこなすために統計知識を勉強する、Rを勉強することで統計をも勉強すると考えて、1から勉強することにしました。

本日勉強した関数

table( ):

度数分布表を作成する

length( ):

データの個数をカウントする。データ数を求めるときに便利

summary():

データの最小値、中央値、平均値、最大値、第一四分位数、第三四位数を算出する。

var ( ) :

データの分散(不偏分散)を算出

sd( ):

データの標準偏差(不偏分散の平方根)を算出

目標はRで検定の関数を書いて、いつでも使えるようにしておくことですが、目標はまだまだ遠いです。

Rによる効果量の重み付き平均算出の関数

丹後(2002)の標準化された平均値の差(hedge’s gと同じ式)のメタ分析
母数モデル(fixed models)に基づいて作成(pp.83-84)

ES<-function(Ne,Me,SDe,Nc,Mc,SDc){

# サンプルサイズの合計

N<-Ne+Nc

# プールされた標準偏差

SD<-sqrt(((Ne-1)*(SDe^2)+(Nc-1)*(SDc^2))/(N-2))

# Hedge’s gを算出

g<-(Me-Mc)/SD

# 各研究の標準誤差

SE<-sqrt((N/(Ne*Nc))+((g^2)/(2*N)))

# 各研究におけるgの信頼区間(上限)

CIH<-g+1.96*SE

# 各研究におけるgの信頼区間(下限)

CIL<-g-1.96*SE

# 各研究の重み

w<-1/(SE^2)

#結果をデータフレーム化

DF<-data.frame(Ne,Me,SDe,Nc,Mc,SDc,g,w,CIH,CIL)
}

使用方法

Ne:実験群のサンプル
Me:実験群の平均値
SDe:実験群の標準偏差
Nc:統制群のサンプル
Mc:統制群の平均値
SDc:統制群の標準偏差
ans<-ES(Ne,Me,SDe,Nc,Mc,SDc)
ans

ansと入力するとデータフレームが表示される(はず)

効果量の重み付き平均(gwとする)と95%信頼区間の算出

# 文字wにESのデータフレームのw行の値を格納

w<-ans$w

# 文字gにESのデータフレームのg行の値を格納

g<-ans$g

# Hedge’s gの重みつき平均を算出

gw<-sum(w*g)/sum(w)

# g(重み付き平均)の信頼区間(上限)

CI.high<-gw+(1.96/sqrt(sum(w)))

# g(重み付き平均)の信頼区間(下限)

CI.low<-gw-(1.96/sqrt(sum(w)))

均質性の検定

# 均質性の検定

Q <-sum(w*((g-gw)^2) )

で算出されるはず。

備忘録を兼ねて記録しています。

参考文献

丹後俊郎(2002).『メタ・アナリシス入門ーエビデンスの統合をめざす統計手法』東京:朝倉書店

2011 5/18追記

標準誤差(SE)の計算式に誤りがありました。修正済み。

Rによるデータフレームの読み込みとベクトルへの変換

外部ファイル(csvファイル)を統計ソフトRに読み込み、特定の列をベクトルに変換する方法。R mac版の場合。備忘録を兼ねて。

① 作業ディレクトリを設定する。

R→環境設定→起動と進むと、「ディレクトリ」を設定する場面が出てくるので、任意の場所に決める。指定したディレクトリの中にcsvファイルを保存する。ディレクトリを設定しないとファイルを読み込むことができない(と思われる)ので、注意。

② csvファイルを読み込む:

csvファイルの1行目が項目名の場合はheader=TRUEを記述する

x<-read.csv(“sample.csv”, header=TRUE)      # xという文字に、sample.csvを読み込む

csvファイルの1行目から数値が入力されている場合はheader=Falseを記述する

x<-read.csv(“sample.csv”, header=F)

例:sample.csv(以下参照)を読み込む場合

———sample.csv———————-

participants  pre-test  post-test

1                  20            40

1                  18             38

1                  16             55

1                  14             40


1列目が項目名になっているため、header=TRUEを記述する。

> x         # xを読み込む
participants pre.test post.test
1            1       20        40
2            2       18        38
3            3       16        55
4            4       14        40

と表示されます。

③ 特定の行の数値をベクトルに変換する

上記のデータフレームからpre.testの数値だけを表示させる時には

x$pre.test   # xの中のpre.testの行を表示させる

と入力。すると

> x$pre.test
[1] 20 18 16 14

と返してきます。このプログラムと、as.vector()関数(( )内のものをベクトルに変える関数)を用いて、以下のように記述

> pre<-x$pre.test     # preという文字の中にsample.csvファイルのpre.testの行の数値を格納
> pre            # preと入力すると、以下の[1]のように返してくる。
[1] 20 18 16 14
> as.vector(pre)      # preに格納されている数値をベクトルに変換
[1] 20 18 16 14
> is.vector(pre)      # preがベクトルであるかを調べる関数
[1] TRUE

以上で終了。

この方法を知っているだけで、エクセルファイルでまとめたデータをRで入力して分析するのにとても楽になるような気がします。まだ試していませんが、c( )の中に一つ一つデータを手打ちで入れていくよりはよほど簡単なのは明白です。