004703
統計学関連なんでもあり

統計学に関する質問とか意見などなんでも書き込んでください。回答があるか,回答が正しいかは保証できません。
レポートや宿題の答えを求める記事や,その他不適切と判断した記事は,管理者の独断で即刻削除します。
ローマ数字,丸付き数字などのフォントセット依存文字および半角カタカナは使わないでください。
記事を引用する際には,適切に引用してください(全文引用はしないでください)。
問題が解決した(結局解決しなかった)場合は,その旨報告していただければ,コメントした人や他の読者に参考になるでしょう。


[トップに戻る] [利用上の注意] [ワード検索] [過去ログ] [統計学自習ノート] [管理用]
latest article: No. 22654, 2018/12/05(Wed) 18:18
おなまえ
タイトル
コメント
URL
添付ファイル
暗証キー (英数字で8文字以内)

自分の投稿記事を修正・削除する場合 --> 処理 記事No 暗証キー


名義尺度の最小二乗フィッティングで追加制約がある場合の標準誤差の計算方法
  投稿者:小西 2018/12/05(Wed) 18:18 No. 22654
長文失礼致します。
以下についてお分かりになる方がいましたらご教授頂けますでしょうか?

現在、交互作用のない2因子名義尺度に対する最小二乗フィッティングを行うプログラムを作成しているのですが、特殊な制約を加えた場合の標準誤差の計算方法がわからず困っております。
目標は2種類の加工装置の各個体が製品のばらつきにどの程度寄与しているかを推定するため、平均値(切片)、ばらつきがどの個体に起因しているか(回帰係数)、それら係数の信頼区間を求めることです。
例えば装置Aが3水準(A1,A2,A3)、装置Bが4水準(B1,B2,B3,B4)あるとして、AとBの組み合わせごとに製品の測定値が与えられます。ただし、いくつかの組み合わせが欠損値となることもあり得ます。

基本的にはJMPの計算結果を参考にしているのですが、測定値の偶然性や欠損値などによってデータが線形従属になった(解に自由度がある)場合、JMPでは一部の水準を除外(係数を0に固定)する計算方法を採用しています。
この方法では入力データの水準の並び順が変わっただけで除外される水準が変化してしまいますし、特定の水準を特別扱いしていることにもなるので、作成中のプログラムではこの点を改良したいと考えています。
※JMPにおけるランク落ちの対処法の詳細は
 http://www.jmp.com/japan/support/help/13/flm-sls.shtml
 の下部にある「1次従属性がある場合」のサブトピックで説明されています。

そこで、特定の水準を特別扱いせずに解を一意に決定する方法として、各個体の回帰係数(切片は除く)の2乗和を最小にする制約を設けることを考えました。
(理論的根拠をもとに考えたわけではないので、統計的に不適切だったりより良い方法があるかもしれません)
そして、等式制約付き最小二乗問題を解くことで、そのような回帰係数を求めることもできました。
しかし、信頼区間を求めるのに必要な標準誤差をこの制約下でどのように計算すればよいかがわかりません。※入力データによっては自由度が不足して標準誤差を求められないケースがあることは承知しております。
この場合の標準誤差はどのように計算すればよいでしょうか?(下記詳細もご参照ください)

■参考:JMPと同等な計算方法
・入力データ例(製品測定値と使用装置) ※線形従属になる一例
 11 A1 B1
 13 A1 B2
 12 A2 B1
 15 A2 B2
 14 A3 B3
 10 A3 B4

 1列目の測定値ベクトルをy、データ数をnとします。
 ※この例では各組合せの測定値が高々1個になっていますが、実際には制限はありません。

・各因子水準をコード化(因子内の回帰係数の和は0に制約)
     A13 A23
 A1 →  1  0
 A2 →  0  1
 A3 →  -1  -1
 これを行列Aとします。

     B14 B24 B34
 B1 →  1  0  0
 B2 →  0  1  0
 B3 →  0  0  1
 B4 →  -1  -1  -1
 これを行列Bとします。

・コード化した入力データ
 切片 A13 A23 B14 B24 B34
  1  1  0  1  0  0
  1  1  0  0  1  0
  1  0  1  1  0  0
  1  0  1  0  1  0
  1  -1  -1  0  0  1
  1  -1  -1  0  -1  -1

・線形従属列の除外
 コード化した入力データの1列目、1〜2列目、1〜3列目、1〜4列目、1〜5列目、1〜6列目と順にランクを調べると6列目が増えたときにランク落ちしていることがわかるので、6列目を削除します。
 ※より大きな入力データの場合は後続の列も同様に繰り返します。

 切片 A13 A23 B14 B24
  1  1  0  1  0
  1  1  0  0  1
  1  0  1  1  0
  1  0  1  0  1
  1  -1  -1  0  0
  1  -1  -1  0  -1
 これを入力データ行列X、Xのランクをrとします。

・削減した係数を復元する行列
 ブロック行列
 (1 0 0)
 (0 A 0)
 (0 0 B)
 から上記で除外した線形従属列(この例ではB34列)を同様に除外した行列をZとします。

