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

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

[トップに戻る] [利用上の注意] [ワード検索] [過去ログ] [統計学自習ノート] [管理用]
latest article: No. 22685, 2019/02/21(Thu) 10:08
おなまえ
タイトル
コメント
ＵＲＬ
添付ファイル
暗証キー (英数字で8文字以内)
自分の投稿記事を修正・削除する場合 --> 処理 記事No 暗証キー

独立変数の水準サイズについて
　　投稿者：畠山 2019/02/20(Wed) 17:28 No. 22683
はじめて質問させて頂きます。

ロジスティック回帰分析で各独立変数の水準サイズはどの程度あれば良いのでしょうか？

「作業者」という独立変数が3つの水準（スタッフA、スタッフB、スタッフC)を持ち、
それぞれの数がスタッフA:250、スタッフB:500、スタッフC:1,000であった場合、得られた結果は各水準のサイズの差が大きく影響しますでしょうか？

あるreviewerから、この数の差が大きいので結果が信用ならないとのコメントを頂きました。
全体のサンプルサイズは独立変数×10以上であり、十分かと思うのですが、水準ごとの数についてわからず質問させて頂きました。

お忙しいところ恐縮ですが、宜しくお願い申し上げます。
Re: 独立変数の水準サイズについて
　　投稿者：青木繁伸 2019/02/20(Wed) 19:52 No. 22684
査読者は，水準ごとのサンプル数が 2，4 倍とかなり異なることを言っているのでしょう

この違いが恣意的なものなら問題になるでしょうが，各水準での抽出割合を同じにして標本を抽出すると，結果として母集団の構成を反映することになります。このような標本抽出は，むしろ好ましいことでしょう。逆に母集団で 1:2:4 の割合で存在するものを，抽出割合を無視して 1:1:1 でサンプルを取って分析したらその方がよほど問題だと思いますが。

このようなことは事前に抽出数を決定しない変数には常に存在するものです。サンプルを無作為抽出した後，いろいろな変数ごとに集計すると水準ごとのサンプルサイズには相当な違いが出てくるでしょう。全ての変数で，ある要素を持っているものと持たないものの割合が 1:1 になることなんてないでしょう。変数によっては数倍から数十倍になるものもあるでしょう。

どうも，査読者は古典的な実験で，対照群と各処理群のサンプルサイズを全て同じにすべしという例数設定と同じだと誤解しているのではないか？

例えば，社会調査のような非実験系の統計で，平均値の差の検定のような単純な場合でも，データを無作為抽出した後にいろいろな属性で再分類して検定を行うような場合，例数がかなり異なることはよくあることです。「そのような場合に検定は行えないし，行ったとしても信用ならない」なんて，誰もいいませんよね。

平均値の差の検定だろうがロジスティック回帰分析だろうが，「水準ごとのサンプルサイズが違うからだめだ」ということにならないのは明らかでしょう。
Re: 独立変数の水準サイズについて
　　投稿者：畠山 2019/02/21(Thu) 10:08 No. 22685
青木先生

ご返答ありがとうございます。
今回は臨床データを後方視的に検討したものだったので、恣意的なものではありませんでした。

各水準のサンプル数を意図的に変えたわけではない旨を上手く査読者の先生に伝えたいと思います。
お忙しいところご返答下さり誠にありがとうございました。


Rのreadlineによる入力待ちについて
　　投稿者：桜とお城 2019/01/18(Fri) 13:02 No. 22677
　いつも参考にさせて頂いております。
　　
　ある処理の結果を見て，変数に所定の数値をキーボードから入力してから，次の処理をするようにしたいので，以下の例をやってみました。

x=readline("input>")
x=as.numeric(x)
x*2

一行ずつ実行するならうまく行くけど，スクリプトで3行一括で実行すると，入力待ちになりません。

