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

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


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

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


Rの一般化線形回帰のモデル式について
  投稿者:Y 2022/05/26(Thu) 19:29 No. 23187

↓この0が極端に多い計数値データについて、

0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0

これを、負の二項分布に近似させ、同分布のパラメータ(μ,k)のkをRを使って推定したいです(μは試料平均でよいとのことなので)。
Rのglm.nbで算出されるthetaをこのkにあてればよいと聞きました。

この関数については、モデル式が y~x1+x2+x3+1 のような目的変数と説明変数からなるモデルは豊富に紹介されているのですが、本件のように説明変数がない(?)使い方については事例が見つからずこまっています。
1例のみ glm.nb(y~1, とする事例を見つけたのですが、本件もこのモデル式でよいのでしょうか。
もし、よいとすれば、この y~1 という式はどのような意味を記述しているのでしょうか。
ご教示いただけると幸いです。

Re: Rの一般化線形回帰のモデル式について
  投稿者:aoki 2022/05/26(Thu) 23:39 No. 23188
いわゆる null model ですね。

独立変数がないので,予測値は,「どんな場合であっても全データの平均値」とするものです。

> glm(y ~ 1)

Call: glm(formula = y ~ 1)

Coefficients:
(Intercept)
0.07964601769912

Degrees of Freedom: 112 Total (i.e. Null); 112 Residual
Null Deviance: 16.28318584071
Residual Deviance: 16.28318584071 AIC: 105.7703175934
となりますが,Intercept = 0.07964601769912 ということは,どのような独立変数であっても(実際は独立変数はないのですが)予測値は 0.07964601769912 だ!という,どうでもよいモデルということです。

実際,
y = c(0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0)
の平均値は
> mean(y)
[1] 0.07964601769911507
ですから。モデルもへったくれもないです。

また,glm であろうが lm であろうが,(お飾りの部分はなくなりますが)主要な結果は同じになります。
> lm(y ~ 1)

Call:
lm(formula = y ~ 1)

Coefficients:
(Intercept)
0.07964601769912


このデータ,どういう条件下で得られたものですか?
「独立変数がない」なんて...

Re: Rの一般化線形回帰のモデル式について
  投稿者:Y 2022/05/27(Fri) 09:14 No. 23190
詳細なご回答、誠にありがとうございます。
根本的に間違ったことをしていたようです。

このデータは、ある製品の製造に使用する水の細菌検査の結果です。数値は一定量の水から検出された細菌の個数です。毎日1サンプルを検査し、時系列にならべたものです。

当業界では、この種のデータは、負の二項分布、ゼロ過剰ポアソン分布等に従うので、これらの中でよくフィットするものを選択し、推定された分布の例えば99.9パーセンタイル点を管理基準値にする・・・といったことが行われています。
添付文献の方法(SASによるパラメータ推定)をそのままやろうとしているのですが、SASがないのでRでの演算方法を調べていたのです。この部分以外は数式が記載されているので、エクセルで処理できました。
また、この文献の事例データは、今回の0ばかりのデータとは違うため、そのままではダメだろうと思ってはいたのですが・・・


Ljung-Box統計量について
  投稿者:統計の子猫 2022/03/05(Sat) 18:41 No. 23184
お世話になっております。本サイトをたびたび拝見し、勉強させていただいております。
私の考えが正しいかどうか、お伺いしたく質問をさせていただきます。

ある時系列の自己相関の集まりがゼロと異なるかどうかを調べる統計的検定の1つに、Ljung-Box(リュング・ボックス)統計量があります。
その式は、
  n(n+2) \sum_{j=1}^h\frac{\hat{\rho}^2_j}{n-j}
です。ここで、nは自己相関を推定するデータ数、\hat{\rho}は時間差 j における自己相関の値、hは検定する時間差です。
ある時系列の自己相関がゼロということは、その時系列は線形的には、独立で同一な分布を持つ乱数であるという扱いであると理解しております。

Ljung-Box統計量では自己相関を用いますが、自己相関ではなく相互情報量を用いた場合でも、同じ議論が出来るように思いますが、如何でしょうか。
もし、ある時系列が、独立で同一な分布を持つ乱数として扱えるのであれば、言い換えれば、無相関な時系列であれば、その相互情報量も理論上
(理想的には)ゼロになると思います。
そうであれば、Ljung-Box統計量の式の自己相関を相互情報量に入れ替えて、ある時系列の相互情報量の集まりがゼロと異なるかどうかを調べる
ことできる統計量として、そのままこの式を使っても良いように思います。

先生のお考えをお聞かせください。よろしくお願いいたします。

Re: Ljung-Box統計量について
  投稿者:aoki 2022/03/06(Sun) 21:24 No. 23185
残念ながら,私にはわかりかねます。
然るべき場所でお問い合わせください。

Re: Ljung-Box統計量について
  投稿者:統計の子猫 2022/03/06(Sun) 22:02 No. 23186
お忙しい中、返信、有り難うございます。
再度、自分で考えてみます。


小数以下を揃える
  投稿者:コロン 2021/12/12(Sun) 12:22 No. 23172
お世話になっております。

早速ではございますが,小数点を揃える方法をご教示いただけますでしょうか。

mean1 <- round(mean(d[method=="〇" & personal=="〇", 3]), digits=2)
mean2 <- round(mean(d[method=="△" & personal=="△", 3]), digits=2)
mean3 <- round(mean(d[method=="■" & personal=="■", 3]), digits=2)

mean1とmean3は小数点以下が3桁以上となるため小数第2まで表示させることはできているのですが,mean2が30(.00)となり,mean1やmean3と同じように小数第2まで表示されているかと思ったら,小数点以下は表示されていませんでした。mean2の30を30.00にするにはどのようにしたら良いのでしょうか。

formatでやってみたのですが,文字列として処理されてしまっているようで,その後の四則演算でエラーが起きてしまいました。

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

Re: 小数以下を揃える
  投稿者:aoki 2021/12/12(Sun) 14:07 No. 23173
コピペして動く形で提示してください。d なんか不要で邪魔なだけです。なんの情報ももっていません。

以下のようなことなんでしょう?

> round(31.2345, digits=2)
[1] 31.23
> round(30, digits=2)
[1] 30
> round(31.45678, digits=2)
[1] 31.46

> formatでやってみたのですが,文字列として処理されてしまっているようで,その後の四則演算でエラーが起きてしまいました。

についても,実際に書いたプログラムを示してください。

sprintf は文字列に変換するものですから,その後四則演算するとエラーになるのはあたりまえですが,そもそも,その後四則演算するのなら,小数点以下2桁で丸めたのでは誤差を含むことになるので,意味がないですね。そのまま計算して,最終結果の表示のときに文字列として表示するものでしょう。

> sprintf("%.2f", 31.2345)
[1] "31.23"
> sprintf("%.2f", 30)
[1] "30.00"
> sprintf("%.2f", 31.45678)
[1] "31.46"

Re: 小数以下を揃える
  投稿者:コロン 2021/12/12(Sun) 14:30 No. 23174
青木先生

申し訳ございませんでした。再度ご教示くださいませ。

使用データは次の通りで,これをdとしました。

method personal score
黒板 低 18
黒板 低 6
黒板 低 29
黒板 中 20
黒板 中 12
黒板 中 31
黒板 高 40
黒板 高 49
黒板 高 31

以下はplot等を使ったグラフとなります。このグラフに平均点を桁数を揃えて表示させたいのです。

x <- c(1, 2, 3)
mean4 <- round(mean(d[method=="黒板" & personal=="低", 3]), digits=2)
mean5 <- round(mean(d[method=="黒板" & personal=="中", 3]), digits=2)
mean6 <- round(mean(d[method=="黒板" & personal=="高", 3]), digits=2)

allmean2 <- c(mean4, mean5, mean6)

plot(x, allmean2, xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(0.5, 3.5), ylim=c(0, 50), type="b", pch=2, lwd=2, lty=2, col=3)

text(x=1, mean4-3, mean4)
text(x=2, mean5+3, mean5)
text(x=3, mean6+3, mean6)

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

Re: 小数以下を揃える
  投稿者:aoki 2021/12/12(Sun) 15:14 No. 23175
No. 23173 の最後に書いていることそのままですが,以下でよいでしょう。

text(x=1, mean4-3, sprintf("%.2f", mean4))
text(x=2, mean5+3, sprintf("%.2f", mean5))
text(x=3, mean6+3, sprintf("%.2f", mean6))

なお,最近 R を書くことがメッキリ減って,
d[method=="黒板" & personal=="低", 3]
で出るエラーの理由がわからなくて焦りました。

Re: 小数以下を揃える
  投稿者:コロン 2021/12/12(Sun) 20:54 No. 23176
青木先生

お休みのところ、ありがとうございます。text内で対応するわけですね。

ところで、d[method==.....のところなのですが、私のRではエラーが出ませんでしたが、何か問題がございますでしょうか。もしよろしければご指導いただけますでしょうか。

Re: 小数以下を揃える
  投稿者:aoki 2021/12/12(Sun) 22:45 No. 23177
エラーが出ない??

はてさて。

method personal score
黒板 低 18
黒板 低 6
黒板 低 29
黒板 中 20
黒板 中 12
黒板 中 31
黒板 高 40
黒板 高 49
黒板 高 31

というデータということで,わざわざ,普通のエディタを立ち上げ,そのようなデータファイルを作りました。Mac の terminal でやれば

foor [1] > cat test.dat

method personal score
黒板 低 18
黒板 低 6
黒板 低 29
黒板 中 20
黒板 中 12
黒板 中 31
黒板 高 40
黒板 高 49
黒板 高 31

ですね。

これを,reda.table で読み込んで,

> d <- read.table("test.dat", header=TRUE)

表示してみると,ちゃんと読み込めています

> d
1 method personal score
2 黒板 低 18
3 黒板 低 6
4 黒板 低 29
5 黒板 中 20
6 黒板 中 12
7 黒板 中 31
8 黒板 高 40
9 黒板 高 49
10 黒板 高 31

で,あなたが書いたようにやってみると

> d[method=="黒板" & personal=="低", 3]
Error in `[.data.frame`(d, method == "黒板" & personal == "低", 3) :
object 'method' not found

エラーになります。

しばし熟考後,以下のようなんだろうなと

> d[d$method=="黒板" & d$personal=="低", 3]
[1] 18 6 29

え〜〜と,tydyvers とか dplyr とか,ほかのナントカライブラリ,を使えばできるんでしょうか?時代遅れなので,pure な R しか使えないんです。

わざわざデータファイルを作らないといけないとか,余分なライブラリを読んでおかないといけないとか(少なくとも,私は dplyr とか tidyverse その他は,使わない),(追試してみようとする他のユーザの方のためにも)他の余分な前提条件が必要ないように提示するのがよろしいかと。(便利ならば,そのようなライブラリーを使うように奨励すればよいとは思いますが)l

つまるところは,例えば,

> d <- structure(list(method = c("黒板", "黒板", "黒板", "黒板",
"黒板", "黒板", "黒板", "黒板", "黒板"), personal = c("低",
"低", "低", "中", "中", "中", "高", "高", "高"), score = c(18L,
6L, 29L, 20L, 12L, 31L, 40L, 49L, 31L)), class = "data.frame", row.names = c(NA,
-9L))

> mean4 <- round(mean(d[d$method=="黒板" & d$personal=="低", 3]), digits=2)
> mean4
[1] 17.67

と提示すればよろしいのかなと。

逆に,
mean4 <- round(mean(d[method=="黒板" & personal=="低", 3]), digits=2)
で mean4 がなぜ求まるのかな〜〜と思ったわけです。(時代遅れ?)

Re: 小数以下を揃える
  投稿者:aoki 2021/12/12(Sun) 23:17 No. 23178
長くなりましたので,総括します。

潜在回答者に負担をかけないような,具体的な質問の仕方をお願いします。

それは,あなたのためでもあるのです。

Re: 小数以下を揃える
  投稿者:コロン 2021/12/13(Mon) 09:37 No. 23179
青木先生

お忙しい中,ご説明くださりありがとうございます。

はい,私の元々のコードではなぜかエラーは表示されませんでした。

先生が書いてくださったコードを見て,(関係はないのかも知れませんが)VBAの以下を思い出しました。

ws1.Range(ws1.cells(), ws1.cells())

です。

ありがとうございました。

Re: 小数以下を揃える
  投稿者:aoki 2021/12/13(Mon) 10:58 No. 23180
どこか前の方で,

> method <- c("黒板", "黒板", "黒板", "黒板", "黒板", "黒板", "黒板", "黒板", "黒板")
> personal <- c("低", "低", "低", "中", "中", "中", "高", "高", "高")

とかやってませんでしたか?
そのあとで d <- dataframe(method, personal, score) としたとか。

> d[method=="黒板" & personal=="低", 3]
[1] 18 6 29

[ ] の中の method も personal も データフレーム中の method,personal ではありません。

Re: 小数以下を揃える
  投稿者:aoki 2021/12/13(Mon) 11:06 No. 23181
どこか前の方で,

> method <- c("黒板", "黒板", "黒板", "黒板", "黒板", "黒板", "黒板", "黒板", "黒板")
> personal <- c("低", "低", "低", "中", "中", "中", "高", "高", "高")

とかやってませんでしたか?
そのあとで d <- dataframe(method, personal, score) としたとか。

> d[method=="黒板" & personal=="低", 3]
[1] 18 6 29

[ ] の中の method も personal も データフレーム中の method,personal ではありません。


A群とB群のANCOVA
  投稿者:当直中 2021/12/05(Sun) 20:58 No. 23170
A群(200人)とB群(800人)の血糖値の差を年齢、性、BMIなどなどを調整したうえで、ANCOVAで違いがあるかを検討していました。
A群とB群ではもともと色々背景に違いがあるのでこれらをマッチさせて比較するように言われました。そこでプロペンシティスコアでマッチさせた上で、A群とB群の血糖値を比較しようとおもっています(1対1、1対2など?)。
マッチングした後はA群とB群の対象者に「対応」ができるとおもうのですが、ANCOVAでこの2群の血糖値の差を検定する場合、この「対応」を考慮しないといけませんか?A群(200人)とB群(800人)の血糖値の差をみたときと同じプログラムでよいですか?
SASを用いて解析しています。

Re: A群とB群のANCOVA
  投稿者:aoki 2021/12/06(Mon) 22:15 No. 23171
マッチングは難しいです。
差をもたらす変数をマッチングに使ってしまうと,当たり前ですが,差がなくなります。


片側検定か両側検定か?
  投稿者:50代のおじさん 2021/11/05(Fri) 11:55 No. 23166
先生のHPの中にある演習問題で以下の問題があります。

問題3 新しく開発された製品の品質は従来品に比べて優れていることを示したい。a,b いずれの検定手法をとるべきか解答欄に記入し,送信ボタンをクリックしなさい。

答えはaの「片側検定」なのですが,本当でしょうか?片側検定の方が両側検定より有意になりやすいので,一般的に検定が厳しい両側検定を用いる事が主流と聞いた事があります。実際は,どうなのでしょうか?

Re: 片側検定か両側検定か?
  投稿者:50代のおじさん 2021/11/05(Fri) 11:58 No. 23167
参考先です。これを読んで,両側検定をした方が良いのかなと思いました。

https://www.dynacom.co.jp/product_service/packages/snpalyze/sa_t2_oneside-twoside.html

Re: 片側検定か両側検定か?
  投稿者:aoki 2021/11/06(Sat) 12:58 No. 23168
αエラーだけではなく,βエラーも考える必要があります。
http://aoki2.si.gunma-u.ac.jp/lecture/Kentei/beta-error.html

いずれにせよ,理論的に方向性があるならば片側検定を使うべきです。

「片側検定の方が両側検定より有意になりやすいので,一般的に検定が厳しい両側検定を用いる」なら,有意水準を下げてやればよいだけです。つまり,有意水準 0.025 の片側検定と,有意水準 0.05 の両側検定は同じなのですから。

Re: 片側検定か両側検定か?
  投稿者:50代のおじさん 2021/11/09(Tue) 14:56 No. 23169
ご返信ありがとうございました。また,ご回答が遅くなり,申し訳ありませんでした。αエラー,βエラーを知らないので,まずはそこから勉強したいと思います。もし,まだ分からない点がありましたら,再度,質問させて下さい。


比率の差の検定
  投稿者:50代のおじさん 2021/11/03(Wed) 15:05 No. 23158
もう独学で20年位,統計学を勉強していますが,未だに基本が良く分かっていない者です。初めて投稿します。

比率の差の検定の質問です。

以下の様なデータがあったとします。
    賛成    反対   合計
男   22     50    72
女   51     30    81
合計  73     80    153

比率 0.3013699 0.6250000

賛成の中の男性の割合は0.3013699,反対の中の男性の割合は0.6250000であった。この比率に有意差はあるだろうか?

Rを用いて以下のコマンドを打ちました。

> prop.test(c(22,50),c(73,80))

Rの結果を以下に示す。

2-sample test for equality of proportions with continuity
correction

data: c(22, 50) out of c(73, 80)
X-squared = 14.774, df = 1, p-value = 0.0001212
alternative hypothesis: two.sided
95 percent confidence interval:
-0.4861744 -0.1610859
sample estimates:
prop 1 prop 2
0.3013699 0.6250000

P値が0.05未満なので,この2つの比率には差があると言える。

次に片側検定を行う(イェーツの連続の補正を外していますが気にしないで下さい)。

> prop.test(c(22,50),c(73,80),correct=FALSE,alternative="less")

Rの結果を以下に示す。

2-sample test for equality of proportions without continuity
correction

data: c(22, 50) out of c(73, 80)
X-squared = 16.047, df = 1, p-value = 3.09e-05
alternative hypothesis: less
95 percent confidence interval:
-1.000000 -0.198212
sample estimates:
prop 1 prop 2
0.3013699 0.6250000

P値は0.0000309と0.05未満なので,前者の比率(0.3013699)は後者の比率(0.6250000)よりも小さいと言える。

この様な片側検定は,χ二乗検定では出来ないと考えて良いでしょうか?Rでchisq.test関数を持ちてトライしたのですが,片側検定が出来るオプションは無い感じでした。

以上,よろしくお願いします。

Re: 比率の差の検定
  投稿者:aoki 2021/11/03(Wed) 20:22 No. 23160
結論から先にいうと,prop.test は面倒見がよいのに対して,chisq.test はソンナノシッタコッチャネーと素っ気ないということです。

prop.test() のソースリストでは,p 値(PVAL)を決めるのは以下の箇所です(注釈は私が付け加えました)。

if (alternative == "two.sided") # 両側検定の場合
PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
else { # 片側検定の場合
if (k == 1) # 母比率の検定の場合
z <- sign(ESTIMATE - p) * sqrt(STATISTIC)
else z <- sign(DELTA) * sqrt(STATISTIC) # 二群の比較の場合
# 検定統計量の平方根をとり
# 符号は DELTA > 0 のとき正,DELTA < 0 のとき負
PVAL <- pnorm(z, lower.tail = (alternative == "less"))
# p 値の設定
# alternative == "less" のとき,下側確率
}

# 例示されたデータについて分析してみます

tbl <- matrix(c(22, 50, 51, 30), byrow=TRUE, ncol=2)
print(tbl)

[,1] [,2]
[1,] 22 50
[2,] 51 30

# まずは prop.test() の結果です

result.prop.test <- prop.test(c(22, 50), c(73, 80), correct=FALSE, alternative="less")
print(result.prop.test)

2-sample test for equality of proportions without continuity
correction

data: c(22, 50) out of c(73, 80)
X-squared = 16.047, df = 1, p-value = 3.09e-05
alternative hypothesis: less
95 percent confidence interval:
-1.000000 -0.198212
sample estimates:
prop 1 prop 2
0.3013699 0.6250000

# プログラム中での各変数の値の確認をします

DELTA <- unname(result.prop.test$estimate[1] - result.prop.test$estimate[2])
print(DELTA) # 比率の差(符号付き)

[1] -0.3236301

STATISTIC <- unname(result.prop.test$statistic)
print(STATISTIC) # 検定統計量(χ2分布に従う)

[1] 16.04666

print(sqrt(STATISTIC)) # χ2分布統計量の平方根は標準正規分布統計量

[1] 4.005828

print(sign(DELTA)) # 符号の決定

[1] -1

z <- sign(DELTA) * sqrt(STATISTIC)
print(z) # 符号付きの標準正規分布統計量

[1] -4.005828

# ★★ 後で参照するためのブックマーク

alternative <- "less"
PVAL <- pnorm(z, lower.tail = (alternative == "less"))
print(PVAL) # 片側 p 値

[1] 3.090026e-05

prop.test からの結果。上の PVAL と一致していることを確認。

print(result.prop.test$p.value)

[1] 3.090026e-05

さて,ここからが回答です。
prop.test() に対して,chisq.test() は,prop.test() と違い,実に素っ気ない。

prop.test() がやってくれたことは,chisq.test() では自分でするしかない。

result.chisq.test = chisq.test(tbl, correct=FALSE)
print(result.chisq.test) # chisq.test の結果(片側検定なんて知らないから!と,素っ気ない)

Pearson's Chi-squared test

data: tbl
X-squared = 16.047, df = 1, p-value = 6.18e-05

片側検定の場合,自分でやるしかない(prop.test() はやってくれたのに)。

z <- sqrt(result.chisq.test$statistic) # 標準正規分布に従う統計量を自分で計算する
print(z)

X-squared
4.005828

p.value <- pnorm(-z) # 「適切に」 z の符号を付けること
print(p.value)

X-squared
3.090026e-05

上のブックマーク ★★ と同じ結果になったことを確認する。

Re: 比率の差の検定
  投稿者:50代のおじさん 2021/11/04(Thu) 09:33 No. 23161
大変,丁寧なご回答,ありがとうございます。

私は,Rのプログラミングを習得していないので,まずは,先生の回答を理解する事から始めたいと思います。ですので,少々,お時間を下さい。

Re: 比率の差の検定
  投稿者:50代のおじさん 2021/11/04(Thu) 14:17 No. 23162
プログラミングの箇所で質問です。

alternative <- "less"
PVAL <- pnorm(z, lower.tail = (alternative == "less"))

となっていますが,二行目は単に
PVAL <- pnorm(z, lower.tail = TRUE)

と同じ事でしょうか(計算させると同じ数値「3.090026e-05」となります)?その場合,何故,alternative <- "less"と定義して,PVAL <- pnorm(z, lower.tail = (alternative == "less"))とややこしく表記されるのかが分かりませんでした。もしかしたら,下側確率を求める事を強調されているのでしょうか?

Re: 比率の差の検定
  投稿者:aoki 2021/11/04(Thu) 22:01 No. 23163
> alternative <- "less"
> PVAL <- pnorm(z, lower.tail = (alternative == "less"))

と書いたのは,関数の引数で,「alternative <- "less"」と指定したのを引用しただけですので,引数で何が指定されたかを暗喩すれば,

> PVAL <- pnorm(z, lower.tail = (alternative == "less"))

だけでよいのです。

Re: 比率の差の検定
  投稿者:50代のおじさん 2021/11/05(Fri) 08:34 No. 23164
お答えありがとうございます。

Rで,

alternative == "less"

と入力すると,

[1] TRUE

とアウトプットされたので,納得しました。

引き続き,プログラムを読み進みます。また,質問させて下さい。

Re: 比率の差の検定
  投稿者:50代のおじさん 2021/11/05(Fri) 09:48 No. 23165
Rのプログラム,理解出来ました。プログラムの中にそれぞれの命令に対して注釈を付けて頂き,それで理解する事が出来ました。ありがとうございました。


質問項目文末が異なる比較
  投稿者:コロン 2021/11/01(Mon) 15:20 No. 23154
お世話になっております。

論文を読んでいて気になることが出てきたため,ご意見を頂ければ幸いに存じます。

事前事後のアンケートなのですが,事前では「〜できそうですか?」と問い,事後では「〜ができましたか」(* 〜 は同じ内容)と問い,4件法で間隔尺度と見なして平均値を出し,対応のあるt検定をしています。みている内容が異なるのにこのような検定は可能なのでしょうか?

もし不可能であればですが,このような場合,なにかよい方法はございますか。

Re: 質問項目文末が異なる比較
  投稿者:aoki 2021/11/01(Mon) 19:41 No. 23156
確かに,予測値(希望値)と実測値という点では異なるということになるかも知れませんが,このような違いは多かれ少なかれあるのではないでしょうか?
たとえば,投薬前後での症状(検査値)も,投薬有り・無しという点では異なりますが,それは,薬効を見るのであるから問題ないということになるのでは?
或いは,投薬前に主治医がこの患者の薬効はこの程度だろうと予測して,実際の投薬後の薬効(実測値)と比較するということになると,コロンさんの提示した状況と同じようになるのではないでしょうか?

Re: 質問項目文末が異なる比較
  投稿者:コロン 2021/11/02(Tue) 07:29 No. 23157
青木先生

お返事いただきありがとうございます。ご提示いただきました2つの例で納得いたしました。

ありがとうございました。


2016/02/17 〜 2021/09/13 の記事を書庫へ移動
  投稿者:aoki 2021/10/08(Fri) 14:05 No. 23151
書庫 048 2016/02/17 --- 2021/09/13
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc048/index.html


信頼区間
  投稿者:ビギナー 2021/09/11(Sat) 20:15 No. 23140
信頼区間について、,鉢△里海箸わかりませんので、ご教授いただけないでしょうか?
よろしくお願い申し上げます。

.イドラインに「パラメータが正規分布に従わない場合には、ノンパラメトリック法で求めた90%信頼区間を判定に用いてもよい。」と記載されていたのですが、信頼区間はノンパラメトリック法でも求めることができるのでしょうか?

▲イドラインに対称信頼区間、最短非対称信頼区間の用語が出てくるのですが、これらの用語の説明がありませんので、意味がわかりませんでした。ご教授いただけないでしょうか?

Re: 信頼区間
  投稿者:aoki 2021/09/19(Sun) 22:48 No. 23142
気づくのが遅れましたが。

どの「ガイドライン」でしょうか?
URLなどをお示し下さい。
そのお方に直接質問するのがよいとは思いますが...

一般的には,統計量の分布について,特定の分布(例えば正規分布,t分布など)パラメトリック分布に従う場合は分布関数に基づいて信頼区間を計算します。あるいは経験的分布関数しかわからない場合でも,経験的分布における下側p%,上側q% などは計算できるので信頼区間は求めることができます(ブートストラップなどで)。

「対称信頼区間、最短非対称信頼区間」ググってご覧になるとよいかと思いますが。
といって,見てみると,全く見つかりませんね。

"asymmetry confidence interval proportion" などでググってみるとよいかも知れません。
いくつか出てきますが,読む気が起きません...歳ですね。

対称信頼区間は,「点推定±α」のように求められる信頼区間
非対称信頼区間は,「点推定-β, 点推定+γ」のように求められる信頼区間
また,非対称信頼区間はていぎにより何通りもありうるので,その中で信頼区間がもっとも狭いものを最短非対称信頼区間と言うのでしょう。
対称信頼区間は,
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/bohiritu-conf.html
の正規分布に近似する方法
非対称信頼区間は同じリンクの
F 分布に基づく方法(二項分布による方法)
さらには,ポアソン分布に基づく方法とか
さらにさらに,多種多様な信頼区間の定義とか。
そのなかで,信頼区間が最短のものが最短非対称信頼区間なのでしょう。




Re: 信頼区間
  投稿者:ビギナー 2021/09/21(Tue) 11:06 No. 23149
青木先生

ご多忙中のところ、ご回答いただき、本当にありがとうございました。
ガイドラインとは「後発医薬品の生物学的同等性試験ガイドラインに関する質疑応答(Q&A)について」のことでございます。以下にURLをお示しいたします。
https://www.pmda.go.jp/files/000234569.pdf

,亮遡笋砲弔ましてはガイドラインp13のQ-40の(A)に記載がございます。
△亮遡笋砲弔ましてはガイドラインp14のQ-42の(A)に記載がございます。

信頼区間とは正規分布又はt分布の性質を利用して、平均値±SEで求めるものと思っていましたが、ガイドラインでは「ノンパラメトリック法で解析した信頼区間を用いて・・・」と記載されておりましたので、分布を用いずに(中央値を用いる?)どのようにして信頼区間を求めることができるのかわかりませんでした。統計書を約40冊見ましたが、このことについて記載している本はございませんでした。対称信頼区間、最短非対称信頼区間の用語についても調べましたが、統計書に記載されておらず、ネットで検索しましたが、見つかりませんでした。
それで困り果てて青木先生にご質問したしだいでございます。
青木先生のご説明で対称信頼区間、最短非対称信頼区間の用語の意味がわかりました。
経験的分布関数という関数の存在を初めて知りました。経験的分布関数について、勉強いたします。統計の知識がまた増え、大変有難く思っております。
本当にありがとうございました。


c-indexの比較結果の解釈に関して
  投稿者:Boz 2021/09/20(Mon) 12:02 No. 23144
がん治療後の生存率予測に関する研究を行っています。3つの生存予測モデルを作成し、それぞれのc-index(C統計量)を算出したところ、下記の結果となりました。
モデルA:c-index=0.65(95%CI, 0.58-0.73)
モデルB:c-index=0.71(95%CI, 0.64-0.78)
モデルA+B:c-index=0.70(95%CI, 0.63-0.77)
※モデルA+BはモデルAとモデルBを組み合わせたものです。

続いて、各モデルのc-indexに有意な差があるかをRのパッケージ(compareC)を使って検討したところ、下記の結果となりました。
モデルA+B vs モデルA:P値=0.033
モデルB vs モデルA:P値=0.12

ここでお聞きしたいのは、モデルBの方がモデルA+Bよりも若干c-indexが高いにも関わらず(0.71と0.70)、モデルB vs モデルAには有意差が出ず、モデルA+B vs モデルAだけに有意差が出た理由の考察に関してです。『今回の母集団ではモデルBとモデルAのc-indexの差が、たまたま大きかった』という考察でよろしいでしょうか?

基本的な質問かもしれず恐縮ですが、ご教示いただけると幸いです。
よろしくお願い致します。

Re: c-indexの比較結果の解釈に関して
  投稿者:aoki 2021/09/20(Mon) 13:12 No. 23146
c-index の差の検定はよく知りませんが,一般的には,使用する変数の個数が違うモデル同士の比較は微妙ではないでしょうか?また,欠損値除去のせいで分析対象数が違う場合も関連するかも。
少なくとも,「たまたま大きかった」というようなことではないと思います。

Re: c-indexの比較結果の解釈に関して
  投稿者:Boz 2021/09/20(Mon) 14:40 No. 23148
青木先生

早々のご返信、深く感謝申し上げます。

>欠損値除去のせいで分析対象数が違う場合も関連するかも。
今回欠損値はなく、分析対象数はすべてのモデルで同じです。他に考え得る理由はございますでしょうか?

>一般的には,使用する変数の個数が違うモデル同士の比較は微妙ではないでしょうか?
これは、変数の個数が多いモデルの方が一般的に予測能が良くなりやすい、という意味でよろしいでしょうか?

>c-index の差の検定はよく知りませんが
Rのパッケージ(compareC)の元論文9ページ目を確認しましたところ、
z score =(c-indexの差)÷(c-indexの差の分散の平方根)を算出し、その値が分布の95%信頼区間外にあれば有意となるようです。
ですので、分子である『c-indexの差』は似たような値でも、分母である『c-indexの差の分散の平方根』に違いがあったため、今回2つのP値に違いが出たと予想しています。

元論文:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4314453/pdf/nihms643928.pdf


χ二乗検定の結果の説明法について
  投稿者:清水 2021/09/20(Mon) 00:05 No. 23143
お世話になっております。例えば、Aという転帰について、いくつかの項目があった場合に、χ二乗検定で有意と選択された因子1があったとします。
その場合に、転帰Aと因子1が関連があった、という書き方をしている人がいました。関連というのはあいまいで、避けたほうが良いと習った記憶がありますが、出典がはっきりしません。例えば連続変数動詞であれば、相関係数を求めて、「相関がある」と述べられると思いますが、このような場合にはどのように表記するべきでしょうか。

そもそも項目が複数個あった場合には、名義ロジスティック分析をするべきではないかという気がしてきました。
お忙しいところ恐縮です。何卒よろしくお願い致します。

Re: χ二乗検定の結果の説明法について
  投稿者:aoki 2021/09/20(Mon) 12:51 No. 23145
結論みたいなもの

英語ページ
クロス集計表レベルでは relationship, 連続変数では correlation
日本語ページ
クロス集計表レベルでは連関, 連続変数では相関
連関という用語は若干古くさい感じがするものの,関連よりは専門用語っぽい

連関は,線形相関,曲線相関も含め,あらゆる「関連性」を指すと思います

以下,参考

クロス集計表レベルでは relationship, 連続変数では correlation
https://bookdown.org/josiesmith/qrmbook/association-of-variables.html
https://www.displayr.com/what-is-a-crosstab/

relation, relationship
https://www.questionpro.com/cross-tabulation.html
https://humansofdata.atlan.com/2016/01/cross-tabulation-how-why/
https://www.qualtrics.com/experience-management/research/cross-tabulation/

タイトルではcorrelation, 本文中では relatonship
https://www.surveycake.com/en/featureinfo?f=cross-tabulation

以下のような記述が大勢

A crosstab is a table showing the relationship between two or more variables. Where the table only shows the relationship between two categorical variables, a crosstab is also known as a contingency table.

In this example, the two variables can both be viewed as being ordered. Consequently, we can potentially describe the patterns as being positive or negative correlations (negative in the table shown). However, where both variables are not ordered, we can simply refer to the strength of the correlation without discussing its direction (i.e., whether it is positive or negative).

連関 association 属性間に相互関係が存在することを表す言葉。関連という言葉も同じ意味で使われる。
連関表 contingency table
因果関係でなく,連関(association)関係としてとらえることが適切
連関関係の場合には行パーセント,列パーセントの両方で解釈できる
ファイ係数は,2個の質的変数の連関(相関)の強さを数値化した指標である
負の連関(関連),正の連関(関連)
強い連関 → 関連性がある 弱い連関 → 関連性がない

https://www.ibm.com/docs/ja/spss-statistics/23.0.0?topic=option-crosstabs
クロス集計表の説明に「相関」という用語は出てこない。「連関」が使われている。

> そもそも項目が複数個あった場合には、名義ロジスティック分析をするべきではないかという気がしてきました

その通りです。
クロス集計表の結果から変数を選択したりというのは誤りですが,
結果の解釈をする段階で,クロス集計表での所見が役に立つことはあります。

Re: χ二乗検定の結果の説明法について
  投稿者:清水 2021/09/20(Mon) 14:33 No. 23147
青木先生

 早速ご連絡ありがとうございます。非常に勉強になりました。英語論文ではrelationshipを使ってみたいと思います。
 今後ともどうぞよろしくお願い致します。

 清水拝

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