・最小二乗問題Xa=yを解く
 解ベクトルaから切片,A1,A2,B1,B2が決まり、Zaですべての係数が求まります。
 平均二乗誤差:MSE=(y'y-a'X'Xa)/(n-r)
 標準誤差:SE=sqrt(MSE*diag(Z*inv(X'X)*Z'))
 ※実際にはQR分解を利用してこれらと等価な計算を行います。

■詳細:検討中の2乗和を最小化する計算方法
※以下ではI_kはk次単位行列です。

・各因子水準のコード化(ここでは因子内の回帰係数の和を0に制約しません)
 A=I_3, B=I_4

・コード化した入力データ
 切片 A1 A2 A3 B1 B2 B3 B4
  1  1  0  0  1  0  0  0
  1  1  0  0  0  1  0  0
  1  0  1  0  1  0  0  0
  1  0  1  0  0  1  0  0
  1  0  0  1  0  0  1  0
  1  0  0  1  0  0  0  1
 これを入力データ行列X、Xのランクをrとします。

・等式制約付き最小二乗問題を解く
 Ca=0       ※切片以外の二乗和を最小化
 s.t. Da=0    ※因子内の回帰係数の和を0に制約
    X'Xa=X'y  ※正規方程式で最小二乗誤差に制約

 ここでCは以下のブロック行列
 (0 I_7)

 Dは以下の行列
 (0 1 1 1 0 0 0 0)
 (0 0 0 0 1 1 1 1)

 また、以下のブロック行列をEとします。
 ( D )
 (X'X)

 これにより解ベクトルaは求められますが、標準誤差SEの計算方法がわかりません。

 素人考えでJMPの計算式の形だけまねて、
 平均二乗誤差:MSE=(y'y-a'X'Xa)/(n-r)
 標準誤差:SE=sqrt(MSE*diag(pinv(X'X))) ※pinvは疑似逆行列
 としてみましたが、JMPで列の除外がない入力データの場合にJMPの結果と一致しませんでした。
 JMPの計算ではdiagの中の式(共分散行列?)にDa=0相当の制約条件が反映されているため、同様にDa=0相当およびEの零空間で二乗和を最小にする制約条件をdiagの中に反映する必要がある気がしています。


その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/01(Sat) 10:59 No. 22639
いろいろ枝葉末節にこだわりすぎた感がありました。
サンプルサイズが小さい場合に精度が低いのはやむを得ない(当然)なので,結局は,対数正規分布のパラメータとそれを対数変換した正規分布のパラメータの関係式から解を求めるしかないようですね。
対数正規分布の標本平均,標本分散から,非線形連立方程式で正規分布のパラメータを求め,それに基づき信頼限界を計算し,それを指数変換する。

元データなしで,信頼限界を求める
getCL = function(MEAN, VAR, n) {
func = function(p) {
return(c(exp(p[1]+p[2]/2)-MEAN, exp(2*p[1]+p[2])*(exp(p[2])-1)-VAR))
}
library(nleqslv)
ans = nleqslv(c(2,2), func) # 対数正規分布の平均値と分散から
mean = ans$x[1] # 対数をとったときの正規分布の平均値と
sd = sqrt(ans$x[2]) # 標準偏差を非線形連立方程式を解いて求める
cat(sprintf("mean = %f, sd = %g\n", mean, sd))
cl.normal = mean+c(1, -1)*qt(0.025, n-1)*sd/sqrt(n-1) # 正規分布目盛りでの信頼限界
cl = exp(cl.normal) # 対数目盛りに変換
cat(cl)
}
対数正規分布にしたがうデータから,標本平均・標本分散 MEAN = 3.390011 VAR = 37.060
が求められている(元データはない)

これから対数をとった正規分布の標本平均,標本標準偏差を求める
mean = 0.500339, sd = 1.20041
となる
それに基づいて,信頼区間を求めると
0.667087 4.077615

サンプルサイズが小さいと rlnorm で生成されたデータの平均値と分散は,母数とはかなりかけ離れたものになっている

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:Saito 2018/12/03(Mon) 10:29 No. 22642
コメントありがとうございます。
しかし、一つ前のスレッドでも書きましたが、こちらの手元にあるのは平均値とその標準「誤差」であり、元データの平均値と標準「偏差」ではありませんでした。そのため、微妙に私のやりたいこととずれていることに気づきましたので、再度ご教授頂ければ幸いです。何度も何度も申し訳ありません。

まず、元データの平均値と標準偏差が使える場合に合うことを示します。

> set.seed(1)
> n <- 100000
> u <- 1
> sig <- 0.5
>
> #生物の密度データがあるとする(非負)。平均密度とその信頼区間を知りたい。
> y <- exp(rnorm(n, u, sig))
> hist(y)
>
> #yというデータがあるもとで、正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
>
> res2 <- summary(lm(log(y)~1))
> lo2 <- exp(res2$coef[1]-1.96*res2$coef[2])
> up2 <- exp(res2$coef[1]+1.96*res2$coef[2])
>
> MEAN = mean(y)
> VAR = var(y)
> lo2;up2
[1] 2.706802
[1] 2.723691
> getCL(MEAN, VAR, n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906>

上記のように、MEAN=mean(y), VAR=var(y)として、元データyの平均値と分散が使えるという状況であれば、getCLとlo2, up2がほぼ一致することを確認しました(なお、sigが大きいと合わないようなので小さくし、小サンプルであることの問題も省きたいため、サンプルサイズnは大きくしてあります)

ところが、私の手元にあるのは、yの平均値と標準偏差ではなく、解析結果であるyの平均値とその標準誤差となっています。つまり、res$coef[1]とres$coef[2]が手元にあるということでして、mean(y)とvar(y)の値が手元にあるわけではありませんでした・・・。私の説明が拙かったようで、本当に申し訳ありません。
確認してみましたが、当然ながら、こちらとは合いません。

> MEAN = res$coef[1]
> VAR = res$coef[2]^2
> lo2;up2
[1] 2.706802
[1] 2.723691
> getCL(MEAN, VAR, n)
mean = 1.124605, sd = 0.0016874
3.07897 3.079034>

何度もお尋ねするのは心苦しいのですが、res$coef[1]とres$coef[2]が手元にあり、mean(y)とvar(y)の値(とn)が手元にない条件で、lo2, up2に準じるものは出せるでしょうか。何度も申し訳ありません。

どうぞよろしくお願いいたします。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/03(Mon) 11:23 No. 22643
> 私の手元にあるのは、yの平均値と標準偏差ではなく、解析結果であるyの平均値とその標準誤差となっています

標準誤差=標準偏差/sqrt(サンプルサイズ) ですけど?
標準偏差=標準誤差*sqrt(サンプルサイズ)
分散=標準偏差^2

でしょう?

> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
lm する必要なんてなくて,
res$coef[1] = mean(y)
res$coef[2] = sd(y)/sqrt(n)
ですよね。

簡単なところを複雑にやっているので,頭が混乱します

> set.seed(1)
> n <- 100000
> u <- 1
> sig <- 0.5
> #生物の密度データがあるとする(非負)。平均密度とその信頼区間を知りたい。
> y <- exp(rnorm(n, u, sig))
> hist(y)
> #yというデータがあるもとで、正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
> res <- summary(lm(y~1))
> lo1 <- res$coef[1]-1.96*res$coef[2]
> up1 <- res$coef[1]+1.96*res$coef[2]
> lo1
[1] 3.068823
> up1
[1] 3.089189
> res$coef[1]
[1] 3.079006
> mean(y)
[1] 3.079006
> res$coef[2]
[1] 0.005195524
> sd(y)/sqrt(n)
[1] 0.005195524
> (res$coef[2]*sqrt(n))^2
[1] 2.699347
> var(y)
[1] 2.699347
>

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/03(Mon) 11:39 No. 22644
つづき

> mean(y)
[1] 3.079006
> var(y)
[1] 2.699347
> n
[1] 1e+05

以下で mean, var を使っているけど yは必要ない。平均値と分散を引数として与えているだけ。

> getCL(mean(y), var(y), n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906

以下のように必要なのは数値だけ
ただ,ここで示したように有効数字7桁では正確ではないので,コンピュータ上で持っている正確な値を使っただけ。
あなたは,y は持っていないが,計算された平均値と分散(をサンプルサイズで割ったものの平方根である標準偏差)を持っているので,その二つの数値を指定するだけ。

> getCL(3.079006, 2.699347, 100000)
mean = 0.999331, sd = 0.500551
2.70805 2.724906

標準誤差と分散の関係をちゃんとすれば,res$coef[1], res$coef[2] を使っても,当然同じ結果になる

> res$coef[1]
[1] 3.079006
> (res$coef[2]*sqrt(n))^2
[1] 2.699347


> getCL(res$coef[1], (res$coef[2]*sqrt(n))^2, n)
mean = 0.999331, sd = 0.500551
2.708051 2.724906 当然ながら同じになる

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/03(Mon) 13:49 No. 22647
また枝葉末節なこと(基本的なこと)になりますが,あなたが使っている 1.96 は標準正規分布由来のものだが,本来は t 分布に基づく値でなければならない。

http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean2.html
の「母分散が不明の場合」にあたるので。

いまはテストデータでやっているから n がとても大きく,1.96 で問題ないが,実際に運用するときには問題になる。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:Saito 2018/12/04(Tue) 08:56 No. 22650
コメントありがとうございます。
不勉強で申し訳ありません。理解できていないのですが、標準誤差から標準偏差に戻すときに、nが未知でも戻るのでしょうか?今手元にあるのは、平均値とその標準誤差のみで、サンプルサイズnは不明なのです。getCLの中身を見る限り、何だか無理なような気がしてきています・・・。
もし戻らないとなると、当初の目的を果たすことは不可能でしょうか。

t分布については失念していました。
n=10だと、2.228ですね。ありがとうございます。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:青木繁伸 2018/12/04(Tue) 09:38 No. 22652
> 標準誤差から標準偏差に戻すときに、nが未知でも戻るのでしょうか?

定義式に書かれているとおり。
nが分からないなら無理でしょう。
現代統計学では,サンプルサイズこそ一番大事なものです。

妥当と思われるnを仮定するか,何通りかのnで代替してこの範囲という結果を示すしかないでしょう。

Re: その2 正規分布を仮定した信頼区間を。。。
  投稿者:Saito 2018/12/04(Tue) 11:25 No. 22653
コメントありがとうございます。
nがわからなければ不可能という旨、承知いたしました。
色々とご教授くださりありがとうございました。

以下余談です。
何故こうなるのかわからず、おそらく間違っていると思って相談しませんでしたが、この件について私ではない方(前任者)のメモで

Lower=exp(log(res$coef[1])-1.96*res$coef[2]/res$coef[1])
Upper=exp(log(res$coef[1])+1.96*res$coef[2]/res$coef[1])

として、平均値を対数を取って、なぜかCVを足す(引く)というやり方で正規分布を仮定して得られた平均値と標準誤差(res$coef[1], res$coef[2])を使って、対数正規分布を仮定したときの95%信頼区間(Lower, Upper)になるというメモがありました。やってみても、getCLと全然違う範囲になりますし、CVという単位が消えた値を足すのも奇妙なので、間違っているのでは、と思っています。とはいえ、他の用途でも構いませんので、この式に見覚えがありましたら、ご教授くだされば幸いです。


片側検定の帰無仮説
  投稿者: 2018/12/03(Mon) 09:19 No. 22640
例えば、母比率の差の検定で、
対立仮説が「A≠B」の場合、
帰無仮説は、両側検定の場合、「A=B」と認識しています。
対立仮説が「A>B」の場合、
片側検定の場合の帰無仮説は、
「A=B」でよいのでしょうか?
厳密には「A≦B」と書くべきではないでしょうか?

Re: 片側検定の帰無仮説
  投稿者:青木繁伸 2018/12/03(Mon) 10:25 No. 22641
検定統計量は,帰無仮説が正しいつまり,母数が特定の一つの値であると仮定したときのもの
片側検定でも,帰無仮説は「母数=特定の値」
http://aoki2.si.gunma-u.ac.jp/lecture/Kentei/p-value.html
の図で,「H0の下での統計量の分布」というのは,たとえば A=B のときに計算できるもの。A≦B というようなとき,分布は特定できないでしょう?

片側検定でも,両側検定でも「H0の下での統計量の分布」は同じものです。
どの部分の面積(確率)をP値とするかが違うだけです。

Re: 片側検定の帰無仮説
  投稿者: 2018/12/03(Mon) 12:06 No. 22645
ご回答ありがとうございます。

同様の質問になりますが、
母比率の差の検定のp値の意味についてご教授頂けましたら幸いです。

両側検定も片側検定も、
実際には差がないのに、差があると誤ってしまう確率
ということになるのかと思うのですが、
片側検定の場合、例えば、
実際にはA≧Bなのに、A<Bと誤ってしまう確率
とすることは不適当なのでしょうか?

Re: 片側検定の帰無仮説
  投稿者:青木繁伸 2018/12/03(Mon) 13:15 No. 22646
検定の基本的なところなので,正確に述べる必要があります。
「似たようなこと」をいっているのでよいだろうと言うことにはなりません。
前提と叙述の順序と用語法がとてもだいじ。
======
P 値とは,「帰無仮説のもとで,観察された統計量より極端な値を取る確率」

http://aoki2.si.gunma-u.ac.jp/lecture/Kentei/p-value.html において

両側検定の場合は PL, PU の和

片側検定の場合は,対立仮説の方向によりPL か PU のいずれか

もし,P 値が(たとえば)0.05 より小さい場合は帰無仮説を棄却する。

両側検定の場合は,帰無仮説はたとえば 母比率=0.5 で,対立仮説は 母比率≠0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率≠0.5 と結論する。
しかし,帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P

片側検定の場合は,帰無仮説は同じく 母比率=0.5 で,対立仮説はたとえば 母比率>0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率>0.5 と結論する。
しかし,帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P

最後の行で,「帰無仮説が正しい(母比率≦0.5)」ではないことに留意

Re: 片側検定の帰無仮説
  投稿者: 2018/12/03(Mon) 15:39 No. 22648
ご教授頂きまして誠にありがとうございます。
少しずつわかった気がいてきました。

両側検定の場合は,帰無仮説はたとえば 母比率=0.5 で,対立仮説は 母比率≠0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率≠0.5 と結論する。

帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
これを
帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
としてもよいのでしょうか?

同様に

片側検定の場合は,帰無仮説は同じく 母比率=0.5 で,対立仮説はたとえば 母比率>0.5。
帰無仮説が棄却されれば対立仮説が正しいと考え 母比率>0.5 と結論する。
帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
これを
帰無仮説が正しい(母比率=0.5)のに母比率>0.5 とする確率は P
としてもよいのでしょうか?
それとも
帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
でしょうか?

Re: 片側検定の帰無仮説
  投稿者:青木繁伸 2018/12/03(Mon) 17:16 No. 22649
> 帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
> これを
> 帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
> としてもよいのでしょうか?

「帰無仮説を棄却する」は
「対立仮説を採択する」かつ「対立仮説は母比率≠0.5」なので「母比率≠0.5とする」
と同じですから,よいでしょう。

> 帰無仮説が正しい(母比率=0.5)のに帰無仮説を棄却する確率は P
> これを
> 帰無仮説が正しい(母比率=0.5)のに母比率>0.5 とする確率は P
> としてもよいのでしょうか?

両側検定の場合と同じです。よいでしょう。

> それとも
> 帰無仮説が正しい(母比率=0.5)のに母比率≠0.5とする確率は P
> でしょうか?

「母比率>0.5 とする確率」と「母比率≠0.5とする確率」は違うので,後者は間違い。
そもそも,方向性のある片側対立仮説との比較なのに,なぜ両側仮説が出てくるのか?
むりやり考えると
「母比率>0.5 とする確率」と,逆側の「母比率<0.5 とする確率」を合計したものが 2P ということか。
つまり,片側対立仮説「母比率>0.5」も,「「母比率<0.5」も 統計学的には帰無仮説「母比率=0」から見れば同じく棄却域にあるので,下側確率と両側確率の和だが,今の場合分布は対称なので確率は 2P

Re: 片側検定の帰無仮説
  投稿者: 2018/12/04(Tue) 09:28 No. 22651
ありがとうございました。
だんだん理解できました。


正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法
  投稿者:Saito 2018/11/29(Thu) 17:11 No. 22635
お世話になっております。
標記の通りなのですが、生データがなく正規分布を仮定した際の平均値と標準偏差しかない状況で、対数正規分布を仮定していた場合の信頼区間を算出する方法について、私のコーディングが間違っているだけかもしれませんが、うまく推定できないので、どなたかご教授頂けると幸いです。

下記に例を示します。

> set.seed(1)
> n <- 10000
> u <- 0.5
> sig <- 0.8
> x <- exp(rnorm(n, u, sig))
>
> #1.平均値
> mean(x)
[1] 2.269916
>
> #2.標準偏差
> sd(x)
[1] 2.123908
>
> #3.信頼区間下限値
> mean(x)-1.96*sd(x)
[1] -1.892943
>
> #4.信頼区間上限値
> mean(x)+1.96*sd(x)
[1] 6.432775

このとき、手元にあるのは、1-2の平均値、標準偏差、そして95%信頼区間(の近似値)である、3-4の値のみで、生データであるxは利用不可能な状況であるとします。

上記の信頼区間の算出方法ですと、xが負の値を取りえないのに、マイナスの値をとってしまっていることが確認できると思います。この部分を正の値の下限値を出すようにしたいと考えています。

一応、対数正規分布の信頼区間から正規分布の信頼区間を得る方法が使えるのかと思い、下記のURLにある計算式を打ち込んでみますと、実際の95%点とずれが出るように思われます。
http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf

#推定値
> mean(x)/exp(1.96*sqrt(log(1+(sd(x)/mean(x))^2)))
[1] 0.4797148
> mean(x)*exp(1.96*sqrt(log(1+(sd(x)/mean(x))^2)))
[1] 10.7408

#実際の95%点
> quantile(x, c(0.025, 0.975))
2.5% 97.5%
0.3267348 8.0988751

たぶんこういう風に逆算でやるのはダメ、ということなのだと思うのですが・・・。
試しにlogの世界でやってみるとおおよそ合います。
> quantile(log(x), c(0.025, 0.975))
2.5% 97.5%
-1.523258 2.489656
> mean(log(x)-1.96*sd(log(x)))
[1] -1.490756
> mean(log(x)+1.96*sd(log(x)))
[1] 2.477682

どこかで間違っているのだと思いますが、どなたかmean(x)とsd(x)だけを使って、quantile(x, c(0.025, 0.975))で出力される95%点(負の値を取らない)を出すやり方をご存知の方がおられましたら、ご教授ください。

どうぞよろしくお願いいたします。

Re: 正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法
  投稿者:青木繁伸 2018/11/29(Thu) 19:46 No. 22636
いろいろあるので,順を追って

set.seed(1)
n <- 10000
u <- 0.5
sig <- 0.8
としたとき,
x0 <- rnorm(n, u, sig)
の x0 の平均値と標準偏差は当然ながら 0.5 と 0.8 ではないということ

> x0 <- rnorm(n, u, sig)
> mean(x0)
[1] 0.4947704
> sd(x0)
[1] 0.8098852

正確に指定された平均値と標準差を持つ正規乱数を生成するのは,一度正規化した後に,標準偏差を掛けて平均値を加える(標準化の逆をやる)。

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8

この指数を取って対数正規分布にする。

> x2 <- exp(x1)
> mean(x2)
[1] 2.268486
> sd(x2)
[1] 2.130959

しかし,平均値と標準偏差を計算するのはそもそもおかしい。対数正規分布において計算される平均値と標準偏差のもつ意味は,正規分布のそれらが持つ意味とは全く異なる。たとえば,対数正規分布の平均値±1.96*標準偏差の範囲内にあるデータの割合は95%などではない。

また,mean(x)-1.96*sd(x) を「信頼区間下限値」としているが,それは上で書いたように,あるいは,あなたが最後の方で quantile といっていることからも,データの下から 2.5% 目の推定値でしょう。あなたは「母平均の信頼区間」を得たいのかな?ちがいますよね。

>> 試しにlogの世界でやってみるとおおよそ合います。
>> > quantile(log(x), c(0.025, 0.975))
とありますが,

あなたの対数正規乱数 x は,正規乱数の指数をとったもので,私が上で作った正確な例では x2 です。
> quantile(log(x2), c(0.025, 0.975))
2.5% 97.5%
-1.083910 2.063209
はあたりまえですが,以下と同じですね。 log(x2) == x1, x2 ==exp(x1) ですから。
> quantile(x1, c(0.025, 0.975))
2.5% 97.5%
-1.083910 2.063209

解はあるか?

その前に
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/wrong-sd.html
を読んでみてください。

そこにある図で,正規分布での平均値±標準偏差が対数正規分布では exp(平均値)*/exp(標準偏差) に対応すると書いてあります。

しかし,逆はできないのです。

あなたは図の上側の対数正規分布を対象にしているのですが,対数正規分布で計算した標準偏差は,図の下側の正規分布には変換するすべがないのが分かるかと思います。

解はあると思いますが,私は知らない。

http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf
は,正規分布のパラメータμとσと対数正規分布の平均値と分散の関係式であるが,解析的には成り立つが,実際のシミュレーションデータでは正確には成り立たない。
set.seed を外して何回かやってみると分かるが,x1 の 平均値,標準偏差はいつも指定された値(0.5, 0.8)になるが,それを exp したとたん,x2 の平均値,標準偏差は毎回変わる(つまり,解析的な値とは異なる.)。そのため,解析的な関係式にしたがわない。

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
>
> x2 <- exp(x1)
> mean(x2)
[1] 2.274771 ★これと
> sd(x2)
[1] 2.183789 ●これと

> x1 <- scale(rnorm(n))*sig+u
> mean(x1)
[1] 0.5
> sd(x1)
[1] 0.8
>
> x2 <- exp(x1)
> mean(x2)
[1] 2.266348 ★これ
> sd(x2)
[1] 2.08723 ●これ

Re: 正規分布を仮定した信頼区間を対数正規分布を仮定した信頼区間に直す方法
  投稿者:Saito 2018/11/30(Fri) 11:32 No. 22637
コメントありがとうございます。
色々と間違っていたようで申し訳ありません。
質問のコメントについて考えているうちに、私の質問自体がよくなかったことに気づきましたので、最後のほうで新しく例を作りました。二度手間にしてしまって申し訳ありません。私の質問としては、データではなく信頼区間について聞きたかった(ことに気づきました)ということです。申し訳ありません。大変恐れ入りますが、ご教授くだされば幸いです。

先に頂いたコメントについてご返信させていただきます。

> x0 <- rnorm(n, u, sig)
> の x0 の平均値と標準偏差は当然ながら 0.5 と 0.8 ではないということ

についてですが、おっしゃる通りです。ただ、中心極限定理で、0.5と0.8の周りに正規分布してくれるような方法でもOKです。

> しかし,平均値と標準偏差を計算するのはそもそもおかしい。対数正規分布において計算される平均値と標準偏差のもつ意味は,正規分布のそれらが持つ意味とは全く異なる。

仰る通りです。ただ、そういう出力(正規分布を仮定した平均値と標準偏差)がすでにされてしまって(紙になってしまって)いるので、生データを抜きに、それをなんとかしようということでした。

>また,mean(x)-1.96*sd(x) を「信頼区間下限値」としているが,それは上で書いたように,あるいは,あなたが最後の方で quantile といっていることからも,データの下から 2.5% 目の推定値でしょう。あなたは「母平均の信頼区間」を得たいのかな?ちがいますよね。

すみません、mean(x)-1.96*sd(x)はご指摘の通り信頼区間ではありませんでした。また、この指摘を受けて、私も聞きたいことが聞けていないことに気づきましたので、下で改めてシミュレーションデータを作らせていただきました。二度手間になってしまい、申し訳ありません。

>その前に
>http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/wrong-sd.html
>を読んでみてください。

上に同じです。読んでいくうちに、私の聞きたいことが聞けない質問内容になっていることに気づきました。本当に申し訳ありません。

>http://cse.fra.affrc.go.jp/okamura/memo/lognormal_ci.pdf
>は,正規分布のパラメータμとσと対数正規分布の平均値と分散の関係式であるが,解析的には成り立つが,実際のシミュレーションデータでは正確には成り立たない。

すみません、新しい質問を下記でさせていただく前に、ここについてお聞かせください。解析的に成り立たなくとも、中心極限定理で近い値にはなるのかと思っていたのですが、うまく推定できません。

set.seed(1)
n <- 10000
u <- 0.5
sig <- 0.8

B <- 1000
x <- matrix(exp(rnorm(n*B, u, sig)), n, B)

xlo1 <- apply(x, 2, mean)/exp(1.96*sqrt(log(1+(apply(x, 2, sd)/apply(x, 2, mean))^2)))
xup1 <- apply(x, 2, mean)*exp(1.96*sqrt(log(1+(apply(x, 2, sd)/apply(x, 2, mean))^2)))

xlo2 <- apply(x, 2, quantile, 0.025)
xup2 <- apply(x, 2, quantile, 0.975)

hist(xlo1, xlim=c(0.3, 0.55))
hist(xlo2, add=T, col=2)

hist(xup1, xlim=c(7, 13.5))
hist(xup2, add=T, col=2)

xlo1とxlo2、xup1とxup2はほぼ一致するのかと思ったのですが、どこか誤解しているようです。

以下似たような話ですが、データとパラメータで知りたいことが違っていましたので、質問を修正致します。大変申し訳ありません。

set.seed(1)
n <- 10
u <- 1
sig <- 2

#生物の密度データがあるとする(非負)。平均密度とその信頼区間を知りたい。
y <- exp(rnorm(n, u, sig))
hist(y)

#yというデータがあるもとで、正規分布を誤差分布に仮定したときの平均密度の95%信頼区間
res <- summary(lm(y~1))
lo1 <- res$coef[1]-1.96*res$coef[2]
up1 <- res$coef[1]+1.96*res$coef[2]
> lo1;up1
[1] -1.684872
[1] 22.9317

この信頼区間だと、下限がマイナスの値になってしまっていて、都合が悪いです。
ですので、誤差に対数正規分布を仮定して解析をやり直したものが下記になります。

#yというデータがあるもとで、対数正規分布を誤差分布に仮定したときの平均密度(元のスケール)の95%信頼区間
res2 <- summary(lm(log(y)~1))
lo2 <- exp(res2$coef[1]-1.96*res2$coef[2])
up2 <- exp(res2$coef[1]+1.96*res2$coef[2])
> lo2;up2
[1] 1.345521
[1] 9.318763

私がやりたいこととしては、yを使わず、res$coef[1]とres$coef[2]を使って、lo2とup2を出せないか、ということを考えています。
厳密にlo2とup2に一致しなくとも、正規分布を仮定したときの平均値と標準誤差しかないが、それらから対数正規分布を仮定したときの信頼区間を出す方法、を探しています。
乱数種を替えればyが変わるのでlo1とup1も毎回変わりますが、マイナスのlo1が出る仕組み自体がまずいのでに、それを正の値しか出ない方法に変えたいのです。

どうぞよろしくお願いいたします。


Rによる行列の処理
  投稿者:KUMON 2018/10/24(Wed) 20:13 No. 22626
Rのプログラムに関する質問です。

各成分が0,1 の2値からなるn×m行列があります。
任意の列に対して、i行目が1ならi+1行目も1となるように変換するコードがうまく書けません。ご教授よろしくお願いいたします。

Re: Rによる行列の処理
  投稿者:青木繁伸 2018/10/24(Wed) 20:52 No. 22627
以下のような列があったとき,あなたはどのような結果をお望みですか?

0
1
0
0
1
1
0
0
0
1
1
0
1
0
1
1
0
1
0

Re: Rによる行列の処理
  投稿者:KUMON 2018/10/24(Wed) 21:45 No. 22629
青木先生

ご返信ありがとうございます。

0
1
1
0
1
1
1
0
0
1
1
1
1
1
1
1
1
1
1

という結果です。

Re: Rによる行列の処理
  投稿者:青木繁伸 2018/10/25(Thu) 21:30 No. 22632
以下のようにすることも出来ますが,
a = c(1, 0, 0, 1, 1, 0, 0, 0, 1, 1,  0, 1, 0, 1, 1, 0, 1, 0)
cat(a)
z = (a | c(0, a[-length(a)]))+0
cat(z)


Re: Rによる行列の処理
  投稿者:KUMON 2018/10/30(Tue) 21:07 No. 22633
青木先生

ご回答ありがとうございます。

列ベクトルa の要素を1つずらしたb を使用して(頭は0で最後尾は削除したもの)
ifelse(a+b>0,1,0)と自解していましたが、

先生のコード (a|b)+0 がよくわかりません。解説していただけないでしょうか?


Re: Rによる行列の処理
  投稿者:KUMON 2018/10/30(Tue) 21:25 No. 22634
分かりました!

論理和を使用しているのですね。+0を加えることでTRUE→1、FALSE→0と変換するわけですね。

参考になります。処理スピードが格段に上がりました。ありがとうございました。


相関係数のp値
  投稿者:コロン 2018/10/25(Thu) 17:16 No. 22630
いつもお世話になっております。

毎回稚拙な質問で恐縮なのですが、以下のデータの相関係数行列はもとまるのですが、それぞれのp値を行列表として求めてくれる関数はないのでしょうか。頑張って探したのですが...。

申し訳ございませんが、よろしくお願いいたします。

q1 q2 q3 q4 q5
3 4 2 5 5
1 2 2 4 5
5 4 3 4 1
4 3 5 2 2
1 2 2 1 3
3 3 3 3 3
4 5 5 4 5

Re: 相関係数のp値
  投稿者:コロン 2018/10/25(Thu) 17:57 No. 22631
青木先生

申し訳ございません。自己解決いたしました。先生のサイトの「相関係数行列の計算と検定」を発見いたしました。


Eviewsにおけるイベントスタディ分析について
  投稿者:ゆう 2018/10/23(Tue) 16:24 No. 22625
・題名の通り、Eviewsを使ってどのような手法でイベントスタディ分析を行うのか手順を教えて頂きたいです。

補足
・回帰分析までは出来ます。そこからAR.CARの出し方がわかりません。
どなたかご存知の方よろしくお願いいたします。苦しんでいます。。。

(エクセルのやり方でも大丈夫です。)


複合アウトカム
  投稿者:あおいひげ 2018/10/17(Wed) 17:43 No. 22623
青木先生

ご回答有り難うございます。
正準相関分析・・・参考になりました。
Rを使用するのはかなりハードルが高いと思っています。
再度、研究デザインから検討し直します。

このスピードで、知りたいことに端的にお返事頂いていることに大変感銘を受けました。今後も、お世話になるかと思います。宜しくお願いします。

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/17(Wed) 22:39 No. 22624
今度は,スレッドを維持するように


複合アウトカム
  投稿者:あおいひげ 2018/10/16(Tue) 16:30 No. 22621
青木先生
回答有り難うございます。

複合アウトカムについては、死亡と再入院など、意味合いが異なるイベントを組み合わせることの意義について見直されてきているようです。
 現実的には、列挙したアウトカムの中でも1番重要視しているものをプライマリーアウトカムとし、それ以外はセカンダリーアウトカムとしてオマケぐらいの位置づけで解析をやるしかないのかもしれませんね。
 因みに、正準相関分析は私が使用しているEZRでは処理できそうにないのですが、他のソフトなら可能なものなのでしょうか。
 

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/16(Tue) 17:55 No. 22622
EZR は,R のシュガーコートなので,R の全てに対応しているわけではない
(EZR にどっぷりというのは避けた方がよい)

R には,cancor という関数がある(stats に含まれているので,library とかしないで,すぐ使える)
使用例(example)も見ることができる
> example(cancor)

cancor> ## No test:
cancor> ## signs of results are random
cancor> pop <- LifeCycleSavings[, 2:3]

cancor> oec <- LifeCycleSavings[, -(2:3)]

cancor> cancor(pop, oec)
$cor
[1] 0.8247966 0.3652762

$xcoef
[,1] [,2]
pop15 -0.009110856 -0.03622206
pop75 0.048647514 -0.26031158

$ycoef
[,1] [,2] [,3]
sr 0.0084710221 3.337936e-02 -5.157130e-03
dpi 0.0001307398 -7.588232e-05 4.543705e-06

以下略


複合アウトカム
  投稿者:あおいひげ 2018/10/15(Mon) 18:15 No. 22619
初めて投稿致します。

医学系臨床研究のデザインを検討しています。統計初心者で、いろいろ試行錯誤で苦しんでおりますのでご教授お願い致します。統計ソフトはEZRを使用しております。

集団をある値が高い群(H群)と低い群(L群)の2群に分けました。100〜200例 VS 300〜500例 (集計を開始していないので見積もり概算です) 。
 アウトカムとして、連続変数と二値変数の2種類の変数を全部で10個くらい想定しています。交絡因子は年齢、性別、既往歴、など10種類くらいを検討しています。
 多変量解析になるかと思いますが、複合アウトカムの処理方法が分かりません。アウトカムを一つずつ、重回帰やロジスティックなどを多重性に気を付けながらやるしかないのでしょうか。

曖昧な質問ですみません。
お教え下さいませ。

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/15(Mon) 19:21 No. 22620
結果が複数の変数ということならば,正準相関分析(重回帰分析の一般形)ということになるかな?実際に使った研究例は多くはないかも。


【R】 リスト構造からの読み取り
  投稿者:明石 2018/10/02(Tue) 21:46 No. 22616
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

---------------------------

ベイジアンネット bnlearn パッケージの出力結果のリスト構造からグラフ構造を取得したいと思います。
以下、グラフ描画までのコード例をお示しします。

library(bnlearn) # ベイジアンネットのパッケージ

df <- learning.test # 組み込みのデータセット

res <- hc(df)   # 構造学習 hill-climbing法を適用
plot(res, main = "hill-climbing法による構造学習") # 添付の画像ファイル

str(res)

取得したリスト構造に注釈 銑┐鯢嬪燭靴燭發里髻∨尾にお示しをします。
注釈 銑┐硫媾蠅ら、グラフの親子構造(親は→の元、子は→の先)を出力する
Rプログラムのご教示をいただければ大変に助かります。

誠にお恥ずかしいのですが、まったく手も足も出ません。

所望する出力です。
A,B
A,D
B,E
C,D
F,E

--------------------------------
以下、str(res) の出力結果に、 銑┐涼躰瓩鯢嬪燭靴泙靴拭

List of 3
$ learning:List of 8
..$ whitelist: NULL
..$ blacklist: NULL
..$ test : chr "bic"
..$ ntests : num 40
..$ algo : chr "hc"
..$ args :List of 1
.. ..$ k: num 4.26
..$ optimized: logi TRUE
..$ illegal : NULL
$ nodes :List of 6
..$ A:List of 4             ⇒
.. ..$ mb : chr [1:3] "B" "C" "D"
.. ..$ nbr : chr [1:2] "B" "D"
.. ..$ parents : chr(0)
.. ..$ children: chr [1:2] "B" "D"   ⇒
..$ B:List of 4             ⇒
.. ..$ mb : chr [1:3] "A" "E" "F"
.. ..$ nbr : chr [1:2] "A" "E"
.. ..$ parents : chr "A"
.. ..$ children: chr "E"        ⇒
..$ C:List of 4              ⇒ァ      
.. ..$ mb : chr [1:2] "A" "D"
.. ..$ nbr : chr "D"
.. ..$ parents : chr(0)
.. ..$ children: chr "D"        ⇒
..$ D:List of 4
.. ..$ mb : chr [1:2] "A" "C"
.. ..$ nbr : chr [1:2] "A" "C"
.. ..$ parents : chr [1:2] "A" "C"
.. ..$ children: chr(0)
..$ E:List of 4
.. ..$ mb : chr [1:2] "B" "F"
.. ..$ nbr : chr [1:2] "B" "F"
.. ..$ parents : chr [1:2] "B" "F"
.. ..$ children: chr(0)
..$ F:List of 4              ⇒
.. ..$ mb : chr [1:2] "B" "E"
.. ..$ nbr : chr "E"
.. ..$ parents : chr(0)
.. ..$ children: chr "E"        ⇒
$ arcs : chr [1:5, 1:2] "A" "A" "C" "B" ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "from" "to"
- attr(*, "class")= chr "bn"

何卒どうぞ、よろしくお願いいたします。
いつもありがとうございます。


Re: 【R】 リスト構造からの読み取り
  投稿者:青木繁伸 2018/10/02(Tue) 23:14 No. 22617
この件については,何も知りませんが,
以下のようにすれば,簡単に解が得られのでは?
> res$arcs
from to
[1,] "A" "B"
[2,] "A" "D"
[3,] "C" "D"
[4,] "B" "E"
[5,] "F" "E"
See bn-class for details.

arcs: the arcs of the Bayesian network (a two-column matrix, whose columns are labeled from and to). Undirected arcs are stored as two directed arcs with
opposite directions between the corresponding incident nodes.

Re: 【R】 リスト構造からの読み取り
  投稿者:明石 2018/10/03(Wed) 19:59 No. 22618
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

今回も大変に良い勉強をさせていただきました。
心から御礼を申し上げます。


ゼロでの割り算
  投稿者:青木繁伸 2018/09/09(Sun) 11:10 No. 22612
中澤さんが,2x2 より大きい行列で R の mcnemar.test の結果が NaN になる場合があることについて,書いている。

http://minato.sip21c.org/bulbul2/20180908.html
http://minato.sip21c.org/bulbul2/20160913.html
http://minato.sip21c.org/bulbul2/20160911.html

「対角成分の和がゼロの組合せを計算からスキップできる拡張マクネマー検定の関数」として
mct <- function(x) {
L <- NROW(x)
X <- 0
N <- 0
for (i in 1:(L - 1)) {
for (j in (i + 1):L) {
s <- (x[i, j] + x[j, i])
if (s > 0) {
N <- N + 1
X <- X + (x[i, j] - x[j, i]) ^ 2 / s
}
}
}
return(c( X2 = X, df = N, p = (1 - pchisq(X, N))))
}
を挙げている。
しかし,「このやり方が正しいという理論的根拠を見つけられないので悩ましいところだが」とあるのだが。

「N23とN32がゼロ,という状態だと,(N23+N32)が分母になる項で"Division by zero"エラーがでてしまって,カイ二乗値の計算結果自体がNaNになってしまう」そこで「分母がゼロになる組合せをスキップする」ということを書いている。

この記述は正しくない。N23=N32=0のとき,つまり(0-0)/(0+0) = 0/0 の計算結果は,NaN である("Division by zero"ではない。Python だと,ZeroDivisionError: division by zero になる)。何かに NaN を足すと NaN になり,それ以降 NaN に何を足しても NaN のまま。ということで,R の mcnemar.test は NaN を返すと,中澤さんも書いている。

N23=N32=0 のときの計算結果が NaN になった場合は足し算をスキップすればよいわけで,N23=N32=0 の場合は計算をスキップしても同じではある。またその場合には自由度のカウントはしない。

以下のような関数を書いてみた。
mct2 <- function(x) {
i <- col(x)[lower.tri(x)]
j <- row(x)[lower.tri(x)]
a <- mapply(function(i, j) (x[i, j] - x[j, i]) ^ 2 / (x[i, j] + x[j, i]), i, j)
chisq <- sum(a, na.rm = TRUE) # NaN も除かれる
df <- sum(!is.nan(a)) # !is.na(a) でもよい
return(c(X2 = chisq, df = df, p.value = pchisq(chisq, df, lower.tail = FALSE)))
}
「ゼロ割」は避けるのではなく,積極的に正しく対処するべきということか?(同じようなことではあるが)

蛇足ではあるが,自由度を減ずるのは,カイ二乗分布の再生性から明らか(自由度 1 のカイ二乗分布にしたがう確率変数 k 個の和は自由度 k のカイ二乗分布にしたがう)。

R の mcnemar.test を修正するとすれば,
STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)])

STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
PARAMETER <- sum(!is.nan(y[upper.tri(x)]^2/x[upper.tri(x)]))
にする(前の方の PARAMETER <- r * (r - 1)/2 は不要)。
プログラムにおいて,特殊な条件の発生を見逃すと命取りというところ。
しかしさすがに,sum(y[upper.tri(x)]^2/x[upper.tri(x)]) とは,にくいねえ。
> x = matrix(c(c(7, 2, 5, 3, 11, 2, 0, 3, 4, 0, 5, 9, 3, 15, 14, 10)), 4, 4)
> mct(x)
X2 df p
15.428836864 5.000000000 0.008678852
> mct2(x)
X2 df p.value
15.428836864 5.000000000 0.008678852

> mcnemar.test(x) # 元のまま

McNemar's Chi-squared test

data: x
McNemar's chi-squared = NaN, df = 6, p-value = NA

> mcnemar.test2(x) # 修正したもの

McNemar's Chi-squared test

data: x
McNemar's chi-squared = 15.429, df = 5, p-value = 0.008679
なお,「"Division by zero"を避けるために全部のセルに0.01を足す」というのは明らかに間違い(自由度が減らない)。

Re: ゼロでの割り算
  投稿者:青木繁伸 2018/09/09(Sun) 12:49 No. 22613
2 0
0 5

というのと

2 1
1 5

というのは,0-0 = 1-1 = 0 という情報を与えるという意味では同じなので,自由度は減らさなくてもよいか

STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
だけの変更で,これは実質的には NaN の場合に 0 を加えるというのと同じであるから
PARAMETER はそのまま PARAMETER <- r * (r - 1)/2

とすれば,mcnemar.test(x+0.01) もありか。ただし,加える数は他のまともなセルのことも考えるとできるだけ小さい方が良いだろうとは思う。例えば +1e-15 とか

Re: ゼロでの割り算
  投稿者:中澤 2018/09/09(Sun) 21:30 No. 22614
記事ありがとうございます。
確かにNaNが正しいですね。不正確な記述でした。
自由度の扱いは悩ましいですね。スキップでなくて分子がゼロなのでゼロとみなすということですね。勉強になりました。

Re: ゼロでの割り算
  投稿者:中澤 2018/09/19(Wed) 17:55 No. 22615
http://minato.sip21c.org/bulbul2/20180915.html
に書きましたが,カテゴリ数が3以上の場合の拡張マクネマー検定と同様の帰無仮説を検定する方法としてBhapkarの検定というのがあって,irrパッケージのbhapkar()で実行できました。いま書いている論文ではこれを使おうと思います。


パネルデータの欠損値補完について
  投稿者:羽島 2018/09/04(Tue) 09:21 No. 22609
青木先生

羽島と申します、この度はお世話になります。

Rを使って以下のようなデータを分析しており、その中で住所コードの欠損値を補完する必要が出てきました。LOCFのような手法での補完では精度に不安があるため、個人IDを参照することで欠損値を補完したいのですが、どのようにプログラムを書けば良いでしょうか?

お忙しいところ大変恐縮ですが、ご教示頂けましたら幸いです。

個人ID  年   住所コード

100   2015    810
100   2016    NA
100   2017    810
201   2015    NA
201   2016    NA
201   2017    666

Re: パネルデータの欠損値補完について
  投稿者:青木繁伸 2018/09/04(Tue) 13:03 No. 22610
あまり簡単なプログラムは書けないかな
d <- data.frame(ID = c(100, 100, 100, 201, 201, 201, 330, 330, 330),
year = c(2015, 2016, 1017, 2015, 2016, 2017, 2015, 2016, 2017),
code = c(810, NA, 810, NA, NA, 666, NA, 543, NA))
tbl <- unique(na.omit(d[c("ID", "code")]))
for (i in seq_len(nrow(d))) {
if (is.na(d$code[i])) {
d$code[i] = tbl[which(d$ID[i] == tbl[,1]), 2]
}
}
d

Re: パネルデータの欠損値補完について
  投稿者:羽島 2018/09/04(Tue) 18:15 No. 22611
青木先生

お返事の方ありがとうございます。

元のデータは3万件ほどの個人データがあり手作業での修正を覚悟しておりましたが、
先生のお力で無事補完することが出来ました。

自身の勉強不足を恥じるばかりですが、この度は本当にありがとうございました。


【R】 二重繰り返し処理の高速化
  投稿者:明石 2018/09/01(Sat) 10:42 No. 22606
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

---------------------------

記事 No.22590 「繰り返し処理の高速化」では、ベクトル化を用いた、劇的な高速化のご教示をいただきました。
ありがとうございました。

今回は、二重繰り返し処理です。

例示データは、先生のサイトからお借りいたします。

http://aoki2.si.gunma-u.ac.jp/R/my-cor.html
単相関係数,偏相関係数,重相関係数を計算する

x <- matrix(c( # 5ケース,4変数のデータ行列例(ファイルから読んでも良い)
1, 5, 6, 4,
2, 14, 5, 3,
3, 3, 4, 2,
4, 2, 6, 6,
3, 4, 3, 5
), ncol=4, byrow=TRUE)

ここでは、相関係数行列を計算していますが、
今回、私がやりたいのは、ケース(行ベクトル)同士のコサイン類似度行列です。

以下、プログラムを作成しました。

x <- matrix(c( # 5ケース,4変数のデータ行列例(ファイルから読んでも良い)
1, 5, 6, 4,
2, 14, 5, 3,
3, 3, 4, 2,
4, 2, 6, 6,
3, 4, 3, 5
), ncol=4, byrow=TRUE)

nr <- nrow(x)
mx <- matrix(0,ncol=nr,nrow=nr)
diag(mx) <- 1

for(i in 1:(nr-1)) {
for(j in (i+1):nr) {
s <- (x[i,] %*% x[j,]) / (sqrt(sum(x[i,]^2)) * sqrt(sum(x[j,]^2)))
mx[i,j] <- s
mx[j,i] <- s
}
}

> mx
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0000000 0.8438196 0.9183979 0.8735555 0.8992005
[2,] 0.8438196 1.0000000 0.7847512 0.5725026 0.7829858
[3,] 0.9183979 0.7847512 1.0000000 0.9132886 0.9081355
[4,] 0.8735555 0.5725026 0.9132886 1.0000000 0.9229730
[5,] 0.8992005 0.7829858 0.9081355 0.9229730 1.0000000

私が扱いたいデータでは、
ケース数は、数千
属性数は、数百
ですので、高速化処理ができれば大変にありがたいです。

ご教示をいただければ、大変に助かります。
どうぞよろしくお願いいたします。

Re: 【R】 二重繰り返し処理の高速化
  投稿者:青木繁伸 2018/09/02(Sun) 11:31 No. 22607
二重 for-loop の内側で使われる
x[i,] %*% x[j,]
は,行ベクトルの積和ですから
for-loop なしの
z <- x%*%t(x)
です。

同じく,分母は
sqrt(sum(x[i,]^2))という積和の平方根ですから x%*%t(x) の対角要素の平方根です
d <- sqrt(diag(z))
二重 for-loop で使われる
sqrt(sum(x[i,]^2)) * sqrt(sum(x[j,]^2))
は,
d2 <- outer(d, d)
です。
以上をまとめて,ループなしで答えが求まります。
z <- x%*%t(x)
d <- sqrt(diag(z))
d2 <- outer(d, d)
mx2 <- z / d2
5000x500の行列だと,私のマシンではfor-loop プログラムは約400秒,ベクトル演算を使ったプログラムは約10秒で,およそ40倍速くなります

Special Thanks !(Re: 【R】 二重繰り返し処理の高速化)
  投稿者:明石 2018/09/02(Sun) 12:24 No. 22608
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

劇的な改善方法をご教示をいただき、誠にありがとうございます。

表層的ではなく、計算の本質を見抜くことの重要性を教えていただきました。

心から、心より、御礼を申し上げます。
本当にありがとうございました。


適合度のexact testの結果の書き方について
  投稿者:黒岩 2018/08/24(Fri) 19:13 No. 22603
青木先生

いつもお世話になっております。
下記についてご教示頂けると幸いです。

gftで適合度のexact testを行うと、
カイ二乗値、自由度、P値、正確なP値、が出力されますが、
レポートなどでこの結果を書くときの
適切な記述方法を教えて頂けないでしょうか?

単純にこの4つを書き並べれば良いのでしょうか?
P値をはっきり区別するために、
例えば漸近有意確率P、正確有意確率P、などと
書き分ける必要はないのでしょうか?
あるいは、正確なP値のみではいけないのでしょうか?

Re: 適合度のexact testの結果の書き方について
  投稿者:青木繁伸 2018/08/24(Fri) 21:39 No. 22604
「適合度の正確検定による有意確率 p < 0.xxx」でよいと思いますが...

Re: 適合度のexact testの結果の書き方について
  投稿者:黒岩 2018/08/25(Sat) 04:32 No. 22605
青木先生

早速の回答、ありがとうございました!


NaNを1に置き換えたい
  投稿者:赤岳 2018/08/21(Tue) 17:30 No. 22600
いつもお世話になっております。
次の文をRで指定しましたが、logの値がマイナスになり、エラーとなります。
logの中の乱数がマイナスやゼロの場合、1を返すようにしたいのですが、うまく文をかけません。
ご教示いただけないでしょうか。
*****************
> N <-100
> P<- rnorm(N, 5, 1)
> dP <- c(1:N)*0
> EP <- c(1:N)*0
> count<- 0
> for (i in 1:N) {
+ dP[i] <- P[i] - log(rnorm (1, 1, 0.5))
+ EP[i] <- exp (dP[i])/exp(P[i])
+ if (EP[i] >0.5)
+ count <- count+1
+ }
log(rnorm(1, 1, 0.5)) で警告がありました: 計算結果が NaN になりました
if (EP[i] > 0.5) count <- count + 1 でエラー:
TRUE/FALSE が必要なところが欠損値です
> count/N
[1] 0.41

Re: NaNを1に置き換えたい
  投稿者:青木繁伸 2018/08/21(Tue) 22:23 No. 22601
本筋ではないところをコメント
dP <- c(1:N)*0 は dp <- numeric(N) のほうが適切でしょう。

rnorm(1, 1, 0.5) が 0 以下のとき,log(rnorm(1, 1, 0.5)) を 1 にするのですか?

for ループを使っている場合は,特に何の工夫も必要ないです
N <-100
P <- rnorm(N, 5, 1)
dP <- c(1:N)*0
EP <- c(1:N)*0
count<- 0
for (i in 1:N) {
x <- rnorm (1, 1, 0.5) # ここから
if (x <= 0) {
x <- 1
} else {
x <- log(x)
}
dP[i] <- P[i] - x # ここまで追加(変更)
EP[i] <- exp (dP[i])/exp(P[i])
if (EP[i] >0.5)
count <- count+1
}
count/N


for ループをなくしてみましょう。ベクトル演算します。
N <- 100
p <- rnorm(N, 5, 1)
x <- rnorm(N, 1, 0.5) # まとめて N 個生成
x <- ifelse(x <= 0, exp(1), x) # x <= 0 のとき log(x) を 1 にするなら,x = exp(1) にしておく
dP <- P - log(x) # x が正なら log(x),正でなければ log(exp(1)) = 1 になる
EP <- exp(dP) / exp(P) # EP の計算
mean(EP > 0.5) # count = sum(EP > 0.5); count/N は,mean(EP > 0.5)


注: ifelse(x <= 0, 1, log(x)) としても,一度全ての x に log(x) を適用して,x <= 0 のときに 1,そうでないときに log(x) の値を取るように動作するので,エラーが出る。
そこで,上のようにへんな(?)ことをする。

Re: NaNを1に置き換えたい
  投稿者:赤岳 2018/08/21(Tue) 22:53 No. 22602
青木先生
ご教示ありがとうございました。
> mean(EP > 0.5) # count = sum(EP > 0.5); count/N は,mean(EP > 0.5)
このmeanがこんな風に使われるのは初めて知りました。
また、確かにifelse(x <= 0, 1, log(x)) のようにしましたが、エラーメッセージはクリアーできませんでした。
解決しました。ありがとうございました。
(p <- rnorm(N, 5, 1)は、P <- rnorm(N, 5, 1)と修正しました。)


直接法と群間差の検定
  投稿者:M 2018/08/18(Sat) 21:48 No. 22593 HomePage
青木先生  皆様

お世話になります。いろいろ調べたのですが、自信がないので、質問させていただきたいと存じます。

添付のサイトのP182 表88のような検定を行いたいのですが、検定方法がわかりません。

質問1: 直接法とは何でしょうか?

表88には、注釈に、注1)年齢階級(20−29歳,30−39歳,40−49歳,50−59歳,60−69歳,70歳以上の6区分)と世帯員数(1人,2人,3人以上世帯の3区分)での調整値。割合に関する項目は直接法,平均 値に関する項目は共分散分析を用いて算出。 む

と書いてあります。

また、群間差の検定について、注2)世帯の所得について,多変量解析(世帯の所得額を当該世帯員に当てはめて,割合に関する項目はロジスティック回帰分析,平均値に関する項目は共分散分析)を用いて600万 円以上を基準とした他の2群との群間比較を実施。 とあり、

質問2: 群間差の検定は、200万円未満VS600万円以上、200-600未満VS600万円以上の2回、群間差の検定を行えばよろしいでしょうか?T検定ではご法度なので、躊躇しております。  

また、報告書の冒頭には、
(1)所得と生活習慣・食生活に関する分析  年齢(20−29歳,30−39歳,40−49歳,50−59歳,60−69歳,70歳以上 の6区分)と世帯人数(1人,2人,3人以上世帯の3区分)で調整を行った後, 割合に関する項目は多変量ロジスティック回帰分析,平均値に関する項目は共分 散分析を行った

群間比較について,年齢調整を行い分析を行ったものについては,割合に関す る項目はCochran-Mantel-Haenszel検定,平均値に関する項目は,共分散分析を 行った。

とあります。

お手数をお掛け致しますが、どうぞよろしくお願い致します。

Re: 直接法と群間差の検定
  投稿者:青木繁伸 2018/08/19(Sun) 06:12 No. 22594
サイト情報がありません(家マークでのリンクではなく,本文中に記載してください。)

というわけで,詳細は分かりませんが,直接法については以下のごとく。

直接法とは年齢調整死亡率の計算に用いられる方法。ほかに,間接法もある。
死亡率に限らず,率の調整に使われることもある。
直接法についての説明は,
http://www.takenet.or.jp/~hayakawa/u-tanqa04.htm
http://www.kenko.pref.yamaguchi.lg.jp/img/kenko21/download/h26map/h26kaisetsu-new.pdf
など,
標準化死亡比の検定は,
https://www.niph.go.jp/soshiki/07shougai/datakatsuyou/data/h22smr-kinki.pdf
などを参照のこと。

ここで質問しても詳細な(適切な)回答が得られるとは限りません。
著者に直接お問い合わせになるのが筋でしょう。

Re: 直接法と群間差の検定
  投稿者:M 2018/08/19(Sun) 14:43 No. 22596
青木先生 

お返事、どうもありがとうございました。
また、リンクが添付できず大変失礼いたしました。

リンクは、
https://www.mhlw.go.jp/bunya/kenkou/eiyou/dl/h26-houkoku.pdf

ページ 182 表88 です。

特に、群間差の検定を行いたいのです。この表88では、運動習慣のない者(%)と表示してあり、P値で群間の検定結果が表示されています。
どうやって検定したものか、知りたいのです。

表の注ではなく、報告書の冒頭に「割合に関す る項目はCochran-Mantel-Haenszel検定」
とありますので、Cochran-Mantel-Haenszel検定かなと思うのですが、そうするとすべて2値にしなければならないでしょうか?

例えば 運動習慣 ある・ない× 性別 男・女 

層として、収入200万円未満・600万円以上、
     200-600未満・600万円以上
     200万円未満・200-600未満

の3回、Cochran-Mantel-Haenszel検定を行い、有意確率で、層に設定した収入の違いに関係なく有意に差があるかを検定すればよろしいでしょうか?

T検定では、こういう組み合わせにした検定はご法度のはずで、躊躇しております。

うまく説明できたか不安なのですが、厚労省が回答してくれるかわかりませんので、再質問させていただきました。お手数をお掛け致しますが、どうぞよろしくお願い致します。

Re: 直接法と群間差の検定
  投稿者:青木繁伸 2018/08/19(Sun) 18:32 No. 22597
厚生省は発注元?で,実際は「立研究開発法人医薬基盤・健康・栄養研究所は,厚生労働省に提出された調査票について,入力・集計・作表を行った。」わけで,研究者レベルの質問には答えてくれるかも。
主要な成果については論文として投稿される(されている)かも。この年のデータのものでなくても,以前のものにも同じ分析手法が使われていると思うので,もう少し丁寧な説明や参考文献も記載されているかも。

さて,182ページ表88ですが,

「注2)世帯の所得について,多変量解析(世帯の所得額を当該世帯員に当てはめて,割合に関する項目はロジスティック回帰分析,平均値に関する項目は共分散分析)を用いて600万円以上を基準とした他の2群との群間比較を実施。」
ですから,ロジスティック回帰分析による結果でしょう(詳細はわかりません)。

「割合に関する項目は Cochran-Mantel-Haenszel検定」というのは,たとえば,50ページの表8の脚注に同じ記述がありますよね。層は年齢階級なのかな(よくわからない)。

統計解析過程についての詳細な説明がないので,よくわかりません。

Re: 直接法と群間差の検定
  投稿者:青木繁伸 2018/08/19(Sun) 20:54 No. 22598
解析過程をもっとも正確に記述できるのはプログラミング言語なので,全部とはいわず代表的なデータ解析を行ったプログラムを掲載するのは優れた試みだと思いますね。
R なんかは,最適でしょう。まさか Excel でやってますということはないと思いますけど。

Re: 直接法と群間差の検定
  投稿者:M 2018/08/20(Mon) 13:02 No. 22599
青木先生

お返事、どうもありがとうございました。

何度読んでもよくわからないのは、解析過程がわかるよう書かれていないからのようですね。

青木先生がおっしゃるようにデータ解析の過程をある程度、読者が理解できるように書く方がデータとしても大切だと改めて思いました。

貴重なお時間をどうもありがとうございました。心から感謝です。


【R】 繰り返し処理の高速化
  投稿者:明石 2018/08/14(Tue) 18:06 No. 22590
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

---------------------------
受験者の生年月日と試験日から、試験日での年齢を計算したいと思います。

受験者データはデータフレームdatです。
生年月日の変数名はbirthです。日付型です。

試験日の変数名はtestです。日付型です。
test <- as.Date("2018-9-13")

以下の繰り返し処理で、受験者の試験日での年齢yearを計算できることを確認しました。
for (i in 1:nrow(dat)) {
dat$year[i] <- length(seq(dat$birth[i], test, by = "year")) - 1
}

受験者データが大規模ですので、繰り返し処理を高速化したいと思います。

ご教示いただければ大変に助かります。
何卒どうぞ、よろしくお願いいたします。

Re: 【R】 繰り返し処理の高速化
  投稿者:青木繁伸 2018/08/15(Wed) 02:12 No. 22591
50000 人分のテストデータを作ってやってみました
test <- as.Date("2018-9-13")
dat <- data.frame(birth=as.Date(Sys.Date()-1:50000)) # 50000 人分のテストデータ

system.time({
for (i in 1:nrow(dat)) {
dat$year[i] <- length(seq(dat$birth[i], test, by = "year")) - 1
}
})
ユーザ システム 経過
11.711 6.689 18.254

年齢の計算は,簡単なのです。試験年から生まれ年を引くだけ。
ただし,「生まれ月が試験月より後(大きい)なら」,
または,「生まれ月が試験月と同じで,生まれ日が試験日より後(大きい)なら」
1を引く。

わざわざ as.Date となっているけど,そこから年月日(スカラー y, m, d)を取り出す。
試験年月日は yyyy, mm, dd の整数で与える(as.Date から取りだしてもよいけど)

age2 <- function(x, yyyy = 2018, mm = 09, dd = 13) {
y <- as.integer(substr(x, 1, 4))
m <- as.integer(substr(x, 6, 7))
d <- as.integer(substr(x, 9, 10))
val <- yyyy - y
if (m > mm || (m == mm && d > dd)) {
val <- val-1
}
val
}
system.time({
for (i in 1:nrow(dat)) {
dat$year2[i] <- age2(dat$birth[i])
}
})

でも,これでは for を含むし,年月日の取り出しコストが掛かるのでかえって遅い。

ユーザ システム 経過
13.662 6.667 20.312

一応答えは合っているようだ。

> all(dat$year == dat$year2)
[1] TRUE

for を避けるためにベクトル演算を行うようにする。
関数には dat$birth ベクトルを渡す。
年月日の取り出し,それに基づく年齢(val)の計算も全部がベクトル演算。

age3 <- function(birth, yyyy = 2018, mm = 09, dd = 13) {
y <- as.integer(substr(birth, 1, 4))
m <- as.integer(substr(birth, 6, 7))
d <- as.integer(substr(birth, 9, 10))
val <- yyyy - y
ifelse(m > mm | (m == mm & d > dd), val-1, val)
}
system.time(dat$year3 <- age3(dat$birth))

27 倍ほど速くなった。

ユーザ システム 経過
0.424 0.023 0.448

答えも合っているようだ。

> all(dat$year == dat$year3)
[1] TRUE

なお,val を計算している 2 行を(ifelse を使うのが気持ち悪いからといって)

yyyy - y - (m > mm | (m == mm & d > dd))

にしても,ほとんど速度に違いはない。

年月日が3つのベクトルで与えられれば as.Date から取り出す必要がないのでもっと速い。
birth を "2018-08-15" とわざわざ入力しているのか,2018, 8, 15 という3つの整数を入力しているのかによる。
後者なら,その3つの整数から文字列を作って as.Date していることになる。
3つの整数値を入力してそれをそのまま使うなら,人間もコンピュータも無駄なことをしなくてもよいということになる。

御礼(Re: 【R】 繰り返し処理の高速化)
  投稿者:明石 2018/08/15(Wed) 08:04 No. 22592
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

多面的にご教示いただき、
特に、for を避けるためにベクトル演算を行うようにするところは、
大変に勉強になります。

毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
心から御礼を申し上げます。
ありがとうございました。


Rの代入演算子について
  投稿者:波音 2018/08/09(Thu) 12:18 No. 22587
素朴な疑問ですが、ここのところ青木先生が代入演算子を

<-

ではなく、

=

と書かれるようになったのには、何か意味があるのでしょうか?

※書名は忘れましたが = ではなく <- の方がいいと書いてあった記憶があります。

最近になってソースコード表現の定石も変わってきているのでしょうか、、、

Re: Rの代入演算子について
  投稿者:青木繁伸 2018/08/09(Thu) 15:00 No. 22588
ソースコード表現の定石は変わっていないと思います
「GoogleのRスタイルガイドは、割り当てのために「=」を禁止する」ということもあり,
? "<-" で表示される R のマニュアルでも,The operator <- can be used anywhere, whereas the operator = is only allowed at the top level (e.g., in the complete expression typed at the command prompt) or as one of the subexpressions in a braced list of expressions. とありますので,<- を使うのが妥当であることに変わりありません。

= を使うようになったのは些細な理由からですが,= を使うようになっての感想は,何の不便もなく,タイプ数も(わずか1文字の違いですが)少なく,他の言語を使ってきた経験からも = の方が違和感がないなあというだけです。

公の場では <- を使うようにしたいと思います。

Re: Rの代入演算子について
  投稿者:波音 2018/08/10(Fri) 12:01 No. 22589
回答ありがとうございます。

私はRが経験年数的には一番長いので <- に慣れてしまっていますが、多言語習得ユーザーとしては(明示的な問題が生じない限り) = でも不都合はなさそうと思います。

※代入演算子とは別に、最近のR使いの方はパイプを多用した記述をされるので、コードを読むのに一苦労することが多いです。


Dataを部分的に取り出す方法
  投稿者:赤岳 2018/08/07(Tue) 22:31 No. 22584
いつも勉強させていただいています。
Rの基本操作だと思いますが、わかりません。
今、次のようなデータが1000個あるとします。
x<- rnorm(1000, 15.5, 3.3)

このデータの95%区間にあるデータを取り出したいのですが、どのように取り出せばよいのでしょうか。
たとえば、
> y<- quantile(x, c(0.025, 0.975))
とすると、yのlength=2となり、950個のデータが取り出せません。
どなたか、yが950個のベクトルになるようなコマンドを教えてもらますでしょうか。

Re: Dataを部分的に取り出す方法
  投稿者:青木繁伸 2018/08/08(Wed) 05:48 No. 22585
x <- rnorm(1000, 15.5, 3.3) のベクトル x から,10 以上, 15 以下のデータを取り出すときはどのようにしますか?

x[x >= 10 & x <= 15] ですね

では y <- quantile(x, c(0.025, 0.975)) から,下限と上限を取り出すのはどのようにしますか?

y[1], y[2] ですね

それをあわせれば

x[x >= y[1] & x <= y[2]] です

Re: Dataを部分的に取り出す方法
  投稿者:赤岳 2018/08/08(Wed) 08:26 No. 22586
青木先生、
ありがとうございました。
ここ2日くらいそのような関数というかコマンドばかり探していました。
今後ともよろしくお願いします。


ノンパラ、パラメトリックの決め方
  投稿者:コロン 2018/07/27(Fri) 07:51 No. 22581
いつもお世話になっております。

基本的なことで恐縮なのですが、ご教示いただけますか。

手元に、テスト得点データがあります。実験群、統制群の2グループで、サンプルサイズはそれぞれ、30、25です。平均点の差の検定を行いたい場合は、普通はt検定を行うのではと考えます。

これを、サンプルサイズが小さいという理由で、ノンパラメトリック検定で行うことは大丈夫でしょうか。なぜこのような質問をしたかと言いますと、ある論文の査読をしているのですが、サンプルサイズが小さいとありますが、GPowerで必要な被験者数を求めると上のサイズは小さいとは言えず、テスト得点ですので、母分布は正規分布をなしていると考えます。また、本論文では、正規分布をなしていない根拠として、手元のデータに対して正規分布をなしているかどうかの検定を行い、正規分布をなしていないとも言っています。

そもそもとしまして、ノンパラにするかどうかは手元のデータが正規分布をなすかどうかで決めるものなのか、以前、この掲示板のどこかで指摘ぐあったような記憶があります。

どういう視点で、ノンパラか否かを決めるべきか、また、決定のプロセスの中でやってはいけないこと、ぜひご教示ください。

Re: ノンパラ、パラメトリックの決め方
  投稿者:青木繁伸 2018/07/28(Sat) 20:49 No. 22582
> サンプルサイズが小さいという理由で、ノンパラメトリック検定で行う

あまり説得力はないでしょう
しかし,

> GPowerで必要な被験者数を求めると上のサイズは小さいとは言えず

とはいっても,マン・ホイットニーの U 検定の検定力は,条件を満たしているときの平均値の差の検定(t 検定)の 95% ほどもあるので,劣っている検定手法を採用したと批難するのはちょっと厳しすぎるかも。

> 手元のデータに対して正規分布をなしているかどうかの検定を行い、正規分布をなしていない

検定の多重性の問題が出てくるでしょう
しかし,理論的には正規分布に従うと考えられるものの,標本データでそれを確認できないとすれば,正規性の検定を行っても批難はされないでしょう

ノンパラメトリック検定は種々の仮定を置かないので,ノンパラメトリック検定で有意ならばややこしいことを避けることができるという長所があります。

> 決定のプロセスの中でやってはいけないこと

「データを見てから検定法を選択すること」はやってはいけないことと言えるでしょう。

Re: ノンパラ、パラメトリックの決め方
  投稿者:コロン 2018/07/29(Sun) 08:08 No. 22583
青木先生

とても詳細なご説明をありがとうございます。勉強になりました。

今後ともよろしくお願いいたします。

失礼いたします。


R 半角文字→全角文字への変換
  投稿者:明石 2018/07/25(Wed) 19:23 No. 22578
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

-------------------------

半角、全角、大文字、小文字が混在する文字列を、全角大文字に統一したいのです。

例示ですが、

moji <- "ボーイング787 ボーイング787 boeing787 Boeing787 boeing787 BOEING787"

"ボーイング787 ボーイング787 BOEING787 BOEING787 BOEING787 BOEING787"

大文字に置換は、toupper(moji)を使うえばよいことは分かりましたが、
誠に情けないのですが、半角→全角の関数を見つけることができません。

ご教示いただけましたら大変に助かります。
お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

Re: R 半角文字→全角文字への変換
  投稿者:青木繁伸 2018/07/25(Wed) 19:53 No. 22579
文字の変換には基本的には chartr 関数を使えば良いです
例題の場合には,大文字変換に toupper を使えば,少し簡単に書けます
> moji = "ボーイング787 ボーイング787 boeing787 Boeing787 boeing787 BOEING787"
> moji = toupper(moji)
> chartr("A-Z0-9", "A-Z0-9", moji)
[1] "ボーイング787 ボーイング787 BOEING787 BOEING787 BOEING787 BOEING787"

> moji = "ボーイング787 ボーイング787 boeing787 Boeing787 boeing787 BOEING787"
> chartr("a-zA-Za-z0-9", "A-ZA-ZA-Z0-9", moji)
[1] "ボーイング787 ボーイング787 BOEING787 BOEING787 BOEING787 BOEING787"


Re: R 半角文字→全角文字への変換
  投稿者:明石 2018/07/26(Thu) 09:22 No. 22580
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

テキスト分析で、語のゆらぎを吸収する際に必要な前処理であり、大変に助かりました。

つまらない質問で、ご高名な青木先生に失礼にならないか、
小心者ですので、躊躇していましたが、
勇気をふりしぼってご相談して良かったと思います。

毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。


対応のある2標本の母分布の違いの分析について
  投稿者:HG 2018/06/21(Thu) 12:51 No. 22574
対応のある2標本の母分布に違いがあるかを検証したいのですが、
どのような分析手法を採用すべきでしょうか。

2標本コルモゴロフ・スミルノフ検定は独立2群の比較なので
この場合は使えないと考えております。

ご教示のほど、よろしくお願いいたします。

Re: 対応のある2標本の母分布の違いの分析について
  投稿者:青木繁伸 2018/06/21(Thu) 21:07 No. 22575
「母分布に違い」というのは,位置の母数の違いなのか,散布度の違いなのか?それとも,両方含めての違いなのか?

統計手法の選択ガイド
http://aoki2.si.gunma-u.ac.jp/FlowChart/Tutorial.html
をやってみるとどうなりました?

Re: 対応のある2標本の母分布の違いの分析について
  投稿者:HG 2018/06/22(Fri) 13:22 No. 22577
青木先生
コメントありがとうございます。

「母分布に違い」とは散布度(分布の形)の違いを意図しておりました。説明不足であり、
申し訳ありません。なお、データは離散型の比率尺度データです。

統計手法の選択ガイド
http://aoki2.si.gunma-u.ac.jp/FlowChart/Tutorial.html
では、「検定→その他」で適合度の検定のページ
http://aoki2.si.gunma-u.ac.jp/lecture/tekigoudo.html
にたどり着きますが、対応のある2群の散布度(分布の形)に違いがあるかを検証する分析が
見つけられず、質問させていただきました。

また、
http://aoki2.si.gunma-u.ac.jp/FlowChart/kentei_heikinti_taiouari.2.html
のページで紹介されている検定は平均値(代表値)の比較が目的であり、
散布度(分布の形)の違いを検証するものではないと思っており、判断に迷っております。


R 複数のNA判定に基づく代入
  投稿者:明石 2018/06/16(Sat) 19:06 No. 22571
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
毎々、ご丁寧なご教示をいただき、誠にありがとうございます。
改めて、御礼を申し上げます。

青木先生にご教示いただきたいことが出てきました。
何卒どうぞよろしくお願いいたします。

添付ファイルでお示しをします。 ご覧ください。

表1から表2を作成したいと思います。

表1のデータフレームは、以下のようになっています。
・socreは、すべてNA
・score1とscore2の組み合わせは、以下の3ケースのいずれか
  ・score1がNA,score2は数値
  ・score1が数値,score2はNA
  ・score1,score2両方ともにNA

score1、score2で、数値がセットされている側の値をscoreにセットする。

データフレームを dat として、Rコードを書きました。

jun <- which(! is.na(dat$score1))
dat$score[jun] <- dat$score1[jun]

jun <- which(! is.na(dat$score2))
dat$score[jun] <- dat$score2[jun]

score1、score2を、各々、NAの判定をして、scoreに代入しています。

これらを1つにまとめる書き方があれば、ご教示をお願いいたします。

お手数をおかけいたします。


Re: R 複数のNA判定に基づく代入
  投稿者:青木繁伸 2018/06/16(Sat) 21:35 No. 22572
いろいろやり方はあるでしょうが,もっとも真っ正直に書くならば
> dat = data.frame(
+ score=NA,
+ score1=c(NA, NA, 2.174, 6.042, NA, 2.656, NA, NA, NA, NA, 3.153, 4.344),
+ score2=c(1.481, 4.113, NA, NA, 2.674, NA, 2.500, 2.276, 4.162, NA, NA, NA))

> dat$score = ifelse(is.na(dat$score1) & !is.na(dat$score2), dat$score2, ifelse(!is.na(dat$score1) & is.na(dat$score2), dat$score1, NA))
> dat
score score1 score2
1 1.481 NA 1.481
2 4.113 NA 4.113
3 2.174 2.174 NA
4 6.042 6.042 NA
5 2.674 NA 2.674
6 2.656 2.656 NA
7 2.500 NA 2.500
8 2.276 NA 2.276
9 4.162 NA 4.162
10 NA NA NA
11 3.153 3.153 NA
12 4.344 4.344 NA
トリッキーな書き方ならば,出現する可能性のある数値よりも小さい値を dummy にわりあてて,以下のようにすることも可能
> dummy = -99999 # 出現する可能性のある数値が正ならばこのような値を仮定する
> dat[is.na(dat)] = dummy
> dat$score = ifelse(dat$score1 > dat$score2, dat$score1, dat$score2)
> dat[dat == dummy] = NA
> dat
score score1 score2
1 1.481 NA 1.481
2 4.113 NA 4.113
3 2.174 2.174 NA
4 6.042 6.042 NA
5 2.674 NA 2.674
6 2.656 2.656 NA
7 2.500 NA 2.500
8 2.276 NA 2.276
9 4.162 NA 4.162
10 NA NA NA
11 3.153 3.153 NA
12 4.344 4.344 NA
どちらが効率がよいか?

Re: R 複数のNA判定に基づく代入
  投稿者:明石 2018/06/16(Sat) 21:54 No. 22573
青木先生 様;

お忙しいところを失礼いたします、明石と申します。

今回も、ご丁寧なご教示をいただき、大変に良い勉強をさせていただきました。

2つの方法とも理解できましたが、ifelse()関数を使わせていただきます。

今回も助けていただき、誠にありがとうございました。


サンプリング誤差について
  投稿者:統計迷い人 2018/06/07(Thu) 11:20 No. 22563
サンプリング誤差について質問です。

サンプリング誤差について、視聴率での以下のサイトでの考え方や計算方法は分かったのですが、回答選択肢が3以上ある場合、例えば、非常に好き、やや好き、どちらでもない、やや嫌い、非常に嫌い、などでも、誤差について同様に考えて良いのでしょうか?

https://toukeigaku-jouhou.info/2015/08/25/audience-rate/


よろしくご教示のほど、お願いいたします。

Re: サンプリング誤差について
  投稿者:鈴木康弘 2018/06/10(Sun) 09:53 No. 22564
 サンプリング誤差って何かな?と思ったら、母集団の比率の95%信頼区間の求め方のことですか。
 選択肢がたくさんあっても同じだと、僕は思います..?

Re: サンプリング誤差について
  投稿者:統計迷い人 2018/06/10(Sun) 11:55 No. 22565 HomePage
鈴木康弘 先生

ご回答、ありがとうございます。

視聴率や知名度、購入経験のように、yes no いずれかから一つ選ぶ場合は二項分布であって、それが正規分布に近似することから、95%の信頼区間推定が出来るとテキストにあったのですが、

非常に好き、やや好き、普通、やや嫌い、非常に嫌いとかの、
多項分布の場合、誤差をどう考えたらよいかご教示下さい。

例えば、サンプル数600、非常に好き10%、やや好き15%、ふつう30%
やや嫌い20%、非常に嫌い25%の場合、
95%信頼区間で

非常に好き、非常に好き+やや好きは

どうなるのでしょうか?

よろしくお願いいたします。

Re: サンプリング誤差について
  投稿者:青木繁伸 2018/06/10(Sun) 15:45 No. 22566
各カテゴリーは排他なので,
「非常に好き」と「それ以外」で,sqrt(0.1*0.9/600)
「非常に好き+やや好き」と「それ以外」で,sqrt(0.25*0.75/600)
というように考えるということです。

Re: サンプリング誤差について
  投稿者:統計迷い人 2018/06/12(Tue) 10:16 No. 22570 HomePage
青木教授、アドバイス、大変、ありがとうございました。


信頼性係数 KR-20
  投稿者:やまと 2018/06/10(Sun) 18:02 No. 22567
よろしくお願いします。

RでKuder RichardsonのKR20を計算するためのスクリプトは公開されているのでしょうか。教えていただけますでしょうか。

Re: 信頼性係数 KR-20
  投稿者:青木繁伸 2018/06/10(Sun) 21:13 No. 22568
検索しましたか?

https://rdrr.io/github/cddesja/validateR/man/kr20.html
などもあります。

ここに書かれていることが難しければ,ともかく以下のようにすればよいでしょう。

devtools がまだインストールされていないなら,以下の 1 行を一度だけ
install.packages("devtools")

インストルされていた(インストールした)なら,以下の 2 行を一度だけ
library(devtools)
install_github("cddesja/validateR")
そして,kr20 を使うときには R を立ち上げるたびに,以下の 1 行を 一度だけ
library(validateR)
さて,準備が整ったら,データを以下のように準備する
例では 7 人の被検者に 6 項目の変数を 0/1 で準備
x = matrix(c(
1,1,0,1,1,0,
1,0,0,1,1,1,
0,1,1,0,1,1,
0,0,0,1,1,0,
0,1,0,1,0,1,
0,0,0,1,1,0,
0,0,1,0,0,0), byrow=TRUE, ncol=6)
その後,
kr20(x)
とすれば,結果は
>  kr20(x)
[1] 0.1382488
となる。

====

なお,やっていることは簡単なので,自分で以下のような関数を書けば,面倒なことは一切なしで結果が得られる。
計算方法につては Wikipedia にさえ書かれているが,例えば,
https://www.radford.edu/.../KR21%20&%20KR20%20examples.xls
には,Excel での計算例が書かれているので,計算方法を理解すれば以下のような関数を書くのは容易である(R がそこそこできればではあるが)
KR20 = function(x) {
cm = colMeans(x)
k = ncol(x)
k/(k-1)*(1-sum(cm*(1-cm))/var(rowSums(x)))
}

> KR20(x)
[1] 0.1382488

Re: 信頼性係数 KR-20
  投稿者:やまと 2018/06/11(Mon) 07:38 No. 22569
青木先生

大変詳細な説明をいただきありがとうございます。

無事解決できました。


t検定で処理していいでしょうか
  投稿者:コロン 2018/05/27(Sun) 07:54 No. 22556
いつもお世話になっております。

2クラス(実験群と統制群)で、授業を行う前と行った後の平均値の差を見たいのですが、分散分析(被験者間×被験者内)で分析をしたところ、交互作用は非有意、主効果がそれぞれの要因で有意、となりました。

解釈に苦しむところがあり、t検定でpreの平均を調べたところ、差がない、という結果となりました(等分散性の検定もOK)。

この結果を利用して、事前の能力が同じだと仮定して、postのみをt検定で処理してもよろしいでしょうか?実際に行ったところ、有意でした。

アドバイスをお願いいたします。

Re: t検定で処理していいでしょうか
  投稿者:鈴木康弘 2018/05/30(Wed) 11:29 No. 22557
 単に授業を行った後の2クラスの比較をしたいならt検定でいいんでしょうが、
2要因分散分析ではいけないんですか?

Re: t検定で処理していいでしょうか
  投稿者:コロン 2018/05/30(Wed) 19:18 No. 22558
鈴木康弘先生

お返事、ありがとうございます。簡単に説明させていただきますと以下のようになります。

・(要因A)被験者間要因(実験群と統制群)×(要因B)被験者内要因(preとpost)の2要因分散分析(ともに2水準)
・preにはなく、postにおける平均値に差があることが実験仮説
・分散分析の結果、交互作用がなく、要因AおよびBそれぞれに主効果があった(グラフに書きますと、ともに右肩上がりの平行線)

これを、実験仮説に即してうまく解釈できず...つまり、要因Aの主効果に差があると言うことは、「実験群の方が統制群よりも平均が高い」ということですよね?これだと「postで差が出ている」と解釈できないのではないかと考えた次第です。なので、t検定で行ったところ、うまい具合に出てくれたので、上のような質問をさせていただいた次第です。

もしよろしければ再度コメントをいただけますと幸いに存じます。

Re: t検定で処理していいでしょうか
  投稿者:鈴木康弘 2018/06/03(Sun) 07:05 No. 22562
 授業前に有意差はありませんでした、授業後には有意差はありました、と言うことはできるでしょう。

 でも、分散分析での交互作用はどうなのか?と問われたら、ありませんでした、と白状するしかないのでは。


3 群以上の比率の差を対比較で検定する方法
  投稿者:Bob 2018/05/31(Thu) 21:56 No. 22559
青木先生
お世話になっております。

統計の初歩を学習している者でございます。
以下、お教え頂ければ幸甚に存じます。

Q1:二値尺度において、3群以上を比較する場合、その手法は、2群の比較と同様、フィッシャー法でよいと理解していますが、その理解であっていますでしょうか?

Q2:この検定で有意さがついたと前提です。その後個々の検定をするには何法を使用したらよろしいのでしょうか?
数値型でしたら、ダン法やチューキー法かと存じますが比率の差については存じません。

大変初歩的なことで恐縮でございます。種々調べましたがよく分かりませんでした。
お教え頂ければ幸いでございます。

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:青木繁伸 2018/06/02(Sat) 18:42 No. 22560
フィッシャーの正確検定を使うなら,ボンフェローニ法を適用するかな?

カイ二乗検定を使うならば
K 群の比率の差の検定・多重比較
ライアンの方法
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Ryan.html
テューキーの方法
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Tukey.html

Re: 3 群以上の比率の差を対比較で検定する方法
  投稿者:Bob 2018/06/02(Sat) 19:59 No. 22561
青木先生

ご多忙の折、お教え頂き有難うございました。
とても助かりました。
今後もしっかり統計の勉強をしていきたいと思います。

取り急ぎお礼まで


予測区間の式
  投稿者:hikon 2018/05/16(Wed) 14:29 No. 22552
予測区間の式について教えて頂けましたら幸いです。

検索したところ、下記サイトがあり、式がすべて違い、値も違います。
どれを使えばよいのかわかりません。
コメント頂けましたら幸いです。
https://qiita.com/ksksk/items/75ba95337ccdb32e7cb1
https://ameblo.jp/cabay/entry-12184838325.html
http://heartland.geocities.jp/ecodata222/ed/edj1-1-4.html
https://www.monodukuri.com/gihou/article/946

それはそれとして、今回やりたいことは、
重回帰分析をした後に予測区間をプロットすることです。
お勧めのやり方をご教示頂けましたら助かります。
Rは少し使えます。

お手数おかけいたしますがよろしくお願いいたします。

Re: 予測区間の式
  投稿者:青木繁伸 2018/05/16(Wed) 16:58 No. 22553
2 番目のページと 3 番目のページは,標本平均値の信頼区間と予測区間の話です
1 番目のページと 4 番目のページは,単回帰分析における信頼区間と予測区間の話です

1 番目と 4 番目は,基本的に同じですが,

1 番目のページの式で tn-p-1((1_0.95)/2) の部分は,(1_0.95)/2 は,下に「自由度 n-p-1 の t 分布の (1+0.95)/2(1+0.95)/2 パーセンタイル」と書かれているので誤解を生じるのでしょう。1 番目のページの結果は正しいです。

4 番目のページの式は,t(n-2, α) を「α は有意水準で 5%,(1-α) が信頼水準で,ここへ 95% が入る。」と紛らわしいことを書いているが,信頼水準を 95% にするには t(n-2, α/2) でなければならない。つまり,上側確率が 0.025 である t 値を求める必要がある。αを使ってしまうと間違った答えが出る。

> 重回帰分析をした後に予測区間をプロットする

予測区間の求め方については,
http://sayalabo.com/labo4.html
に説明があります。

ただ,プロットは難しいでしょうね。

Re: 予測区間の式
  投稿者:hikon 2018/05/18(Fri) 08:39 No. 22554
早速のご回答ありがとうございます。
自分の勉強不足を痛感いたしました。
昨日一日勉強いたしました。
まだ理解できていませんが、取り急ぎお礼とご報告をさせて頂きます。
理解が進んだ後、また質問させて頂くこともあると思います。
今後ともよろしくお願いいたします。


ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データであ...
  投稿者:tosh 2018/05/09(Wed) 11:39 No. 22537
ロジスティック回帰分析を行う際、説明変数として用いたい変数が、3種類以上の値を含む名義尺度データであった場合に、それを複数のダミー変数に置き換える方法が、いろいろあることについて悩んでいます。

インターネットで調べられた限りですが、おもに次の(1)〜(3)の方法があるようでした。
これらはどのように使い分ければよいでしょうか。

-----------------

説明変数Aに、"x","y","z" の3種類の名義が含まれている場合、

(1)ダミー変数 { d1, d2, d3 } を作り、変数Aが、

"x" なら { 1, 0, 0 }
"y" なら { 0, 1, 0 }
"z" なら { 0, 0, 1 }

 に置き換える方法。

(2)ダミー変数 { d1, d2 } を作り、変数Aが、

"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { 0, 0 }

 に置き換える方法。

(3)ダミー変数 { d1, d2 } を作り、変数Aが、

"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { -1, -1 }

に置き換える方法。

-----------------

上に挙げたうち、どれか1つの方法を紹介しているサイトはいくつか見つかるのですが、それぞれの方法の比較や、どう使い分ければよいか、という点についての情報がなかなか見つからず、質問させていただきました。

また、この問題が体系的に解説されている書籍などで、もしおすすめのものがありましたら、ご紹介いただければ幸いです(それぞれの方法の結果の解釈、オッズ比への逆変換の方法、ステップワイズ法で変数を増減する場合にダミー変数はどう扱うのか、などについてもよく知りたいと思っています)。

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:後医は名医 2018/05/09(Wed) 12:52 No. 22538
少なくとも(1)は多重共線性でダメなのでは

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/10(Thu) 12:31 No. 22541
ご返信ありがとうございます。

多重共線性は、互いに相関が強い複数の変数が説明変数の中に含まれていると、起こる問題ですよね。
(1)の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…?


(追記)
すみません、次のように修正しました
修正後「 d1 が 1 なら d2 + d3 は必ず0に決定されてしまうので、」
修正前「 d1 が 1 なら d2 とd3 は必ず0に決定されてしまうので、」

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:青木繁伸 2018/05/10(Thu) 22:50 No. 22542
> 1)の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…?

そういうこと

R では,変数Aが名義尺度であるということを
FA = factor(A)
としておけば,
glm(result ~ FA)
でちゃんと処理してくれます。
変数 A は 2 個のダミー変数で表現できます。
3 個以上の変数を使っても,R では 3 個目からは分析に使われません。他の統計ソフトではエラーメッセージを吐いて死んでしまうこともあるでしょう。
独立な2変数変数であれば,どんな値であっても(0 と 1 でなくても -6.23123 と 578.21234 でもなんでも),その変数に対する係数は見かけは違いますが,その係数を使って計算される予測値は全く同じになります(そうでないと困ります。よい結果が出るような変数値のセットがある!なんてことになると,客観性がなくなります)。
ではなぜ,(3)のような変数を作らない!!かというと,出てくる値が3種類になるからです。0と1だと計算も簡単だしわかりやすい。
それにしても 0,1,-1 を与えるというのは,どんなメリットがあるのかなあ?

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/11(Fri) 16:38 No. 22543
ご返信、ありがとうございます。

(3)の意義については、この方法が載っていたpdf資料の内容を後でアップいたします。
(私のパソコンのDLフォルダの整理が悪く、見つかり次第)

通常は (2) の置き換えを行い、(1)はダメ、ということですね。

さっそく、R で、名義尺度の変数を factor型 にし、glm メソッド にかけてみました。
(スクリプトは、別に投稿します)

このスクリプトは次のようなロジスティック回帰分析を行うつもりで書きました。

・データセット: birthwt 低体重出生とそのリスク因子のデータ
・目的変数: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
・説明変数: 母親の年齢、母親の体重、母親の人種、...etc

母親の人種 race {1:白人 2:黒人 3:その他 } が名義尺度なので、スクリプトの中であらかじめ次のように変換しています。

# race を factor型 に変換する
DF1$race <- factor(DF1$race)

glmメソッド の summary をみると、以下に引用したように race2, race3 という2つの変数が生成されているのを確認しましたが、これらは race { 1:白人 2:黒人 3:その他 } をどのように変換した変数なのかを確認する方法がわかりません。

---------------------------------------------------------------
# 結果表示
summary(out1.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF1)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8946 -0.8212 -0.5316 0.9818 2.2125

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.48062 1.19689 0.402 0.68801
age -0.02955 0.03703 -0.798 0.42489
lwt -0.03397 0.01524 -2.229 0.02580 *
race2 1.27226 0.52736 2.413 0.01584 *
race3 0.88050 0.44078 1.998 0.04576 *
smoke 0.93885 0.40215 2.335 0.01957 *
ptl 0.54334 0.34540 1.573 0.11571
ht 1.86330 0.69753 2.671 0.00756 **
ui 0.76765 0.45932 1.671 0.09467 .
ftv 0.06530 0.17239 0.379 0.70484
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.28

Number of Fisher Scoring iterations: 4
---------------------------------------------------------------


ここでの {race2, race3} は、たとえばですが、内部的に

race=1 なら { 1, 0 }
race=2 なら { 0, 1 }
race=3 なら { 0, 0 }

のように変換されているのだと思いますが、それをどのように確認できるのでしょうか? 

str(out1.glm)

としてみたりもしましたが…

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/11(Fri) 16:40 No. 22544
以下、スクリプトです。
----------------------

###################################################
#
# ロジスティック回帰分析例
#
# https://www.gixo.jp/blog/3200/
# より一部改変
#

library(MASS)

str(birthwt)

# *** 目的変数 ***
# low:新生児の体重が2.5kg未満か否か(2.5kg未満が1)
#
# *** 説明変数 ***
# age:母親の年齢
# lwt:母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke:母親の喫煙有無(喫煙有りが1)
# ptl:母親の早産経験の有無(経験有りが1)
# ht:母親の高血圧の有無(有りが1)
# ui:母親の子宮神経過敏の有(有りが1)
# ftv:妊娠第1期における医師の訪問回数
#
# *** 使わない ***
# bwt:新生児の体重(グラム単位)

# ■ 前処理

# bwtは使わないので除く
DF1 <- birthwt[,1:9]

# lwt をポンド単位からグラム単位に変換する
DF1$lwt <- DF1$lwt * 0.454

# race を factor型 に変換する
DF1$race <- factor(DF1$race)

str(DF1)

# ■ ロジスティック回帰

out1.glm <- glm(
low ~ . ,
family = binomial(link = "logit"),
data = DF1)

# 注: モデル式 low ~ . は、目的変数に low、説明変数に low 以外の全ての変数を使うという意味

# 結果表示
summary(out1.glm)

# 各説明変数の係数を指数変換し、オッズ比を表示
exp(out1.glm$coefficients)

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:青木繁伸 2018/05/11(Fri) 21:30 No. 22545
例としてあげられたデータ birthwt$race は整数型で 1, 2, 3 をとる変数でしたが,factor(birthwt$race) で factor型の変数にすると,1, 2, 3 というレベルを持つ factor 変数になります。
> DF1$race
[1] 2 3 1 1 1 3 1 3 1 1 3 3 3 3 1 1 2 1 3 1 3 1 1 3 3 1 1 1 2 2 2 1
途中省略
[161] 1 3 3 3 2 1 3 3 1 1 2 2 2 3 3 1 1 1 1 2 3 3 1 3 1 3 3 2 1
Levels: 1 2 3 ← ここ

このような factor 変数を多変量解析などに使用すると,一番最初の Level が基準値 0 になります。つまり,必要なダミー変数において,取る値の全てが0
一般的には,変数を factor 変数に変換すると,"辞書順の最初の値が基準値 0 になります。
詳しくは factor 関数のオンラインヘルプを参照してもらいたいですが,ちょっと見以下のような結果になりますが,訳ワカメに陥らないようにしっかり理解してください。

> sex = c("male", "female", "U.K.", "male", "male", "female")
> factor(sex, labels=c(1,2,0))
[1] 2 1 0 2 2 1
Levels: 1 2 0
> factor(sex, levels=c("male", "female", "U.K."))
[1] male female U.K. male male female
Levels: male female U.K.
> factor(sex, levels=c("male", "female", "U.K."), labels=c(1,2,0))
[1] 1 2 0 1 1 2
Levels: 1 2 0
さて,このようなことから,R が作っているダミー変数を自分で作るとすると以下のようになります
> DF2 = DF1
> DF2$race = NULL # factor である race を削除
> Race = birthwt$race # 元のデータは整数型 1, 2, 3
> DF2$Race2 = (Race == 2)+0 # Race が 2 のときのみ整数値 1,それ以外なら 0
> DF2$Race3 = (Race == 3)+0 # Race が 3 のときのみ整数値 1,それ以外なら 0
> head(DF2) # Race = 1 のときは DF2$Race2 = DF2$Race3 = 0
low age lwt smoke ptl ht ui ftv Race2 Race3
85 0 19 82.628 0 0 0 1 0 1 0
86 0 33 70.370 0 0 0 0 3 0 1
87 0 20 47.670 1 0 0 0 1 0 0
88 0 21 49.032 1 0 0 1 2 0 0
89 0 18 48.578 1 0 0 1 0 0 0
91 0 21 56.296 0 0 0 0 0 0 1
> out2.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF2)
> summary(out2.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF2)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8946 -0.8212 -0.5316 0.9818 2.2125

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.48062 1.19689 0.402 0.68801
age -0.02955 0.03703 -0.798 0.42489
lwt -0.03397 0.01524 -2.229 0.02580
smoke 0.93885 0.40215 2.335 0.01957
ptl 0.54334 0.34540 1.573 0.11571
ht 1.86330 0.69753 2.671 0.00756
ui 0.76765 0.45932 1.671 0.09467
ftv 0.06530 0.17239 0.379 0.70484
Race2 1.27226 0.52736 2.413 0.01584
Race3 0.88050 0.44078 1.998 0.04576

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.28

Number of Fisher Scoring iterations: 4
ということで,DF1 を使ったときと同じ結果になりました。

なお,reorder という関数がありますが,この存在を知ると,きっと悪夢を見ます(^_^;)

複雑怪奇な変換規則を理解するよりは,自分で「適切な」ダミー変数を作成する方が幸せかも知れません。

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/12(Sat) 04:06 No. 22547
factor関数をテストしてみると、自動的に定まる Levels の順序が独特ですね…。

> factor(c(100,1,3,1,-5))
[1] 100 1 3 1 -5
Levels: -5 1 3 100 → 数値は昇順なのはまあ分かりますが、

> factor(c("b","a","A","c","A","b","C"))
[1] b a A c A b C
Levels: a A b c C → 文字列は辞書順だが小文字が先…?

最初のLevelがどの値になるか、いまいち確信が持てないので、ご説明にあるように、名義尺度をfactor型に変換する際は、 Levels パラメータ と Labels パラメータをちゃんと与えて、最初のLevelになる値を明示すればいい、と考え、次のように進めました。

> DF3 = DF1
> DF3$race = NULL
> DF3$race = factor(birthwt$race,levels=c(1,2,3),labels=c("白人","黒人","その他"))
> out3.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF3)