　基本的なことかもしれませんが，解決する方法がありませんか？どうぞよろしくお願い致します。
Re: Rのreadlineによる入力待ちについて
　　投稿者：青木繁伸 2019/01/25(Fri) 19:32 No. 22681
長い記事の投稿があって，気づきませんでした。

解決法は以下のようにするとよいかと。

つまり，前後を {　} で囲んで，複文にする。

{
x=readline("input>")
x=as.numeric(x)
x*2
}

ただ，これをそのままコンソールにコピーペーストしては何の改善もありません。

R のエディタにこれを書いて，この部分を選択して「実行」（Mac なら command+return)

> {
+ x=readline("input>")
+ x=as.numeric(x)
+ x*2
+ }
input>4
[1] 8

RStudio でも同様
Re: Rのreadlineによる入力待ちについて
　　投稿者：桜とお城 2019/01/30(Wed) 10:58 No. 22682
青木先生，ありがとうございました。問題なくできました。


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
青木先生

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

取り急ぎお礼まで
Re: 3 群以上の比率の差を対比較で検定する方法
　　投稿者：渡辺 2019/01/18(Fri) 18:10 No. 22678
青木先生

古い投稿への相乗りという形で失礼いたします。

「K群 x 2カテゴリー」でなく、

「K x N」のときのやり方を教えていただけますと幸いです。
Re: 3 群以上の比率の差を対比較で検定する方法
　　投稿者：青木繁伸 2019/01/21(Mon) 21:37 No. 22679
そのようなものは見たことがないですが，「K の全ての分割法」×「N の全ての分割法」で検定して，その全ての可能性を n として，個々の検定を α/n でやるということになるのではないかと思います。まあ，ちょっと，現実離れしているかなあとは思います。
Re: 3 群以上の比率の差を対比較で検定する方法
　　投稿者：渡辺 2019/01/23(Wed) 11:22 No. 22680
遺伝学解析で比較したい遺伝子型が複数あり、表現型のクライテリアも4カテゴリーあるデータを得ています。
やはり個別にやって補正するしかないのですね。
2x4のカイ2乗検定を遺伝子型の各組み合わせでやって、組み合わせ数で割ることにします。

お答えいただきありがとうございます。


GLMによるモデルと信頼区間の図
　　投稿者：ぐれいぷ 2019/01/17(Thu) 11:31 No. 22674
　一般化線形モデルによる解析による予測値と、その信頼区間の求め方についてお教え願います。

　7本の樹（A〜G）について、ある測定値（measure）と比率（30枚の葉のうち症状の出た葉（A）の割合）に関して、一般化線形モデル（ロジスティック回帰）で解析しています。元のデータとモデルによる曲線、さらに95％信頼区間を図示したいと考えています。実際のデータによるＲでの解析は以下の通りです（信頼区間は下側のみ示します）。

　図示した場合、信頼区間の破線がx軸の値の細かさで変動し、4から5で0.001刻みにすると階段状になってしまうのですが、モデルによる曲線と信頼区間を示す適切な方法をご存知でしたらお教えください。よろしくお願いします。

> a
tree measure A N
1 A 4.50 11 30
2 B 4.33 11 30
3 C 4.33 13 30
4 D 4.50 17 30
5 E 4.67 23 30
6 F 4.33 23 30
7 G 5.00 23 30

>
> glm <- glm(cbind(a$A,a$N-a$A)~a$measure,family=binomial)

> glm$coefficients
(Intercept) a$measure
-8.270453 1.900956

