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　追記
•  
   

【図・グラフ】

• irisデータでグラフ作成(1)―散布図（plot）
• irisデータでグラフ作成(2)―層別散布図（plot）修正中
  
• irisデータでグラフ作成(3)―回帰直線（全体，および品種を共変量として）
• irisデータでグラフ作成(4)―層別回帰直線（交互作用の確認グラフ）
• 散布図を描くときのデータの重なりをずらして見やすくする
• フォレストプロットをRのggplot2で描く  

（20201127　出力を画像に変更）

 

データを把握したり可視化したりする際に，一度にたくさんの項目の度数分布を出したい時があります。そんなときに，指定した変数のすべてについて実行する方法があると便利です。もっといいやり方がありそうですが，とりあえず自分で見つけた方法について紹介します。

データとして，psychパッケージのbfiデータを使います。度数分布を出すには，table( )でもよいですが，tabyl( )の方が便利なので，janitorパッケージを使います。

janitorパッケージでtidyな度数分布表やクロス表を作成

まず最初に，格納用の空のリストを作っておきます。resと名前をつけてます。
次に，変えていきたい変数一覧を準備します。ここではavarsと名前をつけました。

"A1","A2"のところは，変数名の位置で指定した

としても同じです。

 

 

これらの結果はリストresの中で別々の要素なので，最後に1つのデータフレームにまとめた方が便利です。そこで，どの質問項目の結果か分かるようにするために，項目名をそれぞれ値とした新たな変数を作ります。これは色々試した結果，以下のようにしないとできなかったので，とりあえずですが，次のようにします

 

 

 

最後に，それぞれのリストの要素を1つのデータフレームにまとめます。

 

 

 

関連情報：
ggplot2で変数だけをloopで変えていきながら同じ設定の図をたくさん描く

有名なタイタニック号の生存についてのデータを使います。
加工せずにすぐに使えるデータとして，effectsパッケージのTitanicSurvivalデータセットを呼び出します。全体のデータの一部分のようなので，少なめでn=1,309となっています。

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()で結果を見ることができます。

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()の中に入れます。

exp(coef(titanic)) #オッズ比

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.display(titanic)

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を足すと，すっきりした結果が出てきます。

                          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

（20170603 追記）
broomパッケージを使う方法もありました。tidyさを追求するならこちらがよいかも。

               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

以下の表にあるようなオッズ比と信頼区間のデータをフォレストプロットの
グラフにしたい場合にどうやるか説明します。

なお，ここでのグラフ作成は，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

Theorell et al., 2015 BMC Public Health 15:738. DOI: 10.1186/s12889-015-1954-4 CC BY 4.0
のFig.2よりsummary odds ratioを抜粋

まずデータをRに取り込むために，表をすべて選択して，excelなどに貼り付けます。
それをコピーした後に


とするとデータが取り込みできます。または，以下のコードを実行しても同様です。


次にパッケージをロードしますが，以下のコードはパッケージが入っているか
どうか確認し，なければインストールもしてくれるので便利です。

当然ggplot2は入ってるよ，という方は単に読みこむだけです。


以下のコードが，グラフ描画の本体です。




このとおり，綺麗な図が作成されます。
全体の見栄えを変えたい場合は以下のようにするとできます。



今後，見栄えの変え方などいろいろ試していきたいと思います。

Rで項目反応理論をやってみるだけならすごく簡単です。
以下のサイトで説明されている通りにコピペすればすぐに
できます。


ここでは，もっとかゆい所に手が届くような部分を
説明していきます。基礎知識は自分がまとめた下記スライド参照。


ここでは，パッケージ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つの図に出力したい場合は次のようにします。




その他の主要な指標は以下の通りです。
テスト情報関数(total information function[TIF]) ※




項目情報関数(item information function[IIF])　（2016/10/15　引数を修正）




というわけで，たいていのことは簡単にできてしまいます。
自分のデータセットで試したい，と思った場合に，恐らくいちばん
引っかかるのが，変数の型がfactorになってないなどの点だと
思います。それは，以下の説明を参考にするとうまくいくかも
しれません。



※standard errorと一緒に描画したい場合

【履歴】
2016/10/15　項目情報関数の引数がテスト情報関数と同じだったのを修正