> summary(out3.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF3)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4806 1.1969 0.40 0.6880
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race黒人 1.2723 0.5274 2.41 0.0158 *
raceその他 0.8805 0.4408 2.00 0.0458 *

この結果が、お示しいただいた out.glm2 の結果と一致するのを見て、やっと気付いたのですが、out1.glm に含まれるダミー変数 race2 , race3 の "2"や"3"は、2:黒人、3:その他 のラベルなんですね。
(No.22543を投稿した時点では、通し番号のようなものかな、と誤解していました)

ダミー変数 { race黒人, raceその他 } は、おそらく

race=1 について { 0, 0 }
race=2 について { 1, 0 }
race=3 について { 0, 1 }

という対応になっていて、

> exp(out3.glm$coefficients)
(Intercept) age lwt smoke ptl
1.62 0.97 0.97 2.56 1.72
ht ui ftv race黒人 raceその他
6.44 2.15 1.07 3.57 2.41

から、

・人種が「黒人」であることは、白人である場合と較べ、低体重児の出生リスクが 3.57倍になる
・人種が「その他」  〃  2.41倍になる

という解釈に行き着くのか、と。

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/14(Mon) 19:48 No. 22550
投稿No. 22537に挙げた(3)の方法
(ダミー変数を -1, 0, 1 などとおく方法)