> ##############ｘ軸の値を0.1ごと
> x <- seq(4,5,by=0.1)
> pred.y <- 1/(1+exp(-(-8.270453+1.900956*x)))
> plot(a$measure,a$A/a$N,xlim=c(4,5),ylim=c(0,1))
> lines(x,pred.y)
>
> ##############信頼区間を求める
> mean <- length(a$N)/sum(1/a$N)
> mean
[1] 30
> lower <- qbinom(0.025,mean,pred.y)/mean
> lines(x,lower,lty=2)
>
>
> ##############ｘ軸の値を0.001ごと
> x <- seq(4,5,by=0.001)
> pred.y <- 1/(1+exp(-(-8.270453+1.900956*x)))
> plot(a$measure,a$A/a$N,xlim=c(4,5),ylim=c(0,1))
> lines(x,pred.y)
>
> lower <- qbinom(0.025,mean,pred.y)/mean
> lines(x,lower,lty=2)
Re: GLMによるモデルと信頼区間の図
　　投稿者：青木繁伸 2019/01/17(Thu) 14:05 No. 22675
qbinom(0.025,mean,pred.y)/mean は，qbinom(0.025,mean,pred.y) 自体が離散値なので，それを定数で割ったところで離散値には変わりない。したがって，「階段状に表示される」のが正確で正常な状態です。
Re: GLMによるモデルと信頼区間の図
　　投稿者：ぐれいぷ 2019/01/17(Thu) 14:39 No. 22676
青木先生
　お忙しい中早速のご回答ありがとうございました。モデルと信頼区間について勉強しなおします。今後ともよろしくお願いします。


オフセット項を含む単回帰モデルのパラメタ推定
　　投稿者：新春 2019/01/10(Thu) 22:28 No. 22670
オフセット項を含む単回帰モデルの切片の標準誤差の推定方法について，ご教示いただければ幸いです．
以下のようなＲのスクリプトを実行したとき，切片の推定値と標準誤差をsummary関数で確認できますが，標準誤差（以下の例では0.1114）をどのように算出しているのかわからず悩んでいます．適切な数式や考え方などが掲載されている資料があれば，ご紹介いただきたくお願い申し上げます．

> n <- 300
> x <- rnorm(n, mean=10, sd=1)
> e <- rnorm(n, mean=0, sd=2)
> a <- 7
> y <- a + x + e
> d <- data.frame(x,y)
> lm.out <- lm(y~1+offset(x), data=d)
> summary(lm.out)

Call:
lm(formula = y ~ 1 + offset(x), data = d)

Residuals:
Min 1Q Median 3Q Max
-6.2107 -1.2351 -0.0889 1.3389 5.1611

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.9713 0.1114 62.59 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.929 on 299 degrees of freedom

> (a.hat <- mean(y) - mean(x))
[1] 6.971308
> y.hat <- a.hat + x
> (sigma.e <- sqrt(sum((y - y.hat)^2)/(n-1)))
[1] 1.929243
Re: オフセット項を含む単回帰モデルのパラメタ推定
　　投稿者：新春 2019/01/13(Sun) 16:46 No. 22671
質問させていただいた件，誤差伝播の法則に基づく次の計算で，summary関数で表示される結果と一致する値が得られました．

> n <- 300
> x <- rnorm(n, mean=10, sd=1)
> e <- rnorm(n, mean=0, sd=2)
> a <- 7
> y <- a + x + e
> d <- data.frame(x,y)
> lm.out <- lm(y~1+offset(x), data=d)
> summary(lm.out)

Call:
lm(formula = y ~ 1 + offset(x), data = d)

Residuals:
Min 1Q Median 3Q Max
-5.4845 -1.3463 0.0459 1.3090 4.8227

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.0122 0.1123 62.42 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.946 on 299 degrees of freedom

> var.xbar <- sum((x - mean(x))^2)/(n*(n-1))
> var.ybar <- sum((y - mean(y))^2)/(n*(n-1))
> cov.xybar <- cov(x,y)/n
> (sigma.a <- sqrt(var.xbar + var.ybar - 2*cov.xybar))
[1] 0.1123454

以上，思いつきましたのでご報告いたします．
より適切な方法があれば，ご指摘いただければ幸いです．
Re: オフセット項を含む単回帰モデルのパラメタ推定
　　投稿者：青木繁伸 2019/01/13(Sun) 21:30 No. 22672
基本的には，lm や lm.fit など，関連するソースを読んでいけば良いと言うことになります。