例えば年齢から20代，30代といった年齢層に分けた変数を作るような時に便利
なのがcut関数です。

(1) 単純な例（0～10の値）

 

 [1]  0  1  2  3  4  5  6  7  8  9 10

ここで，5以下と6以上で分けたいときは以下のようにします。

 

[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にしてみると，

 

 [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にしたい場合は次のようにします。

 

 [1] 1 1 1 1 1 1 2 2 2 2 2
Levels: 1 2

※20160918修正。ラベルを0,1としていたが，ラベルが何であっても内部的には1,2
が当てられているので，それを反映するように修正。確認するには

ちなみに，引数のrightがデフォルトでTRUEになっているので，これを
FALSEにすると，次のようにbreaksを設定する必要があります。

 

 [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も。

 

 [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つの水準に区分してみます。

 

    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でおきかえてもよいです。
左端にはマイナス（-）が付きます。

ラベルを変更する場合は以下のように

    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パッケージを指定してみます。

[1] "ltm"      "MASS"     "msm"      "polycor"  "survival"
[6] "mvtnorm"  "expm"     "sfsmisc"  "Matrix"   "lattice"

これで，ltmを使う際には，他の9個のパッケージも同時にインストールしておいた
方がよさそうだと分かります。

さて，ここからが上記リンク先と異なり，WindowsのPCを想定した
手順となります（自分が理解するまで苦労した所らへんです）。

まず，pkgsオブジェクトに格納されているパッケージをすべてダウンロードします。

つまり


です。実行すると，指定したフォルダにパッケージがダウンロードされ，
以下のような出力が出ます。

      [,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]列の情報をそのまま
コピペできると簡単です。


これで完了です。

※2018/02/10追記　今なら全部tidyverseで解決です。Rメモ（tidyverse）目次　参照

データハンドリングの過程で，コンピューターの負担を軽くしたり
データセットの中身を把握しやすくしたりするために，不要な
変数を削る，または必要な変数だけ残すなどする時がよくあります。

StataやSASでは，それぞれdropやkeepという命令で簡単にできた
のですが，Rでは少し面倒で以下のようにするみたいです。
また，dropはEZR，keepはRコマンダーにそれぞれ対応する箇所があった
ので，その情報も書いておきます。

irisデータでやってみます。特に問題はないのですが，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

  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

  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を追加
 

再検査信頼性で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行はとばしてかまいません。

SEMの算出のためには，データセットをwide→longに変換する
作業が必要です。また，その際に変数名にはルールが求められ，変数名
の最後が「.数字」で終わる必要があります。したがって，データフレーム
にまとめる過程で以下のように定義します。※ 今はpivot_longerを使った方が簡単（2021/10/20追記）

次に，ID番号をつけます。

そして，データセットをwideから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

そして，分散を分解して表示するために，線形混合モデル
を実行します。

そして分散を確認するには次のコマンドです。

 

 Groups   Name        Std.Dev.
 id       (Intercept) 7.1635 
 time     (Intercept) 2.1855 
 Residual             3.1778 

ここで出てくる数値はSDなので，これらを二乗する手計算をして
オブジェクトにそれぞれ格納します

 

4.77641

 

 

10.09841

そしてこれらを足した数の二乗根を求めればSEMが算出
できます。

 

3.843963

もう１つの方法として、ANOVAの結果のmean squaresより計算する
こともできます。

パッケージpsychのICC関数では，以下のようにして
分散分析表の結果も出力できます。

 

            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関数の部分追記

仮想データとして，11項目4件法（1～4点）のストレス症状尺度を，同じ対象者が
間をあけて2回回答する状況をイメージしてみます。
最初の測定がst1，2回目の測定がst2としたサンプルデータを作成します。

st1では平均=21，SD=7の分布になるようにデータを20件作成します。
st2では，st1に範囲-5～10の一様分布の乱数を加えて少しばらつくように
します。それぞれ，set.seedを設定して，繰り返しても同じ乱数が得られる
ようにしています。多少得点高めに系統的にずれる感じを演出。

st1とst2をデータフレームにまとめます。

 

  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関数

 

 

 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関数

 

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で算出
も参照ください。