についてですが、掲載されているのは、ある大学病院で開催された「EBM・統計セミナー」という研修会のPDF資料でした(現在はDLできなくなっていました)。
この資料には、次の2つの変換方法が紹介されており、

ア)基準とするカテゴリ(reference group)を仮定する場合
イ)カテゴリ全体の平均をreferenceとする場合

-------------------
ア)基準とするカテゴリ(reference group)を仮定する場合

例 名義尺度データが{ "薬剤なし", "薬剤P", "薬剤Q" } のような場合

ダミー変数 D1, D2 を

 "薬剤なし" = {0,0}
 "薬剤P"  = {1,0}
 "薬剤Q"  = {0,1}

とおく。ロジスティック回帰の結果、D1, D2 について、

 b1 : D1 に対応する回帰係数
 b2 : D2 に対応する回帰係数

が推定され、

 exp(b1) は"薬剤なし"に対する"薬剤P"のオッズ比、
 exp(b2) は"薬剤なし"に対する"薬剤Q"のオッズ比、

という解釈となる。

--------------------
イ)カテゴリ全体の平均をreferenceとする場合

例 名義尺度データが{ "地域A", "地域B", "地域C" }のような場合

ダミー変数 D1, D2 を

 "地域A" = {1,0}
 "地域B" = {0,1}
 "地域C" = {-1,-1}