ご存じとは思いますが，ソースを見るのは，コンソールに lm などと入力するだけです。

また，実際にどのように解析が進むのかをみるのは debug(lm) などとした後，lm を実行するということでよいでしょう。debug についての事前の知識は必要でしょうが。
Re: オフセット項を含む単回帰モデルのパラメタ推定
　　投稿者：新春 2019/01/14(Mon) 16:12 No. 22673
青木先生，

ソースコードの確認やdebugの活用など，ご指摘ありがとうございます．
丁寧にソースコードを追って，どのような計算が行われているのかを確認したいと思います．
ご助言いただけましたこと，また，このような質問の場を維持していただいていることに感謝申し上げます．


【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：明石 2019/01/02(Wed) 15:18 No. 22662
青木先生　様；

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

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

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

文字列の部分的マッチング charmatch
http://www.okadajp.org/RWiki/?R%E3%81%AE%E6%96%87%E5%AD%97%E5%88%97%E5%87%A6%E7%90%86%E9%96%A2%E6%95%B0#b71cbd6d

以下の例題が示されています。
複数の場合は「0」を返すということも確認しました。

> charmatch("me", c("mean", "median", "mode"))
[1] 0

> charmatch("med", c("mean", "median", "mode"))
[1] 2

この２つを組み合わせて、以下も確認できました。
> charmatch(c("me","med"), c("mean", "median", "mode"))
[1] 0 2

