Rメモ
Rメモ目次
researchmapのバージョンアップによるリンク切れを修正 2020307
参照に便利なように,目次を作りました
【データハンドリング】
- データハンドリングのよく使うプログラム集(自分用)
- StataやSASにある変数のkeepやdropをRで行うには
- EZRでパッケージ中のRデータセットをStataデータセットに変換する
- Rマークダウンが使えるようになるまで
- Rのデータセットをcsvファイルにエクスポートする
- パッケージltmのEnvironmentデータと因子型の説明
- 因子型変数の水準の順序(levels)を変える
- オフラインでRにパッケージをインストールする
- 連続量の変数を分割してカテゴリ変数に区分する
- 変数をloopで変えながら度数分布表を一度にたくさん作る
【解析】
- クロス表で欠損含む頻度や行及び列の割合を表示してカイ二乗検定も
- Rを使って平均値と標準偏差と人数だけから2群のt検定(おまけに効果量)
- 2群のt検定,効果量,95%信頼区間(データセットに変数としてある場合)
- 1要因分散分析のオムニバス効果量
- lavaanで確認的(確証的/検証的)因子分析:confirmatory factor analysis(CFA)
- lavaanで確認的因子分析:CFA(2)完全情報最尤(FIML)で推定
- 再検査(再テスト)信頼性のための級内相関係数(ICC)(3):Rで算出
- 測定誤差:測定の標準誤差(SEM:standard error of measurement)
- 項目反応理論(1):ltmパッケージで段階反応モデル
- ロジスティック回帰でオッズ比と95%信頼区間を出す 20180924 追記
-
【図・グラフ】
変数をloopで変えながら度数分布表を一度にたくさん作る
(20201127 出力を画像に変更)
データを把握したり可視化したりする際に,一度にたくさんの項目の度数分布を出したい時があります。そんなときに,指定した変数のすべてについて実行する方法があると便利です。もっといいやり方がありそうですが,とりあえず自分で見つけた方法について紹介します。
データとして,psychパッケージのbfiデータを使います。度数分布を出すには,table( )でもよいですが,tabyl( )の方が便利なので,janitorパッケージを使います。
janitorパッケージでtidyな度数分布表やクロス表を作成
library(psych)
library(janitor)
まず最初に,格納用の空のリストを作っておきます。resと名前をつけてます。
次に,変えていきたい変数一覧を準備します。ここではavarsと名前をつけました。
avars <- c("A1", "A2")
"A1","A2"のところは,変数名の位置で指定した
としても同じです。
res[[i]] <- tabyl(bfi[[i]])
}
res
これらの結果はリストresの中で別々の要素なので,最後に1つのデータフレームにまとめた方が便利です。そこで,どの質問項目の結果か分かるようにするために,項目名をそれぞれ値とした新たな変数を作ります。これは色々試した結果,以下のようにしないとできなかったので,とりあえずですが,次のようにします
res[[i]]$varname <- avars[i]
}
最後に,それぞれのリストの要素を1つのデータフレームにまとめます。
ロジスティック回帰でオッズ比と95%信頼区間を出す
加工せずにすぐに使えるデータとして,effectsパッケージのTitanicSurvivalデータセットを呼び出します。全体のデータの一部分のようなので,少なめでn=1,309となっています。
summary(TitanicSurvival)
survived sex age passengerClass
no :809 female:466 Min. : 0.1667 1st:323
yes:500 male :843 1st Qu.:21.0000 2nd:277
Median :28.0000 3rd:709
Mean :29.8811
3rd Qu.:39.0000
Max. :80.0000
NA's :263
説明変数を,船室等級,性別,年齢として,応答変数を生還の有無としたロジスティック回帰は以下のコードで実行できます。そしてsummary()で結果を見ることができます。
summary(titanic)
Call:
glm(formula = survived ~ passengerClass + sex + age, family = binomial,
data = TitanicSurvival)
Deviance Residuals:
Min 1Q Median 3Q Max
-2.6399 -0.6979 -0.4336 0.6688 2.3964
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.522074 0.326702 10.781 < 2e-16 ***
passengerClass2nd -1.280570 0.225538 -5.678 1.36e-08 ***
passengerClass3rd -2.289661 0.225802 -10.140 < 2e-16 ***
sexmale -2.497845 0.166037 -15.044 < 2e-16 ***
age -0.034393 0.006331 -5.433 5.56e-08 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1414.62 on 1045 degrees of freedom
Residual deviance: 982.45 on 1041 degrees of freedom
(263 observations deleted due to missingness)
AIC: 992.45
Number of Fisher Scoring iterations: 4
Coefficients:の部分に,各説明変数の回帰係数と標準誤差が出てきます。オッズ比と95%信頼区間については,exp()を使って別途計算しないといけません。coef()とconfint()でそれぞれ値を取り出して,exp()の中に入れます。
Intercept) passengerClass2nd passengerClass3rd sexmale age
33.85457044 0.27787894 0.10130084 0.08226211 0.96619149
exp(confint(titanic)) #オッズ比の95%信頼区間
2.5 % 97.5 %
(Intercept) 18.11476468 65.2744468
passengerClass2nd 0.17763418 0.4303922
passengerClass3rd 0.06453100 0.1565287
sexmale 0.05906841 0.1133112
age 0.95410521 0.9781055
上記のような面倒くさい処理をせずとも,epiDisplayパッケージのlogistic.display()関数を使うと,2変量のオッズ比と調整済みオッズ比と95%信頼区間を一発で出してくれるのですごく便利です。
Logistic regression predicting survived : yes vs no
crude OR(95%CI) adj. OR(95%CI) P(Wald's test) P(LR-test)
passengerClass: ref.=1st < 0.001
2nd 0.45 (0.32,0.63) 0.28 (0.18,0.43) < 0.001
3rd 0.2 (0.15,0.28) 0.1 (0.07,0.16) < 0.001
sex: male vs female 0.08 (0.06,0.11) 0.08 (0.06,0.11) < 0.001 < 0.001
age (cont. var.) 0.99 (0.98,1) 0.97 (0.95,0.98) < 0.001 < 0.001
Log-likelihood = -491.2266
No. of observations = 1046
AIC value = 992.4531
logistic.display(titanic, simplified = TRUE)
そして,調整済みオッズ比の結果の部分だけ欲しい時は,引数にsimplified = TRUEを足すと,すっきりした結果が出てきます。
logistic.display(titanic, simplified = TRUE)
OR lower95ci upper95ci Pr(>|Z|)
passengerClass2nd 0.27787894 0.17859824 0.4323486 1.364063e-08
passengerClass3rd 0.10130084 0.06507439 0.1576943 3.666157e-24
sexmale 0.08226211 0.05941134 0.1139017 3.782932e-51
age 0.96619149 0.95427655 0.9782552 5.556650e-08
1. As a factor: ‘success’ is interpreted as the factor not having the first level (and hence usually of having the second level).
$`levels`
[1] "no" "yes"
$class
[1] "factor"
broomパッケージを使う方法もありました。tidyさを追求するならこちらがよいかも。
tidy(titanic, conf.int = TRUE, exponentiate = TRUE)
term estimate std.error statistic p.value conf.low conf.high
1 (Intercept) 33.85457044 0.326702200 10.780687 4.247041e-27 18.11476468 65.2744468
2 passengerClass2nd 0.27787894 0.225538190 -5.677840 1.364063e-08 0.17763418 0.4303922
3 passengerClass3rd 0.10130084 0.225801920 -10.140129 3.666157e-24 0.06453100 0.1565287
4 sexmale 0.08226211 0.166036501 -15.043949 3.782932e-51 0.05906841 0.1133112
5 age 0.96619149 0.006331001 -5.432511 5.556650e-08 0.95410521 0.9781055
フォレストプロットをRのggplot2で描く
グラフにしたい場合にどうやるか説明します。
なお,ここでのグラフ作成は,Barnkob様のブログ記事を参考にしました。
職場環境要因と抑うつ症状の関連
Conditions | OR | OR_lower | OR_upper |
decision latitude | 0.73 | 0.68 | 0.77 |
job strain | 1.74 | 1.53 | 1.96 |
bullying | 2.82 | 2.21 | 3.59 |
のFig.2よりsummary odds ratioを抜粋
まずデータをRに取り込むために,表をすべて選択して,excelなどに貼り付けます。
それをコピーした後に
とするとデータが取り込みできます。または,以下のコードを実行しても同様です。
OR=c(0.73,1.74,2.82),
OR_lower=c(0.68,1.53,2.21),
OR_upper=c(0.77,1.96,3.59)
)
次にパッケージをロードしますが,以下のコードはパッケージが入っているか
どうか確認し,なければインストールもしてくれるので便利です。
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')
当然ggplot2は入ってるよ,という方は単に読みこむだけです。
以下のコードが,グラフ描画の本体です。
geom_linerange(size=8, colour="#a6d8f0") +
geom_hline(aes(x=0, yintercept=1), lty=1) +
geom_point(size=3, shape=21, fill="#008fd5", colour = "white", stroke = 1) +
scale_y_continuous(limits = c(0.5, 4)) + #Yの範囲(0.5~4まで)
coord_flip() + #縦横置き換え
ggtitle("Odds ratio for occupational factor") #タイトル
p
このとおり,綺麗な図が作成されます。
全体の見栄えを変えたい場合は以下のようにするとできます。
p + theme_classic() #背景,グリッド消す
今後,見栄えの変え方などいろいろ試していきたいと思います。
項目反応理論(1):ltmパッケージで段階反応モデル
以下のサイトで説明されている通りにコピペすればすぐに
できます。
ここでは,もっとかゆい所に手が届くような部分を
説明していきます。基礎知識は自分がまとめた下記スライド参照。
ここでは,パッケージltmを用いて,サンプルデータセットの
Environmentを使った段階反応モデルの基本的な実行と作図を
説明します。まずはltmパッケージを読み込みます。
データセットの説明は
でみられますが,概略を説明すると,
・1990年の英国社会的態度調査(Brook et al., 1991)の291名の回答データ
・変数
LeadPetrol ガソリンからの鉛
RiverSea 河川・海洋汚染
RadioWaste 放射性廃棄物の輸送と貯蔵
AirPollution 大気汚染
Chemicals 毒性化学物質の輸送と廃棄
Nuclear 原子力発電所のリスク
・回答選択肢
"very concerned", "slightly concerned" "not very concerned"
とても関心がある, 少し関心がある, あまり関心がない
です。それでは解析に入っていきます。
まずltmパッケージの便利関数descriptで,データの基本
情報を出力します。
以下,出力の抜粋です。
Sample:
6 items and 291 sample units; 0 missing values
Proportions for each level of response:
very concerned slightly concerned not very concerned
LeadPetrol 0.6151 0.3265 0.0584
RiverSea 0.8007 0.1753 0.0241
RadioWaste 0.7457 0.1924 0.0619
AirPollution 0.6495 0.3196 0.0309
Chemicals 0.7491 0.1924 0.0584
Nuclear 0.5155 0.3265 0.1581
Frequencies of total scores:
6 7 8 9 10 11 12 13 14 15 16 17 18
Freq 96 51 37 27 26 18 13 7 6 6 1 1 2
6つの項目があり,291名のデータがあり,欠損値が0という情報
がまず出てきます。次に,それぞれの項目の回答の割合が出ます。
たとえば,RadioWasteはとても関心を持っている者が約75%いる
ということになります。最後のセクションは,総合得点の度数分布
が出力されています。これらに加え,Cronbachのα係数も出力
されます。
いよいよ段階反応モデルの実行ですが,とても簡単で,
だけです。
Coefficients:
Extrmt1 Extrmt2 Dscrmn
LeadPetrol 0.487 2.584 1.378
RiverSea 1.058 2.499 2.341
RadioWaste 0.779 1.793 3.123
AirPollution 0.457 2.157 3.283
Chemicals 0.809 1.868 2.947
Nuclear 0.073 1.427 1.761
Log.Lik: -1090.404
と基本的な情報が出てきます。
次に,それぞれの項目の特性反応曲線(characteristic response curves [CRCs])
を描画します。
項目RadioWasteだけを見たいという場合は,3番目の変数なので,
もう少し見やすくしてみます。
凡例はlegend=Tで,左側に持ってくるためにcx="left"としました。
全項目のCRCを1つの図に出力したい場合は次のようにします。
plot(fit, legend=T, cx="left")
その他の主要な指標は以下の通りです。
テスト情報関数(total information function[TIF]) ※
plot(fit, type = "IIC", items = 0, plot = T)
項目情報関数(item information function[IIF]) (2016/10/15 引数を修正)
というわけで,たいていのことは簡単にできてしまいます。
自分のデータセットで試したい,と思った場合に,恐らくいちばん
引っかかるのが,変数の型がfactorになってないなどの点だと
思います。それは,以下の説明を参考にするとうまくいくかも
しれません。
※standard errorと一緒に描画したい場合
vals <- plot(fit, type = "IIC", items = 0, plot = T)
plot(vals[, "z"], 1 / sqrt(vals[, "test.info"]), type = "l", lwd = 2,
xlab = "Ability", ylab = "Standard Error",
main = "Standard Error of Measurement")
【履歴】
2016/10/15 項目情報関数の引数がテスト情報関数と同じだったのを修正
連続量の変数を分割してカテゴリ変数に区分する
例えば年齢から20代,30代といった年齢層に分けた変数を作るような時に便利
なのがcut関数です。
(1) 単純な例(0~10の値)
x
[1] 0 1 2 3 4 5 6 7 8 9 10
ここで,5以下と6以上で分けたいときは以下のようにします。
x2
[1] (-1,5] (-1,5] (-1,5] (-1,5] (-1,5] (-1,5] (5,10] (5,10] (5,10] (5,10] (5,10]
Levels: (-1,5] (5,10]
ちゃんと区分されたかもっとわかりやすく確認します。
x2
x (-1,5] (5,10]
0 1 0
1 1 0
2 1 0
3 1 0
4 1 0
5 1 0
6 0 1
7 0 1
8 0 1
9 0 1
10 0 1
なぜbreakesの最初に来る数字が-1なのかというと,ここに
書かれている数字は,それより大きい値から次のカテゴリ
に区分されるという意味だからです。
参考までに間違った例として,左端を0にしてみると,
x2X
[1] <NA> (0,5] (0,5] (0,5] (0,5] (0,5] (5,10] (5,10] (5,10] (5,10] (5,10]
Levels: (0,5] (5,10]
このように,0が含まれず欠損の<NA>になってしまいます。
次にx2に戻り,値のラベルについて自動で(-1,5],(5,10]というラベルが付きますが,
"(" の後に来る数字は入らない,"]"の前に来る数字まで入る,
と覚えるとよさそうです。ラベルを1,2にしたい場合は次のようにします。
x3
[1] 1 1 1 1 1 1 2 2 2 2 2
Levels: 1 2
※20160918修正。ラベルを0,1としていたが,ラベルが何であっても内部的には1,2
が当てられているので,それを反映するように修正。確認するには
ちなみに,引数のrightがデフォルトでTRUEになっているので,これを
FALSEにすると,次のようにbreaksを設定する必要があります。
x4
[1] [0,6) [0,6) [0,6) [0,6) [0,6) [0,6) [6,11) [6,11) [6,11) [6,11) [6,11)
Levels: [0,6) [6,11)
(2) 年齢から年齢層のカテゴリ変数作成
範囲を20歳以上45歳未満として年齢データを50件作成します。再現できる
ようにset.seedも。
age <- floor(runif(50, min = 20, max = 44.9)) #floorは関数は切り捨て,44を出すにはmaxを44.9に
sort(age)
[1] 20 21 21 22 23 23 23 23 25 25 25 26 26 27 27 27 28 29 30 30 30 31 31 31 31 33 33 33 34 34 35 36 36 37 37
[36] 37 38 39 39 41 41 42 42 42 42 43 43 43 43 44
念のため最小値,最大値を確認してみます
Min. 1st Qu. Median Mean 3rd Qu. Max.
20.00 26.25 32.00 32.44 38.75 44.00
それでは,20~29歳,30~39歳,40~44歳の3つの水準に区分してみます。
table(age,agec)
agec
age (19,29] (29,39] (39,44]
20 1 0 0
21 2 0 0
22 1 0 0
23 4 0 0
25 3 0 0
26 2 0 0
27 3 0 0
28 1 0 0
29 1 0 0
30 0 3 0
31 0 4 0
33 0 3 0
34 0 2 0
35 0 1 0
36 0 2 0
37 0 3 0
38 0 1 0
39 0 2 0
41 0 0 2
42 0 0 4
43 0 0 4
44 0 0 1
ちゃんと区分できました。
次の例のようにbreaksの左端と右端は,Infでおきかえてもよいです。
左端にはマイナス(-)が付きます。
ラベルを変更する場合は以下のように
labels=c("age20-29","age30-39","age40-44"))
table(age,agec3)
agec3
age age20-29 age30-39 age40-44
20 1 0 0
21 2 0 0
22 1 0 0
23 4 0 0
25 3 0 0
26 2 0 0
27 3 0 0
28 1 0 0
29 1 0 0
30 0 3 0
31 0 4 0
33 0 3 0
34 0 2 0
35 0 1 0
36 0 2 0
37 0 3 0
38 0 1 0
39 0 2 0
41 0 0 2
42 0 0 4
43 0 0 4
44 0 0 1
(2018/01/21 追記)
tidyverseを使う場合はこちら↓
連続量の変数を分割してカテゴリ変数に区分する(2)
オフラインでRにパッケージをインストールする
いかにしてその環境を整えるかについてのメモです。
↓のサイトの解説にすごくお世話になりました
miniCRANパッケージを使ってパッケージをオフラインインストールする
以下は,Windowsの場合での説明です。
(1)ネットにつながる別のパソコンでパッケージをダウンロードする
当然,パッケージを入手するには,ネットにつながる環境が必要です。
個々のパッケージは,CRANのページからダウンロードできますが,
ちゃんと使えるようにするためには,そのパッケージ単体では不足するケースが多く,
依存関係にあるパッケージも一緒に入手する必要があります。
以下の方法で,依存関係を調べます。上記リンク先の解説と同じく
ここでもltmパッケージを指定してみます。
pkgs<-pkgDep("ltm")
pkgs
[1] "ltm" "MASS" "msm" "polycor" "survival"
[6] "mvtnorm" "expm" "sfsmisc" "Matrix" "lattice"
これで,ltmを使う際には,他の9個のパッケージも同時にインストールしておいた
方がよさそうだと分かります。
さて,ここからが上記リンク先と異なり,WindowsのPCを想定した
手順となります(自分が理解するまで苦労した所らへんです)。
まず,pkgsオブジェクトに格納されているパッケージをすべてダウンロードします。
(例:C:/rpkgs)", type="win.binary")
です。実行すると,指定したフォルダにパッケージがダウンロードされ,
以下のような出力が出ます。
[,1] [,2]
[1,] "ltm" "C:/pkgs/ltm_1.0-0.zip"
[2,] "MASS" "C:/pkgs/MASS_7.3-45.zip"
[3,] "msm" "C:/pkgs/msm_1.6.1.zip"
[4,] "polycor" "C:/pkgs/polycor_0.7-8.zip"
[5,] "survival" "C:/pkgs/survival_2.39-5.zip"
[6,] "mvtnorm" "C:/pkgs/mvtnorm_1.0-5.zip"
[7,] "expm" "C:/pkgs/expm_0.999-0.zip"
[8,] "sfsmisc" "C:/pkgs/sfsmisc_1.1-0.zip"
[9,] "Matrix" "C:/pkgs/Matrix_1.2-6.zip"
[10,] "lattice" "C:/pkgs/lattice_0.20-33.zip"
このあと,[,2]の列にある情報をさらに使います。
ここの部分だけ,テキストエディタの矩形(くけい)選択で切り取ると便利です。
ダウンロードしたパッケージは,USBメモリなどで,インストールしたい
オフラインPCに移します。
(2)ネットにつながってないパソコンでパッケージをインストールする
ここからは,パッケージをインストールしたい,オフラインPC
での作業です。USBメモリから,こちらでも同じくC:/pkgsのフォルダに入れます。
ファイルを置く場所はどこでもよいのですが,分かりやすさのために上記のように
しています。日本語がなるべく入ってない方が失敗しにくいと思います。
そして以下のコードを実行します。(1)で出てきた[,2]列の情報をそのまま
コピペできると簡単です。
"C:/pkgs/MASS_7.3-45.zip",
"C:/pkgs/msm_1.6.1.zip",
"C:/pkgs/polycor_0.7-8.zip",
"C:/pkgs/survival_2.39-5.zip",
"C:/pkgs/mvtnorm_1.0-5.zip",
"C:/pkgs/expm_0.999-0.zip",
"C:/pkgs/sfsmisc_1.1-0.zip",
"C:/pkgs/Matrix_1.2-6.zip",
"C:/pkgs/lattice_0.20-33.zip"),
repos = NULL, conrib.url="win.binary", type = "win.binary")
これで完了です。
StataやSASにある変数のkeepやdropをRで行うには
※2018/02/10追記 今なら全部tidyverseで解決です。Rメモ(tidyverse)目次 参照
データハンドリングの過程で,コンピューターの負担を軽くしたり
データセットの中身を把握しやすくしたりするために,不要な
変数を削る,または必要な変数だけ残すなどする時がよくあります。
StataやSASでは,それぞれdropやkeepという命令で簡単にできた
のですが,Rでは少し面倒で以下のようにするみたいです。
また,dropはEZR,keepはRコマンダーにそれぞれ対応する箇所があった
ので,その情報も書いておきます。
irisデータでやってみます。特に問題はないのですが,testデータセットを新たに
作っておきます。
head(test)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
■ drop
test$Sepal.Width <- NULL
head(test)
Petal.Length Petal.Width Species
1 1.4 0.2 setosa
2 1.4 0.2 setosa
3 1.3 0.2 setosa
4 1.5 0.2 setosa
5 1.4 0.2 setosa
6 1.7 0.4 setosa
EZR(Version 1.20)では
Rコマンダー(バージョン2.0-3)では
■ keep
head(test)
Petal.Width Species
1 0.2 setosa
2 0.2 setosa
3 0.2 setosa
4 0.2 setosa
5 0.2 setosa
6 0.4 setosa
EZR(Version 1.20)では
指定した条件を満たす行だけを抽出したデータセットを作成する
Rコマンダー(バージョン2.0-3)では
■条件に合ったケースをdrop
変数のdropではなく,例えば女性だけのデータセットにするなど、
変数の条件に合ったケースを削除したい場合。
まずirisデータのケース数を確認すると、
[1] 150
であり,変数Speciesの内訳をみると
setosa versicolor virginica
50 50 50
となっていることを確認します。このうち、setosaのケースだけを
dropしたいとします。
これで完了です。またケース数を確認してみます。
[1] 100
このように,setosaの50ケースが全て削除されました。
setosa versicolor virginica
0 50 50
setosaという水準名は残っていますが(因子型変数なので),
数は0となっています。
詳しくはhttp://www.statmethods.net/management/subset.html参照(Quick-R)
変更履歴
2014/03/11 作成
2014/03/13 追加
2015/09/02 ■条件に合ったケースをdropを追加
測定誤差:測定の標準誤差(SEM:standard error of measurement)
再検査信頼性でICCを算出したら,同時に
測定誤差についても算出し報告することが必要です。
測定誤差の情報があれば,尺度の得点の解釈が非常に
分かりやすくなるからです。
COSMINチェックリストマニュアルによれば,測定誤差とは,
「測定された構成概念の真の変化に起因しない,得点の系統的および
ランダムな誤差のこと」をいいます。この尺度特性において求められる
統計指標の一つが,測定の標準誤差(SEM:standard error of measurement)
になります。
これについて包括的に解説している有名な論文は↓になります。
When to use agreement versus reliability measures.
de Vet et al. 2006 J Clin Epidemiol. 2006 59(10):1033-9.
この論文によれば,SEMの算出法には次のように3通りあります。
(1)ICCの計算過程の値を利用する:誤差分散の二乗根
(2)σ√(1-ICC)
(3)SDdiff/√2
しかし,SEM自体にも,ICCのようにagreementとconsistency
の区別があり,再検査信頼性の場合には,系統的誤差を考慮した
agreementが望ましいらしいです。
上記3つの算出法の内,(2)(3)は,consistencyのみ
しか出せないので,必然的に(1)の方法でSEMを算出するしかありません。
論文によっては(2)(3)をこの文脈で使ってる場合もありますが、
おそらく望ましいのは(1)と思われます。
では(1)の算出法についてRを用いた計算を紹介します。
参考にしたのは,奥村先生の公開スライド
「臨床的有意性の書き方」
臨床疫学研究における報告の質向上のための統計学の研究会
第16回研究集会 2014/8/2 (土)
の22枚目からの数ページです。
データセットは,ICCの記事
再検査(再テスト)信頼性のための級内相関係数(ICC)(3):Rで算出
で紹介した同じデータを用意します。ICCの続きで本記事を見ていてすでに
変数を作成している方は以下2行はとばしてかまいません。
set.seed(123);st2<-round(st1 + runif(20, min=-5,max=10))
SEMの算出のためには,データセットをwide→longに変換する
作業が必要です。また,その際に変数名にはルールが求められ,変数名
の最後が「.数字」で終わる必要があります。したがって,データフレーム
にまとめる過程で以下のように定義します。※ 今はpivot_longerを使った方が簡単(2021/10/20追記)
次に,ID番号をつけます。
そして,データセットをwideからlongに変換します。
varying = c("st.1", "st.2"),
direction = "long")
念のため確認,st.1とst.2が縦にくっついて1と2の数字部分が
新しくできたtimeという変数に移動している様子が分かります。
id time st
1.1 1 1 17
1.2 1 2 16
2.1 2 1 19
2.2 2 2 26
3.1 3 1 32
3.2 3 2 33
そして,分散を分解して表示するために,線形混合モデル
を実行します。
fitlong <- lmer(st~(1|id)+(1|time), data=stdat_long)
そして分散を確認するには次のコマンドです。
Groups Name Std.Dev.
id (Intercept) 7.1635
time (Intercept) 2.1855
Residual 3.1778
ここで出てくる数値はSDなので,これらを二乗する手計算をして
オブジェクトにそれぞれ格納します
sigma2time
4.77641
sigma2resid
10.09841
そしてこれらを足した数の二乗根を求めればSEMが算出
できます。
3.843963
もう1つの方法として、ANOVAの結果のmean squaresより計算する
こともできます。
パッケージpsychのICC関数では,以下のようにして
分散分析表の結果も出力できます。
res$summary
Df Sum Sq Mean Sq F value Pr(>F)
subs 19 2141.9 112.7 11.16 1.21e-06 ***
ind 1 105.6 105.6 10.46 0.00437 **
Residuals 19 191.9 10.1
この内,indとResidualsのMean Sqを用います。
(MS(ind)-MS(Residuals))/nの式に当てはめると,
4.775
この値とMS(Residuals)を足した値の二乗根がSEMになるので
3.856812
ということで,線型混合モデルで算出した値とほぼ一致します。
SEMを出したら,SDC:smallest detectable change
についても算出すると,得点を解釈しやすくなります。
SDCは,単純にSEM*√2*1.96をすれば出てきます。介入前後で測定して
この値を超えて得点が変化していれば,統計的に意味のある変化であると
解釈できるということになります。
上記の例を使うと,
10.69054
なので,この尺度の得点が11点を超えて変わっていれば,測定誤差
を超えて変化していると解釈できるということです。
以下,統計解析部分の記載例です
中耳炎-6質問票(OM-6):因子構造と解釈可能性を重視した心理測定学的特性
Heidemann et al. (2013) Health Qual Life Outcomes. 11:201.
It was estimated by computing the square root of the within subject variance
of the patients (SEMagreement = √σbetween measurement + σresidual) [文献].
Variance components were obtained from a multilevel mixed effects model
(restricted maximum likelihood estimates) [文献].
変更履歴
20150817 head()の中が間違っていた部分を修正
20150831 ANOVAの結果を使ったSEMの算出法とSDCの算出法追記
20211020 pivot_longerとICC関数の部分追記
再検査(再テスト)信頼性のための級内相関係数(ICC)(3):Rで算出
仮想データとして,11項目4件法(1~4点)のストレス症状尺度を,同じ対象者が
間をあけて2回回答する状況をイメージしてみます。
最初の測定がst1,2回目の測定がst2としたサンプルデータを作成します。
st1では平均=21,SD=7の分布になるようにデータを20件作成します。
st2では,st1に範囲-5~10の一様分布の乱数を加えて少しばらつくように
します。それぞれ,set.seedを設定して,繰り返しても同じ乱数が得られる
ようにしています。多少得点高めに系統的にずれる感じを演出。
set.seed(123);st2<-round(st1 + runif(20, min=-5,max=10))
st1とst2をデータフレームにまとめます。
head(stdat)
st1 st2
1 17 16
2 19 26
3 32 33
4 21 29
5 22 31
6 33 29
まず相関係数をみてみます。
[1] 0.855905
図はこんな感じ
それでは級内相関係数を出してみます。
代表的なパッケージが2つあるのでそれぞれで。
■パッケージirrのicc関数
icc(stdat, "twoway", "agreement")
#icc(cbind(st1,st2), "twoway", "agreement")
Single Score Intraclass Correlation
Model: twoway
Type : agreement
Subjects = 20
Raters = 2
ICC(A,1) = 0.775
F-Test, H0: r0 = 0 ; H1: r0 > 0
F(19,7.78) = 11.2 , p = 0.000947
95%-Confidence Interval for ICC Population Values:
0.368 < ICC < 0.916
ICC(A,1)=ICC(2,1)ですので,ICC 0.775 95%CI(0.368,0.916)という
結果です。
■パッケージpsychのICC関数
library(psych)
ICC(stdat)
Call: ICC(x = stdat)
Intraclass correlation coefficients
type ICC F df1 df2 p lower bound upper bound
Single_raters_absolute ICC1 0.77 7.6 19 20 1.7e-05 0.51 0.90
Single_random_raters ICC2 0.78 11.2 19 19 1.2e-06 0.37 0.92
Single_fixed_raters ICC3 0.84 11.2 19 19 1.2e-06 0.63 0.93
Average_raters_absolute ICC1k 0.87 7.6 19 20 1.7e-05 0.67 0.95
Average_random_raters ICC2k 0.87 11.2 19 19 1.2e-06 0.54 0.96
Average_fixed_raters ICC3k 0.91 11.2 19 19 1.2e-06 0.77 0.96
Number of subjects = 20 Number of Judges = 2
すべての結果が出てきますが,ICC(2,1)は
Single_random_ratersの欄になります。したがって,
ICC 0.78 95%CI(0.37,0.92)という結果です。
ピアソンの積率相関の結果とは結構違うことがわかります。
関連記事として
再検査(再テスト)信頼性のための級内相関係数(ICC)
再検査(再テスト)信頼性のための級内相関係数(ICC)(2):Stataで算出
も参照ください。