とおく。
 
 b1: "地域A"の回帰係数
 b2: "地域B"の回帰係数
 b3: "地域C"の回帰係数

とし、モデル上で、

 b1 + b2 + b3 = 0

という制約を置いたうえで推定する。
計算上は、b1, b2 しか推定されないが、b3 = - b1 - b2 として計算可能である。
exp(b1)をとればオッズ比が得られる。

--------------------
引用ここまでです。

最後の一文にある、「exp(b1)をとれば」は、「exp(-b1-b2)をとれば」ではないのかな…? という疑問がありますが、原文ママです。

なお、この資料の末尾の参考文献には次が挙げられていました。

青木繁伸(2009)「Rによる統計解析」
大橋靖雄(2011)「医師のための臨床統計学基礎編」
丹後俊郎(1996)「路地すてぃく回帰分析SASを利用した統計解析の実際」
豊田秀樹(1992)「原因をさぐる統計学」
服部環・海保博之(1996)「Q&A心理データ解析」
森實敏夫(2004)「入門 医療統計学」

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/14(Mon) 20:12 No. 22551
No.22545でお示しいただいた、自分でダミー変数を作る方法を真似して、(3)の方法でダミー変数を作成する処理を追加してみました。

###################################################
#
# ロジスティック回帰分析例 2
# https://www.gixo.jp/blog/3200/
# より一部改変