さて、以下の例題を考えます。
v <- c("aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
s1 <- c("aa","bb")
s2 <- c("ba","ab")

charmatch(s1, v)
charmatch(s2, v)

> charmatch(s1, v)
[1] 1 2
> charmatch(s2, v)
[1] 4 3

複数の場合は「0」を返すという理解ですので、なぜ、「0」にならないのか、分かりません。
ご教示いただけましたら助かります。

新年早々、お手数をおかけいたします。
何卒どうぞ、よろしくお願いいたします。
Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：青木繁伸 2019/01/03(Thu) 15:08 No. 22663
一番最初の定義があるからでは？

> 正確なマッチが一部だけのマッチより優先される(つまり、ターゲットの先頭部分に 値が正確にマッチし、更にターゲットのほうがより長いケース)。
> charmatch(c("aa","bb"), c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
+ )
[1] 0 3
> charmatch(c("aa","bb"), c("aa", "aas", "bb", "ab", "ba", "abab", "abba", "baba", "baab"))
[1] 1 3
> charmatch(c("aa","bb"), c("aaw", "aas", "bb", "ab", "ba", "abab", "abba", "baba", "baab"))
[1] 0 3

Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：明石 2019/01/03(Thu) 16:24 No. 22664
青木先生　様；

お忙しいところを失礼いたします、明石と申します。
新年明けましておめでとうございます。

青木先生がご指摘してくださいました箇所
> 正確なマッチが一部だけのマッチより優先される
> (つまり、ターゲットの先頭部分に 値が正確にマッチし、更にターゲットのほうがより長いケース)。

の意味を完全に理解できないでおりました。
大変に失礼いたしました。

---------------------------------------
青木先生、追加のご質問がございます。
何卒どうぞよろしくお願いいたします

> s <- c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
> v <- c("aa", "bb", "ab", "ba")
>
> grep(v[1], s)
[1] 1 2 9
> grep(v[2], s)
[1] 3 7
> grep(v[3], s)
[1] 4 6 7 8 9
> grep(v[4], s)
[1] 5 6 7 8 9

sを対象として、ベクトルvを、grep()関数で部分一致照合させた添字番号を取得したい
と思っていますが、
grep()関数では、第一引数は単体で与える必要があることから、
上記のように、vの要素ごとにgrep()関数を呼ぶことになります。

ベクトル v と s を与えて、上記の結果（4本のベクトル）をリストで得る
Ｒプログラムをご教示いただければ大変に勉強になります。

新年早々、お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。
Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：青木繁伸 2019/01/03(Thu) 18:32 No. 22665
sapply(v, grep, s) または lapply(v, grep, s) でできます
> s <- c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
> v <- c("aa", "bb", "ab", "ba")
> sapply(v, grep, s)
$aa
[1] 1 2 9

$bb
[1] 3 7

$ab
[1] 4 6 7 8 9

$ba
[1] 5 6 7 8 9

> lapply(v, grep, s)
[[1]]
[1] 1 2 9

[[2]]
[1] 3 7

[[3]]
[1] 4 6 7 8 9

[[4]]
[1] 5 6 7 8 9

Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：明石 2019/01/04(Fri) 08:19 No. 22666
青木先生　様；

お忙しいところを失礼いたします、明石と申します。
新年早々に、大変に良い勉強をさせていただきました。
この例題を通じて、今まで不得手であった sapply、lapplyも理解できました。
心から御礼を申し上げます。
Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：明石 2019/01/04(Fri) 17:07 No. 22667
青木先生　様；

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

ベクトル v と s を与えて、部分一致照合の添字を得る方法のご教示をいただき、誠にありがとうございました。
大変に良い勉強をさせていただきました。

s <- c("aa", "aa", "bb", "ab", "ba", "abab", "abba", "baba", "baab")
v <- c("aa", "bb", "ab", "ba")

lapply(v, grep, s)

[[1]]
[1] 1 2 9

[[2]]
[1] 3 7

[[3]]
[1] 4 6 7 8 9

[[4]]
[1] 5 6 7 8 9

取得した添字を用いて、要素を得ることができます。
> res <- lapply(v, grep, s)
> s[res[[1]]]
[1] "aa" "aa" "baab"
> s[res[[2]]]
[1] "bb" "abba"
> s[res[[3]]]
[1] "ab" "abab" "abba" "baba" "baab"
> s[res[[4]]]
[1] "ba" "abab" "abba" "baba" "baab"

---------------------------------------
青木先生、追加のご質問がございます。
何卒どうぞよろしくお願いいたします

ベクトル v と s を与えて、部分一致照合の要素
"aa" "aa" "baab"
"bb" "abba"
"ab" "abab" "abba" "baba" "baab"
"ba" "abab" "abba" "baba" "baab"
を得る方法をご教示いただければ、大変に勉強になります。

お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。
Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：青木繁伸 2019/01/04(Fri) 18:19 No. 22668
まあ，一行で書かなきゃいけないということでもないでしょうし，また，他に良い方法があるかもしれませんが，以下のようにすることでも出来なくはないです
> sapply(sapply(v, grep, s), function(i) s[i])
$aa
[1] "aa"   "aa"   "baab"

$bb
[1] "bb"   "abba"

$ab
[1] "ab"   "abab" "abba" "baba" "baab"

$ba
[1] "ba"   "abab" "abba" "baba" "baab"
Re: 【Ｒ】　文字列の部分的マッチング charmatch
　　投稿者：明石 2019/01/04(Fri) 18:58 No. 22669
青木先生　様；

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

sapply(v, grep, s)で得た添字を、s[]に代入する試みをしていましたが、結局できませんでした。

ご教示いただきました関数を定義する方法を理解できました。

ありがたいご教示に、心から、心より感謝いたします。
ありがとうございました。


【Ｒ】青木先生のAHPプログラムの使用方法
　　投稿者：明石 2018/12/20(Thu) 16:47 No. 22658
青木先生　様；

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

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

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

青木先生が作成されたAHP
http://aoki2.si.gunma-u.ac.jp/R/AHP.html

を使用するためのデータx y の与え方について、以下の例題でご確認をさせてください。

・評価基準（価格、品質、外見）
・代替案（商品Ａ、商品Ｂ、商品Ｃ、商品Ｄ）
・一対比較表
　添付の画像ファイルにお示しをします。

# 各評価基準の重要度

x <- c(1/5, 1/7, 1/3)

# 各代替案の重要度
y <- matrix(c(1/3, 1/5, 1/3, 1/3, 3, 5, 7, 7, 3, 1/3, 1/7, 1/5, 1/3, 3, 5, 3, 5, 3), 6, 3)

x と y の与え方は、よろしいでしょうか？

ご確認させてください。
お手数をおかけいたします。

Re: 【Ｒ】青木先生のAHPプログラムの使用方法
　　投稿者：青木繁伸 2018/12/21(Fri) 10:54 No. 22660
それでよいと思いますが
Re: 【Ｒ】青木先生のAHPプログラムの使用方法
　　投稿者：明石 2018/12/21(Fri) 18:10 No. 22661
青木先生　様；

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

いつもありがとうございます。

y <- matrix(c(1/3, 1/5, 1/3, 1/3, 3, 5, 7, 7, 3, 1/3, 1/7, 1/5, 1/3, 3, 5, 3, 5, 3), 6, 3)
で、行と列のサイズの与え方に不安がありましたので、ご確認をさせていただきました。

お手数をおかけいたしました、おかげさまで助かりました。

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

青木先生に、Ｒプログラムのご要望がございます。

意思決定分野のＡＨＰとＤＥＡ（包絡分析法）について、利用できるＲプログラムを調べております。

ＡＨＰについては、ＣＲＡＮのＲパッケージがありますが、データ作成方法が煩雑で、苦慮していました。

ＡＨＰについては、青木先生のプログラムのおかげで、とても助かりました

ＤＥＡについても、ＣＲＡＮのＲパッケージがありますが、マニュアルの理解で苦慮しています。

ＤＥＡについても、ぜひ、青木先生のＲプログラムの作成、公開をお願いしたいと思います。

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

いつもありがとうございます。


確証的因子分析のモデル適合について
　　投稿者：さとう 2018/12/17(Mon) 22:33 No. 22655
はじめて質問させていただきます。
現在、確証的因子分析を行っています。
重みづけのない最小二乗法で実施しましたらGFI、AGFI、NFIが0.9以上でした。
RMSEAの値が出ないため。なんの値を見ていければ良いでしょうか。
教えていただけないでしょうか？
サンプル数は293　正規性の検定では、正規性はみられませんでした。
よろしくお願いいたします。
Re: 確証的因子分析のモデル適合について
　　投稿者：波音 2018/12/20(Thu) 12:18 No. 22657
出力されているGFIやNFIを見て判断する、でいいと思いますが何か懸念されていることがあるのでしょうか、、

指標は挙げられているのも以外にもいくつかありますが、（個人の経験的に）どの指標を見ても何か結論が大きく変わることはないと思っています。

参考URL
http://daas.la.coocan.jp/sem/fitting_and_testing.htm
Re: 確証的因子分析のモデル適合について
　　投稿者：さとう 2018/12/20(Thu) 22:52 No. 22659
コメントありがとうございます。
現在、共分散構造分析を学んでいる最中でして、学習不足で申し訳ありません。
指標の持つ意味を理解することが大切であり、どの指標を見ても結論は大きく変わらないということが、わかりました。質問して学びになりました。
ありがとうございました。


時間依存性共変量のpropensity score
　　投稿者：徒弟 2018/12/18(Tue) 09:28 No. 22656
お世話になっております。

掲題の件で相談させて下さい。

ある患者コホートで、介入Xが、生存時間に影響するか評価したいと思ってます。
ただ、介入のタイミングが患者毎に異なっております。
すなわち介入は時間依存性共変量ということになります。

介入の有無を時間依存性共変量として、アウトカムを生存時間とするCox回帰で評価を行ってきましたが(結果は有意でした)、介入の有無で年齢等の患者背景が異なっており、propensity scoreマッチングのようなことをしたいと思ってます。

通常プロペンシティスコアは介入を従属変数としたロジスティック回帰で計算するのが一般的ですが、今回の場合介入が時間依存性共変量なので、Cox回帰でpropensity scoreを計算してはどうかと考えてます。

propensity scoreをCox回帰で計算するのは妥当でしょうか？
ご意見・ご指導賜れれば幸いです。


名義尺度の最小二乗フィッティングで追加制約がある場合の標準誤差の計算方法
　　投稿者：小西 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値からなるｎ×ｍ行列があります。
任意の列に対して、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
青木先生
回答有り難うございます。

複合アウトカムについては、死亡と再入院など、意味合いが異なるイベントを組み合わせることの意義について見直されてきているようです。
　現実的には、列挙したアウトカムの中でも１番重要視しているものをプライマリーアウトカムとし、それ以外はセカンダリーアウトカムとしてオマケぐらいの位置づけで解析をやるしかないのかもしれませんね。
　因みに、正準相関分析は私が使用している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群)の２群に分けました。100〜200例 VS　300〜500例　(集計を開始していないので見積もり概算です) 。
　アウトカムとして、連続変数と二値変数の２種類の変数を全部で10個くらい想定しています。交絡因子は年齢、性別、既往歴、など10種類くらいを検討しています。
　多変量解析になるかと思いますが、複合アウトカムの処理方法が分かりません。アウトカムを一つずつ、重回帰やロジスティックなどを多重性に気を付けながらやるしかないのでしょうか。

曖昧な質問ですみません。
お教え下さいませ。
Re: 複合アウトカム
　　投稿者：青木繁伸 2018/10/15(Mon) 19:21 No. 22620
結果が複数の変数ということならば，正準相関分析（重回帰分析の一般形）ということになるかな？実際に使った研究例は多くはないかも。


【Ｒ】　リスト構造からの読み取り
　　投稿者：明石 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)

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

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

所望する出力です。
Ａ,Ｂ
Ａ,Ｄ
Ｂ,Ｅ
Ｃ,Ｄ
Ｆ,Ｅ

--------------------------------
以下、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: 【Ｒ】　リスト構造からの読み取り
　　投稿者：青木繁伸 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: 【Ｒ】　リスト構造からの読み取り
　　投稿者：明石 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万件ほどの個人データがあり手作業での修正を覚悟しておりましたが、
先生のお力で無事補完することが出来ました。

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


【Ｒ】　二重繰り返し処理の高速化
　　投稿者：明石 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: 【Ｒ】　二重繰り返し処理の高速化
　　投稿者：青木繁伸 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: 【Ｒ】　二重繰り返し処理の高速化）
　　投稿者：明石 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には、注釈に、注１）年齢階級（20−29歳，30−39歳，40−49歳，50−59歳，60−69歳，70歳以上の６区分）と世帯員数（１人，２人，３人以上世帯の３区分）での調整値。割合に関する項目は直接法，平均 値に関する項目は共分散分析を用いて算出。 む

と書いてあります。

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

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

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

群間比較について，年齢調整を行い分析を行ったものについては，割合に関す る項目は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ですが，

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

「割合に関する項目は 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
青木先生

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

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

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

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


【Ｒ】　繰り返し処理の高速化
　　投稿者：明石 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: 【Ｒ】　繰り返し処理の高速化
　　投稿者：青木繁伸 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: 【Ｒ】　繰り返し処理の高速化）
　　投稿者：明石 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 検定の検定力は，条件を満たしているときの平均値の差の検定（ｔ 検定）の 95% ほどもあるので，劣っている検定手法を採用したと批難するのはちょっと厳しすぎるかも。

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

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

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

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

「データを見てから検定法を選択すること」はやってはいけないことと言えるでしょう。
Re: ノンパラ、パラメトリックの決め方
　　投稿者：コロン 2018/07/29(Sun) 08:08 No. 22583
青木先生

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

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

失礼いたします。

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

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