library(MASS)
str(birthwt)
#
# low: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
# age: 母親の年齢
# lwt: 母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke: 母親の喫煙有無(喫煙有りが1)
# ptl: 母親の早産経験の有無(経験有りが1)
# ht: 母親の高血圧の有無(有りが1)
# ui: 母親の子宮神経過敏の有(有りが1)
# ftv: 妊娠第1期における医師の訪問回数
# bwt: 新生児の体重(グラム単位)

options(digits = 4)

### DF1 ダミー変数を、Rで自動生成する

# bwtは使わないので除く
DF1 <- birthwt[,1:9]

# lwtをポンド単位からg単位に変換する
DF1$lwt <- DF1$lwt * 0.454

# raceをfactor型に変換する
DF1$race <- factor(DF1$race,levels = c(1,2,3),labels = c("白人","黒人","その他"))
str(DF1)

# ロジスティック回帰
out1.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF1)
summary(out1.glm)
exp(out1.glm$coefficients)

#### DF4 ダミー変数を、白人{1,0} 黒人{0,1} その他{-1,-1} とする方法で自作する

DF4 <- DF1

# ダミー変数に変換する
race = birthwt$race
DF4$race白人 = ifelse(race==1, 1, ifelse(race==2, 0, -1))
DF4$race黒人 = ifelse(race==1, 0, ifelse(race==2, 1, -1))
DF4$race = NULL
str(DF4)

# ロジスティック回帰
out4.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF4)
summary(out4.glm)
exp(out4.glm$coefficients)

# "raceその他" の偏回帰係数、オッズ比
b3 = - coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"]
names(b3) = c("raceその他_coef")

b3_exp = exp(- coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"])
names(b3_exp) = c("raceその他_exp(coef)")

-----------------
スクリプトここまでです。
以下、結果を貼ります。
-----------------

> summary(out4.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF4)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1982 1.1509 1.04 0.2978
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race白人 -0.7176 0.2699 -2.66 0.0079 **
race黒人 0.5547 0.3232 1.72 0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.3

Number of Fisher Scoring iterations: 4

> exp(out4.glm$coefficients)
(Intercept) age lwt smoke ptl ht ui ftv race白人
3.3142 0.9709 0.9666 2.5570 1.7217 6.4450 2.1547 1.0675 0.4879
race黒人
1.7414

> b3
raceその他_coef
0.1629
> b3_exp
raceその他_exp(coef)
1.177

------------
結果は以上です。
資料の説明からすると、

「低体重児の出生確率は、全人種の平均出生確率を基準にすると、
 ・白人の場合、0.4879倍
 ・黒人の場合、1.7414倍
 ・その他の人種の場合、1.177倍」

というような解釈になるのかな、と思いましたが…。
これ以上の詳細は、資料を作った方に直接質問するのが筋かと思いますが、具体的にやってみたことの報告のつもりで投稿しました。

[1] [2] [3] [4] [5] [6]

- J o y f u l  N o t e -
Modified by i s s o