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

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

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

教育心理　重回帰　ロジスティクス
　　投稿者：佐藤 2019/07/01(Mon) 16:11 No. 22770
教育心理分野で、ある実習後の学習意欲を目的変数とし、実習中の環境要因２つ、実習生の個人要因２つを説明変数として重回帰分析をしようと考えていました。

学習意欲に対して正規性の検定を行ったところ正規性が認められなかったため、重回帰分析を適応できませんでした。

1. この場合、ノンパラメトリックの相関分析までしか適応できないのでしょうか？

2. 学習意欲のデータを見たところ二峰性があったため、学習意欲の中央値で2群に分割し、高学習意欲群1と低学習意欲群0にダミー変数化してロジスティクス回帰分析をしても良いのでしょうか？（それともこれはデータの恣意的な変更や情報量の削減という観点で許されないことなのでしょうか？）

統計に明るくなくおかしなことを聞いているかもしれませんが、ご教授頂けますと助かります。
Re: 教育心理　重回帰　ロジスティクス
　　投稿者：青木繁伸 2019/07/01(Mon) 21:10 No. 22771
> 学習意欲に対して正規性の検定を行ったところ正規性が認められなかったため、重回帰分析を適応できませんでした。

重回帰分析において，目的変数が正規分布するかどうかは関係ありません。というか，目的変数がどのような分布をするかによって，適切な分析方法があります。

例えば，一つ前の質問のような，目的変数が 0/1 の二値変数の場合にはロジスティック回帰分析とか，その他にも glm 関数は目的変数が様々な場合に対して分析を行うことができるように対処しています。
binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")
また，単に重回帰分析の場合であっても，従属変数が正規分布に従うことを要件にしてはいません。
Re: 教育心理　重回帰　ロジスティクス
　　投稿者：佐藤 2019/07/02(Tue) 14:18 No. 22772
青木先生

ご返信有難うございます。
不勉強でした重回帰分析は従属変数の正規分布を前提としていないんですね。
ネット上でK-S検定を行って正規性を確認してから、
重回帰分析を行っている例があったので、重回帰分析は正規分布を仮定している統計法と誤解をしていました。

従属変数は間隔尺度のデータなので、無理に2値化してロジスティクス回帰分析をするよりも、
重回帰分析を行なうほうがデータに沿った分析方法のように考えられます。

重回帰分析もう一度勉強し直してみます。

青木先生、今回は誠にありがとうございました。
Re: 教育心理　重回帰　ロジスティクス
　　投稿者：佐藤 2019/08/23(Fri) 11:48 No. 22812
青木先生

以前、丁寧にご回答頂きまことにありがとうございます。
以前の質問の続きなのですが、目的変数、説明変数ともに正規性の無いデータで重回帰分析を行なったところ、学会にて正規性の無いデータで重回帰分析は使用できないとの指摘を受けました。

「すくできる！リハビリテーション統計」という書籍や他医療統計の書籍に
「重回帰分析では独立変数も従属変数も正規分布した数量データ（間隔尺度か比率尺度）である必要があります。」

という記述があります。
以前のやり取りの中で、正規分布しないデータでも重回帰分析が適用できるという、
当方の理解が間違っていたのでしょうか？
Re: 教育心理　重回帰　ロジスティクス
　　投稿者：青木繁伸 2019/08/23(Fri) 14:21 No. 22813
その本の著者も，学会（査読者？）も間違えています。間違ったことを書いているクソページもたくさんありますね。

前に説明したことに基づいて，反論しなかったのですか？

説明変数は正規分布しなくても差し支えありません。ダミー変数を説明変数に使う場合を考えて見ればわかるでしょう。ダミー変数は 0/1 の二値データで，正規分布などしません。カテゴリーデータは重回帰分析に使えない（使うなら数量化 I 類を使いなさい）などという，とんでもない時代錯誤の考えを持っている人もまだいるのだろうか？また，実験計画によるデータの場合，説明変数は特定の値を取ります（薬剤投与量と効果の場合，薬剤投与量は例えば 10mg, 20mg, ..., 100mg などとされるでしょう。だれも薬剤投与量を正規乱数から選ぼうなどとは思わないでしょう）。

目的変数は，正規分布しなくても構いません。正規分布しないといけないのは残差の分布です（さらに，それさえも必須条件ではない）。
誤差が正規分布しないと，予測値の信頼区間を求めることができないからです。誤差（測定誤差も含む）は普通は正規分布に従うものです。なので，統計分析できるのです。

日本語のページは信頼できないということなら，英語のページを検索してはいかがでしょうか。

https://www.statisticssolutions.com/assumptions-of-multiple-linear-regression/
Multiple regression assumes that the residuals are normally distributed.

http://www.restore.ac.uk/srme/www/fac/soc/wie/research-new/srme/modules/mod3/3/index.html
Normally distributed residuals: The residuals should be normally distributed.

https://www.m3.com/open/clinical/news/article/604122/
回帰式の形は、被説明変数のタイプによって変えることができます。被説明変数が連続変量の場合は重回帰分析、2値変数（時間依存性がないイベントの発生または非発生、死亡または生存、治癒または非治癒など）の場合はロジスティック回帰分析、打ち切りのある2値変数の場合はCox回帰分析を用います。

重回帰分析では、残差（回帰分析による予測値と実測値の差分）の分布が正規分布である必要があります。

https://bellcurve.jp/statistics/course/9700.html
回帰モデルを考えるにあたって、誤差 μi にはいくつかの仮定している条件があります。
1. μi の期待値は0である
2. μiの分散は常にσ2 である
3. 異なる誤差 μi、μj は互いに独立である
これらの3条件から、μi は互いに独立に正規分布 N(0, σ2) に従うと仮定されます。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

具体例でも示してみましょう。正規分布すべきは「残差」です。

理論モデル y = 2.5 + 4.8 * x

重回帰モデルは y = a + b * x + ε
εが正規分布

データを用意する

n = 1000
x = runif(n, min=0, max=100) # 独立変数
EPSILON = rnorm(n, mean=0, sd=10) # ε：正規分布
y = 2.5 + 4.8 *x + EPSILON # 誤差を含む従属変数

y は正規分布しない！！

hist(y, main="図1")

予測は十分

plot(y ~ x, pch=19, col="#aa000055", alpha=0.3, main="図2")
ans = lm(y ~ x)
summary(ans)
abline(ans)
text(10, 400, sprintf("y = a + b x\na = %.3f\nb = %.3f", ans$coef[1], ans$coef[2]), pos=4)

誤差（予測誤差）は正規分布する

plot(ans$residuals ~ x, pch=19, col="#aa000055", main="図3")

hist(ans$residuals, main="図4")

図は，クリックすることで拡大表示



どのような統計をつかうのでしょうか
　　投稿者：こういち 2019/08/18(Sun) 03:28 No. 22802
青木先生
以下の疑問に関して教えてください

「疾患A型の患者」および「疾患B型の患者」における左室形態（4つタイプがありNormal, concentric remodeling(CR), concentric hypertrophy(CH), eccentric hypertrophy(EH))の分布を明らかにしました。結果の全体像は表1のようになりました。

ここで私が示したいことはA型ではCRが有意に多く、B型ではCHが有意に多い、
ということです。

そのために私がどのようにしたかというと
表2、表3のような2x2表をつくって（この作り方が正しいのか疑問ですが）
Fisher法で計算をしてそれぞれp=0.0427, p=0.0378という結果を得ました。
まず質問はこの統計処理が正しいのかを教えてください。

次にこれがもし正しいとしたら、下位比較をするという意味で、たとえば有意水準をBonferroniによって0.05/4C2=0.0083以下、とするのでしょうか？

別の考え方として4x2のクロス集計をおこなってP=0.1169を得たとして、
有意差なしでA型とB型の左室の型の分布に差がなし、で終了でしょうか？

以上、教えてください
よろしくお願いいたします

Re: どのような統計をつかうのでしょうか
　　投稿者：青木繁伸 2019/08/19(Mon) 10:11 No. 22803
最初から CR，CH 型の比較だけを行うと決めていたのなら，(CR, CH)×(A,B) の 2×2 分割表の検定をすればよいのですが，今回はそうではなさそうなので，全ての組み合わせでの 2×2 分割表の検定をしてボンフェローニ法を適用する。(CR, CR 以外)×(A,B)，(CH, CH以外)×(A,B)などとするのではないので注意。つまり，NormalとCR, NormalとCH，NormalとEH, CRとCH，CRとEH，CHとEH の 6 通りの比較を行う。
なお，「全体で有意である場合に下位検定として行える」という条件が必要かどうかは諸説ありますが，「この条件は不要」との見解が優勢かもしれない。

Fisher の正確検定でなければ，http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Tukey.html に書いてある方法を使えるでしょう（http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Ryan.html でもよいが，前者の法が通りがよいかな？）。
Re: どのような統計をつかうのでしょうか
　　投稿者：こういち 2019/08/20(Tue) 00:56 No. 22804
青木先生

お返事をいただきありがとうございました。
responseが遅くなりお詫び申し上げます。
先生からいただいた説明を何度も何度も読み直してみましたが
私の説明が悪くて質問の意図が伝えられていないと思いました

私の疑問は
「たとえばＡ型という疾患においてCR, CH型の頻度を比較をする」ということではなく
「4つある形態の一つであるCRという形態を呈する患者が、Ａ型では全体の56%, B型では全体の38%で、これらを比較するときに、どの様な統計を用いれば、CRという形態はA型がB型よりも有意に多いことを示せるか？」というA型とB型の比較のための統計処理がわからない、ということです。

そして同様に上記と独立して、
「4つある形態の一つであるCHという形態を呈する患者が、Ａ型では全体の34%, B型では全体の52%で、これらを比較するときに、どの様な統計を用いれば、CHという形態はB型がA型よりも有意に多いことを示せるか？」というA型とB型の比較ののための統計処理がわからない、ということです。

私のいっている比較そのものがそもそも統計処理にそぐわず、有意に多い、ことは示せないのでしょうか？。
基本的なことが分かっていないとんちんかんな質問かもしれないと、
だんだん心配になってきましたが、どうか教えていただければ、本当にありがたいです。
よろしくお願いいたします。
Re: どのような統計をつかうのでしょうか
　　投稿者：青木繁伸 2019/08/20(Tue) 08:28 No. 22805
No. 22802
> 私が示したいことはA型ではCRが有意に多く、B型ではCHが有意に多い、ということです。

と

No. 22804
> 「たとえばＡ型という疾患においてCR, CH型の頻度を比較をする」ということではなく...CRという形態はA型がB型よりも有意に多いことを示せるか

相反する説明がありますが。まあ，それは行方向のパーセントと列方向のパーセントのどちらを評価対象にするかということで，どちらが原因で，どちらが結果なのかで自ずと決まるものですが。でも，chisq.test でも，Fisher.test でも行・列の違いは無視されます（結果は同じ）。

> 「4つある形態の一つであるCRという形態を呈する患者が、Ａ型では全体の56%, B型では全体の38%で、これらを比較するときに、どの様な統計を用いれば、CRという形態はA型がB型よりも有意に多いことを示せるか？」というA型とB型の比較のための統計処理

ならば，二項検定でよいでしょう。これだと，行ごとの検定になります。CR が A:66, B:19 で帰無仮説では母比率=0.5
> binom.test(66, 85)

	Exact binomial test

data:  66 and 85
number of successes = 66, number of trials = 85, p-value = 3.041e-07
同様に CH では

> binom.test(40,66)

	Exact binomial test

data:  40 and 66
number of successes = 40, number of trials = 66, p-value = 0.1089
ただし，この二つの検定は独立ではないです（独立性の検定という意味ではない）。
CR の 66 も 29 も Normal, CH, EH の度数による制約があります。

(Normal,CR,CH,EH)x(A,B) はそれらの制約のもと，Normal, CR, CH, H それぞれで，A, B の割合に違いがあるかを検定するものですよね。
検定ファミリーの一つに (CR,CH)x(A,B) があり，それは，CR, CH で A, B の割合に違いがあるかを検定する。
> fisher.test(matrix(c(66, 40, 19, 26), 2))

	Fisher's Exact Test for Count Data

data:  matrix(c(66, 40, 19, 26), 2)
p-value = 0.03111
CR では A，CH では B が有意に多い（A では CR，B では CH が有意に多い）といえるわけでしょう。

なお，また最初に戻りますが，「A型ではCRが有意に多く、B型ではCHが有意に多い」ということなら，A, B ごとに一様性の検定
> chisq.test(c(7, 66, 40, 5))

	Chi-squared test for given probabilities

data:  c(7, 66, 40, 5)
X-squared = 86.407, df = 3, p-value < 2.2e-16

> chisq.test(c(2, 19, 26, 3))

	Chi-squared test for given probabilities

data:  c(2, 19, 26, 3)
X-squared = 34, df = 3, p-value = 1.981e-07
でも，これはこの場合は Normal と EH がCR, CH より格段に少ないことが原因。
> binom.test(66, 106)

	Exact binomial test

data:  66 and 106
number of successes = 66, number of trials = 106, p-value = 0.01478

> binom.test(19, 45)

	Exact binomial test

data:  19 and 45
number of successes = 19, number of trials = 45, p-value = 0.3713
こんないろいろなことを考えなきゃいけないけど (CR,CH)x(A,B) でやればよいのでは？
Re: どのような統計をつかうのでしょうか
　　投稿者：こういち 2019/08/21(Wed) 21:43 No. 22806
青木先生

丁寧に教えていただきまして、本当にありがとうございます。
教えていただいたことが恐らくは十分に理解されていないと思いますので
時間をかけて考えてみようと思います。

最後にひとつだけ教えてください、
先生が薦めてくださる「(CR, CH)x(A, B)をやる」、というのは
下のような２ｘ２の表を作成してp=0.0311を求め、
（この表の作り方が正しいのかわかりません）
この場合「A型ではCRが多く、B型ではCHが多い」
（←この言い方が正しいかどうかにも自信がありません）ことが
0.0311>0.008(=0.05/6)であり有意とは言えない

という理解でよろしいでしょうか？
最後に教えていただければ本当にありがたいです

Re: どのような統計をつかうのでしょうか
　　投稿者：青木繁伸 2019/08/22(Thu) 08:13 No. 22807
> A型ではCRが多く、B型ではCHが多い
> 0.0311>0.008(=0.05/6)であり有意とは言えない

でよいと思います。
Re: どのような統計をつかうのでしょうか
　　投稿者：こういち 2019/08/22(Thu) 18:18 No. 22811
青木先生

お忙しい中いろいろ教えていただき
本当にありがとうございました
また、よろしくお願いいたします


3Dグラフの作成方法について
　　投稿者：さくらとお城 2019/08/09(Fri) 14:19 No. 22791
青木先生

　いつもお世話になっております。

　複数のチャンネルで計測した経時データを，plot3Dなどで図示しようとしているが，なかなかうまく行きません。

　データとエクセルで作成したグラフは添付ファイルのようになります。実際は，もっとチャンネルが多いので，ここでは2つだけ表示しています。

　きれいに図示できる方法を教えてください。

Download:22791.zip 22791.zip

Re: 3Dグラフの作成方法について
　　投稿者：青木繁伸 2019/08/09(Fri) 20:20 No. 22792
3次元グラフはあまりよい表示方法ではないでしょう。
添付する図を見ても，二次元グラフ（折れ線グラフ）のほうが二つのデータの相違・同等がよく分かるでしょう。
確かにたくさん重なるとわかりにくいでしょうが，3次元折れ線グラフが重なったら，なにがなんだかわからなくなります。

Re: 3Dグラフの作成方法について
　　投稿者：青木繁伸 2019/08/09(Fri) 20:30 No. 22793
ggplot の例によくあるような，縦にずらっと折れ線グラフを並べるのはどうですか？
一つのグラフに，傾向の似通った複数の折れ線グラフを描けばなおよいかも。

Re: 3Dグラフの作成方法について
　　投稿者：さくらとお城 2019/08/10(Sat) 10:56 No. 22794
青木先生

　ご返事ありがとうございます。添付のようなグラフを考えていました。ちなみに，これが別ソフトで作ったもので，細かい設定などが大変で，Rでもできるかなと思って質問しました。

　先生のggplotの例も良さそうで，試してみます。

Re: 3Dグラフの作成方法について
　　投稿者：さくらとお城 2019/08/22(Thu) 10:22 No. 22808
青木先生

　その後，色々と調べたけど，恥ずかしながら，なかなか先生のような図を描けません。ggplotの例を書いたコードを見せて頂いて宜しいでしょうか？

　どうぞよろしくお願い致します。
Re: 3Dグラフの作成方法について
　　投稿者：青木繁伸 2019/08/22(Thu) 12:50 No. 22809
私自身，ggplot は使わないので，細かな設定は知りません。
以下のようなプログラム例を示しておきますので，追加，修正の上ご利用ください。
ポイントは通常のデータフレームを時間とチャンネルとデータ値の3列から成るデータフレームに変換して，それを用いて折れ線グラフを描くということです。
set.seed(123)

m = 10
n = 50

# n x m 行列：列ベクトルがチャンネルごとの時系列データ

simdata0 = matrix(0, n, m)
for (i in 1:m) {
	simdata0[, i] = cumsum(rnorm(n))
}

# チャンネルごとのデータを，
# time:時間（1:n），value:データ値，チャンネル名の，n x 3 のデータフレームとして rbind で縦に繋いだデータフレームに変換

simdata = NULL
for (i in 1:m) {
	simdata = rbind(simdata, data.frame(time = 1:n, value = simdata0[, i], 
		channel = paste0("x", i)))
}

g = ggplot(simdata, aes(x = time, y = value))
g = g + geom_line(aes(colour = channel), size = 0.5)
g = g + geom_point(aes(colour = channel), size = 1, alpha = 0.5)
g = g + facet_wrap(~channel, ncol = 2) # ncol で何列のグラフにするか
g = g + theme_classic() # または g = g + theme_bw() でも
g
以下のようなグラフになります（部分）

Re: 3Dグラフの作成方法について
　　投稿者：さくらとお城 2019/08/22(Thu) 12:53 No. 22810
青木先生

　早速のご返事ありがとうございます。先生のコードを勉強させていただきます。


統計分析手法の選択につきまして
　　投稿者：こだま 2019/08/14(Wed) 13:13 No. 22795
青木先生

はじめまして。

自身の周りに統計手法を専門的に理解している人間があまりおらず書き込みさせて
いただきます。

私は現在、大学院生（博士後期課程）に在学中で現在GISを用いた緑地の分布に
ついての分析結果学を会論文投稿に向け作業をしているところです。

その中で、現在の緑地の分布要因は何に規定されているのか？を統計手法による
分析を用いて明らかにしたいと考えております。

その際にどの分析手法が良いのかに行き詰っているところです（判別分析なのか、
数量化2類なのかと悩んでいます。）。

データはメッシュで集計しており，メッシュ毎に以下のデータを格納しています。

目的変数：緑地の分布タイプ（別途分析により7タイプに分類済）
説明変数：用途地域（13タイプ），区画整理事業の有無，傾斜角度，標高，駅の有無等

目的変数が質的データ、説明変数が質的データと量的データが混在しているため、単純
に判別分析や数量化2類を使えない？のかと思っています（理解不足なのかもしれませんんが）。

お忙しいところお手数をおかけしますが、どうぞご回答いただけますと幸いです。
Re: 統計分析手法の選択につきまして
　　投稿者：青木繁伸 2019/08/14(Wed) 13:20 No. 22796
> 目的変数が質的データ、説明変数が質的データと量的データが混在しているため、単純
に判別分析や数量化2類を使えない？のかと思っています

の点についてのみ。
質的データはダミー変数として取り扱うことで，量的データと一緒にして判別分析に使用できます。（ちなみに，数量化II類はダミー変数を使う判別分析と等価です）
Re: 統計分析手法の選択につきまして
　　投稿者：こだま 2019/08/14(Wed) 13:58 No. 22797
青木先生

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

判別分析を採用してみたいと思います。

標高等の量的なデータを区分して数量化2類かな？などと考えて
いて、その区分の方法が難しいなと思っていたところでした。

今後、また相談させていただくことがあるかもしれませんが、
どうぞよろしくお願いいたします。
Re: 統計分析手法の選択につきまして
　　投稿者：こだま 2019/08/14(Wed) 17:12 No. 22798
青木先生

度々申し訳ございません。

色々な文献（多変量解析の本やインターネット、判別分析を行っている研究論文）を
調べているのですが、多くが目的変数が2択（例えば、公園がある/ない のような者）
の事例でした。

今回、分析しようとしているように目的変数が7つのタイプがあっても分析手法として
採用しても問題はないのでしょうか。

可能であればで結構ですので、ご教授いただけますと幸いです。
Re: 統計分析手法の選択につきまして
　　投稿者：青木繁伸 2019/08/14(Wed) 21:32 No. 22799
ロジスティック判別も線形判別も，3群以上であっても判別はできます。
どこでもよいですが，判別分析のページを読んでみてください。
一番簡単な2群判別の話しか書いていないのは，くそページです。

7群判別ともなると理論的には可能ですが，実際に分析して得られる解の解釈は難しいかも。（貴方の用意した説明変数で十分に判別できるかどうかということですが）
Re: 統計分析手法の選択につきまして
　　投稿者：こだま 2019/08/14(Wed) 23:06 No. 22800
青木先生

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

論文提出まで時間があるので判別分析含め、その他の多変量解析について
色々な参考書やページを読んで勉強いたします。

説明変数の選択が難しいところですが、こちらは既往文献等を参考にしつつ
選択していきたいと思います。

実際に分析を進めてみて何かありましたら相談させていただきたいと思います。

夜分遅くまでお手数をおかけしました。


【R】リストに格納された複数データフレームのEXCEL保存
　　投稿者：明石 2019/08/05(Mon) 16:12 No. 22787
青木先生　様；

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

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

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

リストに格納された複数データフレームを、
拡張パッケージ openxlsx を用いて、
EXCELの複数のシートに格納したxlsxファイルで保存することを考えます。

以下、irisデータでお示しをします。

data(iris)

dfs <- split(iris, iris$Species)

library(openxlsx)　　#　拡張パッケージ

write.xlsx(dfs, file="iris_Species.xlsx")

setosa、versicolor、virginica　の３つのシートからなる
"iris_Species.xlsx"ファイルは保存されますが、以下のメッセージがでます。

「Note: zip::zip() is deprecated, please use zip::zipr() instead」

私が書いたコマンドには、zip()関数はありませんので、
このようなメッセージが表示された場合、どのように対処すればよいのでしょうか。

ご教示をいただけましたら、大変に助かります。
お手数をおかけいたします。
Re: 【R】リストに格納された複数データフレームのEXCEL保存
　　投稿者：青木繁伸 2019/08/05(Mon) 19:40 No. 22788
write.xlsx から呼ばれる saveWorkbook の中で zip を呼ぶようですが（明示的には見えない），私にはそれ以上追求できませんでした。
deprecated というのは，将来的に廃止されるということで，無効になっているわけではないです。ということで，しばらくはそのまま使えば良いと思います。そのうち，使えなくなったらまずは openxlsx パッケージのバージョンを確認して，新しくなっていたらインストールし直すということになるでしょう。それもできないという事態になったら，パッケージ開発者に連絡するしかなくなるかな。
Re: 【R】リストに格納された複数データフレームのEXCEL保存
　　投稿者：青木繁伸 2019/08/05(Mon) 19:43 No. 22789
Excel シートにする必要性がどの程度あるのか，他の人と共有するためには Excel シートでなければならないとか，いきなり Excel を書き込むのではなく，CSV ファイルとして書き出すのではだめなのか，とか，私なら Excel を使わない方向で行きたいですね。
Re: 【R】リストに格納された複数データフレームのEXCEL保存
　　投稿者：明石 2019/08/06(Tue) 08:05 No. 22790
青木先生　様；

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

今回も有難いご教示をいただき、大変に助かりました。

青木先生がお書きになられたことを、戒めたいと思います

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


【R】リストに格納された複数データフレームのファイル保存（ループを使わない）
　　投稿者：明石 2019/07/29(Mon) 17:25 No. 22784
青木先生　様；

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

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

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

リストに格納された複数データフレームのファイル保存を、
ループを使わない書き方について、ご教示いただければ大変に助かります。

以下、例示いたします。

データフレーム dat を、カラム名"支店名"でファイル分割したとします。

dfs <- split(dat, dat$支店名)

支店名ごとに分割されたデータフレーム群を、csvファイル保存します。

n <- length(dfs)
names <- names(table(dat$支店名))

for (i in 1:n) {
　fn <- paste(names[i],".csv", sep="")
　write.csv(dfs[[i]],fn, row.names=F)
}

ループを回さないで、apply関数群などを用いて、シンプルに書けるようならば、
ご教示をいただきたいと思います。

お手数をおかけいたします。
Re: 【R】リストに格納された複数データフレームのファイル保存（ループを使わない）
　　投稿者：青木繁伸 2019/07/29(Mon) 22:16 No. 22785
iris データセットでやってみましょう

Species で 3 分割

iris2 <- split(iris, iris$Species)

書き出す。junk への代入は，ゴミを表示させないためだけ。

junk <- sapply(names(iris2), function(n) write.csv(iris2[[n]], paste(n, ".csv", sep="")))

あまりわかりやすいものではないですね。
貴方が示した for ループを使うほうがずっと分かりやすい。
Re: 【R】リストに格納された複数データフレームのファイル保存（ループを使わない）
　　投稿者：明石 2019/07/30(Tue) 06:13 No. 22786
青木先生　様；

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

今回も、大変に良い勉強をさせていただきました。

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


因子分析の因子の差について
　　投稿者：コロン 2019/07/16(Tue) 13:38 No. 22778
いつもお世話になっております。

早速ではございますが，「4つの因子」が出てきたという例で説明させて頂きます。

数学力が高いもの，中程度のもの，低いものという3つのグループの間で，因子において差があるかどうかを見る場合なのですが...

（1）因子毎に，一要因分散分析（被験者間）を行い，差があれば多重比較で見ていく
（2）4つの因子をすべて用いて，二要因分散分析（被験者「間」×被験者「内」）を行い，必要に応じて下位検定を行う

どちらが正しいのでしょうか？古い論文を見ていると（1）が圧倒的に多いのですが...。

よろしくお願いいたします。
Re: 因子分析の因子の差について
　　投稿者：青木繁伸 2019/07/17(Wed) 09:06 No. 22779
理論的には (2) が正しい（妥当）と思います

なぜ (1) が圧倒的に多いかは，いくつか理由があるでしょう

サンプルサイズが小さく，(2) に耐えない
モデルが複雑になり，解釈が難しくなる
コンピュータプログラムがなく（利用しにくく），計算できない
「だって，みんな (1) をやっているんだもの」と長いものに巻かれる
そもそも (2) という選択肢を知らない
Re: 因子分析の因子の差について
　　投稿者：コロン 2019/07/17(Wed) 09:36 No. 22780
青木先生

お忙しい中，いつもありがとうございます。

やはり，そうですよね....。かなり検索をしてみたのですが，すべて（1）で処理されていたのです。（2）が適切ではないかなとずっと考えていましたので，すっきり致しました。

ありがとうございました。
Re: 因子分析の因子の差について
　　投稿者：コロン 2019/07/20(Sat) 17:19 No. 22781
青木先生

もう一つ重要なことを聞き忘れていました。

各因子の項目の数が異なるのですが、尺度得点あるいは因子得点を利用して２要因分散分析を行ってもよろしいでしょうか。

よろしくお願いいたします。
Re: 因子分析の因子の差について
　　投稿者：青木繁伸 2019/07/24(Wed) 11:04 No. 22782
項目の数が異なるというのが二元配置分散分析の障害になるのですか？
Re: 因子分析の因子の差について
　　投稿者：コロン 2019/07/25(Thu) 16:34 No. 22783
青木先生

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

【要因１】
因子１　７項目
因子２　５項目
因子３　４項目

【要因２】
学力上位群
学力下位群

項目数が異なるため，各因子の平均値をそのまま解釈して良いのかなと悩んでおりました。ひょっとしたら，過去の研究，項目数が異なるために因子毎に比較していたのかなと考えたり...。

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


固辞，数字の置換
　　投稿者：コロン 2019/07/08(Mon) 13:44 No. 22773
いつもお世話になっております。

基本的な質問で申し訳ございませんが，以下のデータがあるとします。

class score
1 10
1 11
2 13
2 14
3 12
3 15
4 17
4 14
5 13
5 10

classの 3 を 5 に修正して新しいデータフレーム形式で表示させるにはどのようにすればよろしいのでしょうか？
Re: 固辞，数字の置換
　　投稿者：青木繁伸 2019/07/08(Mon) 15:37 No. 22774
データフレームを a とすれば，
a$class[a$class==3] <- 5
で置換できます。
Re: 固辞，数字の置換
　　投稿者：コロン 2019/07/08(Mon) 16:40 No. 22775
青木先生

お忙しい中申し訳ございません。

ネットで検索をしていたら，置換は gsub というサイトが多かったため，gsubでやっていたのですが，何度やっても，データフレーム形式にならず，かつ，すべての数字が５に変換されたので，わからなくなり，投稿した次第です。

ありがとうございました。
Re: 固辞，数字の置換
　　投稿者：青木繁伸 2019/07/08(Mon) 18:07 No. 22776
gsub でもできないことはないですが，gsub は本来は文字列の置換ですから
a$class <- as.numeric(gsub(3, 5, a$class))
as.numeric してやるなんて余計なことをしないといけません。
Re: 固辞，数字の置換
　　投稿者：コロン 2019/07/08(Mon) 18:24 No. 22777
青木先生

as.numericなんですね！gsubでもできました。ただ先生から最初に提示いただきましたものがわかりやすいです。

なんどもありがとうございました。


hclust()の距離行列とメソッドの指定
　　投稿者：波音 2019/06/28(Fri) 15:03 No. 22764
Rのhclust()について教えてください。

引数にward.Dとward.D2が指定できるようになっていますが、例えば「d」という距離行列(ユークリッド距離)が存在したと時にいくつかの方法がネット上や書籍で示されています。

・距離行列をそのまま指定する
hclust(d, method="ward.D")

・距離行列の自乗したものを指定する
hclust(d^2, method="ward.D")

・川端ほか「Rにおる多変量解析入門」オーム社（2019）P374記載
hclust((1/2)*d^2, method="ward.D")

結局のところ「解釈ができれば何でもいい」というくらいの問題なのかもしれませんが、アルゴリズム的（？）にはどういう指定が適切なのでしょうか。

川端ほか「Rにおる多変量解析入門」オーム社（2019）P313にward.Dとward.D2の違いが説明されていますが、どういう距離行列とメソッドの組み合わせ指定をすればよいのでしょうか、、、
Re: hclust()の距離行列とメソッドの指定
　　投稿者：青木繁伸 2019/06/29(Sat) 12:05 No. 22765
川端さんが具体的にどのように説明されているのか分かりませんが，私の理解は
http://aoki2.si.gunma-u.ac.jp/lecture/misc/clustan.html
に示してあります。これは奥野忠一らの多変量解析，多変量解析入門にあったものだと思います。まだパソコンなんてものがなく，コンピュータはFORTRAN,COBOLでプログラムするという時代のものです。もともとは，CLUSTAN というプログラムのアルゴリズムの説明にあるものです。実際，私もこれに基づいてクラスター分析のプログラムを書いて使っていました。
当初の R には今の ward.D が ward として実装されていました。推測するに，「バグ」だったのだと思います。
私の書いた本でも，「ウォード法と重心法を採用する場合にはユークリッド平方距離を使うべきである。ユークリッド平方距離を選択できるオプションがないので，後述の hclust の第 1 引数としてユークリッド距離を二乗したものを指定することを忘れないようにしよう。」と書きました。
その後指摘があったかなんかで，本来あるべき ward.D2 が実装されましたが，過去のしがらみを捨てきれなくて（バグだったとも言えなくて）両方を残したのだと邪推しています。
Re: hclust()の距離行列とメソッドの指定
　　投稿者：波音 2019/07/01(Mon) 10:41 No. 22769
ご回答ありがとうございます。

たしかに青木先生の本にも「hclust の第 1 引数としてユークリッド距離を二乗したものを指定」と書かれており、それ以前にもこのように教わっていたので、私の中ではある意味ルール的にd^2としていましたが、、、

やはり（使う人の立場によっては程度問題ではあると思いますが）アルゴリズムを理解しておくことは重要だと感じさせられました。

少し話はそれてしまいますが、誰でも使えるものというのは「考えなくても（理解していなくても）使えるもの」ともいえるので、ただでさえ色々と誤解の多い統計話題は今後さらに注意しないといけないと思いました。

※パッケージ文化に恵まれた我々は特に注意せねばならないのかもしれません。


ダミー変数を従属変数にした場合の切片
　　投稿者：山本 2019/06/29(Sat) 23:47 No. 22766
はじめて書き込みます。山本です。
講義の関係で重回帰分析を行なっているのですが、ダミー変数を従属変数とした際に切片が負という結果がでました。
ダミー変数は0か1かですので正になるはずなのですが…なぜでしょうか
この結果はどう解釈すればよいのでしょうか
それとも説明変数を変更する必要があるのでしょうか
お願いいたします
Re: ダミー変数を従属変数にした場合の切片
　　投稿者：青木繁伸 2019/06/30(Sun) 07:19 No. 22767
従属変数が 0/1 の二値データであることと，切片が負の値を取ることは何の矛盾もありません。
というか，独立変数のとる値により切片は負の値になることも，[0, 1] の値になることも，1より大きな値になることもあります。
たとえば，
d <- data.frame(y = c(1, 0, 1, 1, 0, 0, 0, 1, 0, 0),
	x1 = c(58, 37, 38, 50, 36, 59, 49, 66, 47, 58),
	x2 = c(54, 53, 50, 46, 46, 65, 37, 69, 37, 44))
の場合は，切片は負の値になります。
> lm(y~x1+x2, data=d)
(Intercept)           x1           x2  
  -0.646540     0.005094     0.015825
ここで，x1 から 250 を引いた変数 x1.250 を作り，それを用いて予測式を計算すると
> d$x1.250 <- d$x1 - 250
> lm(y~x1.250+x2, data=d)
(Intercept)       x1.250           x2  
   0.627073     0.005094     0.015825
ここで，x1 から 435 を引いた変数 x1.435 を作り，それを用いて予測式を計算すると
> d$x1.435 <- d$x1 - 435
> lm(y~x1.435+x2, data=d)
(Intercept)       x1.435           x2  
   1.569546     0.005094     0.015825
それぞれの偏回帰係数は全く同じになっていることに注意。
つまり，独立変数の平行移動をした場合，その影響を打ち消すために定数項が調整されるということ（ちなみに独立変数を定数倍したときは偏回帰係数が定数分の1になる）。

で，今回の問題の本質は，なぜ従属変数が 0/1 の二値データなのに予測値が0より小さくなったり，1より大きくなったりしてしまうのかということ。

それは直線回帰をするからやむを得ないということです。
図を見れば明らかです。

二値データはロジスティック回帰分析を行うのが定石です。
図の赤で示したのがロジスティック曲線です。
> glm(y~x1+x2, data=d, family=binomial(link="logit"))

Call:  glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
    data = d)

Coefficients:
(Intercept)           x1           x2  
   -5.08999      0.02083      0.07214  

Degrees of Freedom: 9 Total (i.e. Null);  7 Residual
Null Deviance:	    13.46 
Residual Deviance: 11.93 	AIC: 17.93

Re: ダミー変数を従属変数にした場合の切片
　　投稿者：山本 2019/06/30(Sun) 12:02 No. 22768
素早く回答していただきありがとうございます。
そういった分析もあったのですね…
よく分かりました


介入研究の因子分析
　　投稿者：春 2019/06/13(Thu) 00:34 No. 22746
青木先生

いつもお世話になっております。
また、教えていただきたいことがあり、ご連絡させていただきました。

現在、先行研究で使われていた尺度（20項目ほど）を用いて調査を行いました。
その先行研究では因子分析を行っていませんでしたが、今回の調査は介入研究なので、尺度を因子分析を行った方が傾向を示しやすいと思うのですが、介入研究の場合の因子分析は、どのように行えば良いのでしょうか？
ちなみに、介入前と介入後の得点で、それぞれ因子分析しましたが、同じように分かれませんでした。

ご助言いただけますでしょうか。
よろしくお願いいたします。
Re: 介入研究の因子分析
　　投稿者：春 2019/06/17(Mon) 15:24 No. 22758
青木先生

いつもお世話になっております。

先日、6月13日にご質問させていただいた者です。
どうも、タイミング悪く入ってしまったようで、お返事を、まだ、いただけておりません。
ぜひとも、ご助言いただけないでしょうか？

よろしくお願いいたします。
Re: 介入研究の因子分析
　　投稿者：青木繁伸 2019/06/17(Mon) 19:45 No. 22759
上の方にあるものばかり見ていたようで，気づきませんでした。

さて，先行研究で使われていた尺度は因子分析で作成されたものではないと言うことでしょうか。

とすると，まずいですね。あなたの調査データで一因子性がないという結果になっても，先行研究を否定するということになってしまいますし，対象が違うから違う因子が得られたのだといわれるかも知れませんし。

介入研究だからといって，因子分析のやり方に特徴や注意点があるとは思えませんが。

介入前後で因子が異なることは想像するに難くないと思いますよ。

実際の因子分析の結果だけでなく，基礎的な統計解析を見ずにああだこうだいうのは難しいです。もう少しはっきりいえば，あなたの持っているデータを自由にいろいろ分析してみてその結果を総合して何かをいうということでなければ（貴方の分析結果を見ただけではなく），なんもいえない。ということです。そういうことをするのは，共同研究者以上の立場ですよね。査読者はある意味そのような立場にあるのでしょうが，実際に「データを見せろ，私が解析してみて貴方の言うことが正しいかどうか示してやる」とは言えない訳です。

データ解析というものは，そういうものです。
Re: 介入研究の因子分析
　　投稿者：青木繁伸 2019/06/17(Mon) 21:33 No. 22760
「十分な分析を行い，以下のような結果を得ました」というとき，どのように十分な分析だったかを「詳細に」示し，それを査読者が十分だったと認めるかどうかということです。

蛇足ではありますが，査読者が十分だったと認め論文掲載となっても，読者（第三者）の十分だったかどうかの判断とは必ずしも一致しない（投稿論文への異議申し立てはあり得る）。

その昔，ある学会誌に「奇病」の投稿があったが，老臨床医が「それは，脚気だろ」と指摘してチャンチャン！となったことがある（そんなに大昔のことではない）

厳しく言えばそういうこと。
Re: 介入研究の因子分析
　　投稿者：春 2019/06/18(Tue) 00:41 No. 22761
青木先生

ご助言ありがとうございます。

そもそも、因子分析されていない尺度を使う事自体、無理があるということですね。
もう1つ、気がかりなことがあります。Nの数が30に満たないので、それも、因子分析するには、少なすぎるかなあと思っておりました。

このまま、因子分析を行うことは、諦めた方が良いのかとも思いますが、以下のように、進めるのは、間違いでしょうか？

介入前の因子分析の結果では、4因子に分かれ、α係数も1つだけ0.7に足りませんでしたが、3つは0.8以上あり、内的整合性は高いと判断できます。
介入後も、介入前に抽出された因子に合わせて、α係数を出して充分な値であれば、その後に介入前と介入後で対応のあるｔ検定をする。

・・・というのは、分析としてはどうでしょうか？
ご意見を伺えると助かります。
Re: 介入研究の因子分析
　　投稿者：青木繁伸 2019/06/18(Tue) 08:38 No. 22762
> Nの数が30に満たないので、それも、因子分析するには、少なすぎるかなあと思っておりました。

そこは最初に書いておいてくれなくては。
圧倒的に少なすぎます。後で述べる対応のある t 検定であっても，十分とは言いがたい。

> その後に介入前と介入後で対応のあるｔ検定をする。

何の対応のある t 検定ですか？先行研究の尺度の合計値ですか？それとも，4因子それぞれに含まれる項目の合計点ですか？

そもそも，因子分析の結果の4因子はそれぞれ何項目が含まれるのか，それぞれの因子の意味づけ（解釈）はできるのか。

もともと20項目あるなら，項目別に検定を行ったらどうなるのか？
Re: 介入研究の因子分析
　　投稿者：春 2019/06/18(Tue) 12:13 No. 22763
青木先生

おはようございます。
ご回答いただいて、ありがとうございます。

項目別の検定は、既に行っていて、効果は確認済みです。
しかし、ご指摘いただいたように、Nの数が少ないので、今回は因子分析は諦めようと思います。
しかし、今後、Nを増やして再度、分析を行いたいと思っているので、このまま、続けさせていただきたいと思います。

対応のあるt検定については、抽出した各因子の下位尺度得点として、介入前、介入後の各平均値を算出し、比較をしようと思っています。

しかしながら、今後、因子分析を行うときには、各因子内の項目の数や、各因子の意味づけ（解釈）に気を付けていきたいと思います。
そして、因子分析を用いた方が良いのか、項目別に検定だけの方が良いのか判断しようと思います。

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


【R】do.call()関数を用いたリスト（複数のデータフレーム）のmerge
　　投稿者：明石 2019/06/16(Sun) 09:50 No. 22754
青木先生　様；

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

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

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

リスト（複数のデータフレーム）のmergeを、forループで順次mergeするのではなく、
do.call()関数を用いて、もし書けるようならば、ご教示をお願いいたします。

do.call(rbind, list)　にならい、
do.call(merge, list)　のように、もし書けたとしたら、どのように書けるのか興味があります。

# tempフォルダ下にある、csvファイルのファイル名を取得

path <- setwd("C:/temp")
names <- list.files(path, pattern="csv")
nn <- length(names)

d <- list(nn)

loop <- 1: nn

for (i in loop) {
d[[i]] <- read.table(names[i], header=TRUE, sep=",", stringsAsFactors=F)
}

#　リスト d に格納されている nn個のデータフレームは、いずれも同じ、mergeの属性 "支店名"をもっています。
#　要素の添字に沿って、完全外部結合を順次、繰り返したいと思います。

for (i in loop) {
if( i == 1) {
dat <- d[[1]]
} else {
dat <- merge(dat, d[[i]], by.x="支店名", by.y="支店名" , all.x = TRUE, all.y = TRUE ,no.dups = F)
}
}

もし、do.call(merge, d) のように書けるようならば、ご教示をいただきたいと思います。

お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。
Re: 【R】do.call()関数を用いたリスト（複数のデータフレーム）のmerge
　　投稿者：青木繁伸 2019/06/16(Sun) 12:34 No. 22756
別に do.call でなくても，簡単に書きたいということなんでしょうか？

Reduce を使うのはいかが？
set.seed(11111)
d <- vector("list", 3)
d[[1]] <- data.frame(a=1:5, b=letters[1:5])
d[[2]] <- data.frame(a=c(1,3,5,6,8), c=LETTERS[1:5])
d[[3]] <- data.frame(a=c(2,4,5,7,9), d=rnorm(5))

merge2 <- function(x, y) merge(x, y, all.x = TRUE, all.y = TRUE, no.dups = FALSE)  
Reduce(merge2, d)
一行で書くと
Reduce(function(x, y) merge(x, y, all.x = TRUE, all.y = TRUE, no.dups = FALSE), d)
なお，for を使う場合も，以下のようにすれば十分簡潔に書けるのでは？
x <- d[[1]] # 結果の初期化は，最初のデータフレームを代入
for (i in 2:length(d)) { # 2 番目以降のデータフレームをマージ
  x <- merge(x, d[[i]], all.x = TRUE, all.y = TRUE, no.dups = FALSE)
}
x
どのように書いても，以下のような結果になると思いますが。
  a    b    c            d
1 1    a    A           NA
2 2    b <NA>  0.003630422
3 3    c    B           NA
4 4    d <NA>  0.798122717
5 5    e    C  0.893397342
6 6 <NA>    D           NA
7 7 <NA> <NA> -1.195893897
8 8 <NA>    E           NA
9 9 <NA> <NA> -1.269878611

Re: 【R】do.call()関数を用いたリスト（複数のデータフレーム）のmerge
　　投稿者：明石 2019/06/16(Sun) 16:41 No. 22757
青木先生　様；

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

先生がお書きになられたように、簡潔に書きたいというのが意図です。
そこで、do.call()関数に着目しました。

幾つかの方法をお示しくださいましたので、よく理解できました。
reduse()関数を初めて知りましたので、これを機会に勉強させていただきます。

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


低値・高値の基準について
　　投稿者：H 2019/06/14(Fri) 11:49 No. 22750
青木先生
はじめまして。

ホルモン剤内服後のホルモン値が妊娠率に影響するかについて調べております。
ホルモン値低値群・高値群に分け、その間のコントロール群と妊娠率を比較したいのですが、どのようにして低値・高値の基準を設定すればいいでしょうか。
他の研究でもホルモン剤内服後のホルモン値が妊娠率に影響するかについての報告がありますが低値と高値に分け検討しているものがなく検討したいと思っています。

他の研究と内服薬の種類が異なるためか他の研究とホルモンの平均値が異なっていたため、第１四分位数以下、第３四分位数以上で低値・高値を分けコントロール群と比較し有意差がでたのですが恣意的と言われました。

多忙のところ申し訳ありませんが、適切な基準の設定があればご教授お願い致します。
Re: 低値・高値の基準について
　　投稿者：青木繁伸 2019/06/14(Fri) 21:04 No. 22751
> 他の研究でもホルモン剤内服後のホルモン値が妊娠率に影響するかについての報告がありますが

そのような報告ではどのような分析手法が使われていますか？その分析手法を踏襲することに何らかの困難があるのでしょうか？

どのような分析手法が使われているのか知らない状態ではありますが，私がおすすめできる分析手法は「ホルモン値を独立変数，妊娠の有無を従属変数として，ロジスティック回帰分析」をすれば良いのではないですか？
独立変数としてはホルモン値の他の要因も考慮できますし。

データはなるべくもとのまま使うのが原則でしょう。何らかの加工（第１四分位数以下、第３四分位数以上で低値・高値を分けコントロール群と比較し）などを行うと，必ず恣意的ではないか？といわれるのは避けがたいことです。では，といって，別の分類法も使って，別の結果が出たら，どっちが正しいのかといわれるし，いい結果が出たものだけを発表すると，結局は「恣意的だろう」といわれる。
Re: 低値・高値の基準について
　　投稿者：H 2019/06/15(Sat) 01:33 No. 22752
返信ありがとうございます。
他の報告ではホルモン値が100毎に測定してあるものや平均値の低値と高値で評価してあるものがあります。ただそれらでは明確な根拠がないため恣意的になるかと思われます。

投降後にロジスティクス回帰を行いましたが有意差を認めませんでした。
低値群で妊娠率が低いため、統計的に評価したいのですが恣意的になってしまいますか。
宜しくお願い致します
Re: 低値・高値の基準について
　　投稿者：青木繁伸 2019/06/15(Sat) 13:15 No. 22753
> 投稿後にロジスティクス回帰を行いましたが有意差を認めませんでした
ということなら，仕方ないでしょう。

有意ではないが医学的には意味があるようなら，奥の手で，もう少しデータを増やすとか。
Re: 低値・高値の基準について
　　投稿者：H 2019/06/16(Sun) 11:50 No. 22755
返信ありがとうございます。
有意差なしでデータを出そうと思います。
ありがとうございました。


サンプルについて
　　投稿者：林檎 2019/06/12(Wed) 17:25 No. 22744 HomePage
青木先生

はじめまして。
サンプルについてわからないことがあり、質問をさせていただきました。

私は、10施設の企業を対象に質問紙調査による横断研究を実施し、解析をしています。
500名の有効回答があり、施設ごとと従属変数にて一元配置分散分析を行ったところ、1施設のみ合計得点平均が有意に高いことが分かり、指導教授より、その施設は特殊であるため抜いたほうがいいのではないかとご指摘を受けました。ちなみにこの施設の対象は20名でした。
世の中の論文をみると、このような有意差を見ている論文はさまざまであり、今回先生のご意見をお聞かせいただけましたらと思っております。

宜しくお願い致します。
Re: サンプルについて
　　投稿者：青木繁伸 2019/06/12(Wed) 22:03 No. 22745
> その施設は特殊であるため抜いたほうがいいのではないかとご指摘

そもそも，調査の目的は何だったのですか？というか，調査をして，何が分かると思っていたのですか？

「どの施設も同じようである」ということを知りたかったのなら，「1施設のみ合計得点平均が有意に高い」というのは，ビックリするような知見ではないですか。その施設は何らかの問題があると。

対象者が20名ということが特殊なのか，それ以外の条件が特殊なのか。いずれにしても，サンプルサイズが20であっても，それでも有意な差であると言うことですよ。

「その施設は特殊であるため抜いたほうがいい」ということなら，最初からその施設は対象としない方が良かったでしょう。不都合な結果が出たからなかったことにしようというのは，はなはだ科学的ではないです。

その施設を除いて，「どの施設も同じようでした」という結果を出して，何の意味があるのでしょうか？

指導教授とよく話し合うことをお勧めします。
Re: サンプルについて
　　投稿者：林檎 2019/06/13(Thu) 08:02 No. 22747 HomePage
青木先生

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

目的としては、重回帰分析を行い、企業の社員のストレスの要因を明らかにすることです。
重回帰分析に入れるための変数を探すこと、属性を明らかにするために一元配置分散分析を実施したら、有意差があったという結果でした。
調査をお願いするときには、似たような企業だと思い調査を行ったのですが他の施設よりストレスが有意に高かったことが分かりました。

指導教授と相談をしてみます。
ありがとうございました。
Re: サンプルについて
　　投稿者：青木繁伸 2019/06/13(Thu) 10:32 No. 22748
> 他の施設よりストレスが有意に高かったことが分かり

それこそ，その企業が他の企業とどこが違うか（ストレスの原因）を知るために使える重要なサンプルではないですか。
Re: サンプルについて
　　投稿者：林檎 2019/06/13(Thu) 14:00 No. 22749 HomePage
青木先生

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

この企業は大切なサンプルであるということを改めて実感致しましたので、先生からのお返事をもとに解析を進めていきたいと思います。

とても丁寧なご対応をいただき、ありがとうございました。


クラスタリングされた個体の割合の比較とライアンの方法について
　　投稿者：PV 2019/06/09(Sun) 19:11 No. 22739
青木先生、

はじめまして、
ご多忙のところ失礼いたします。

クラスタリングされた個体を群間で比較し、割合が有意に異なるかどうか、またどの群間でことなるか、検定を行いたいのですが、なかなか方法が見つからず調べていたところ、青木先生のホームページにたどり付きました。
以下記載いたしました質問 (1), (2) ご教授いただけましたら幸いです。

下の表は対象生物170個体をクラスター解析し、
各群(A〜F）で各クラスター（V1〜V9)に分配された個体数です。
	V1	V2	V3	V4	V5	V6	V7	V8	V9
A群	6	6	4	3	4	9	3	8	6
B群	2	1	8	4	0	1	1	2	3
C群	9	2	1	1	1	0	0	4	1
D群	8	2	3	2	0	3	4	2	3
E群	1	0	1	0	1	0	0	0	0
F群	3	3	4	9	6	4	4	7	10
質問(1)「クラスターに含まれる個体数の傾向は群間で異なるか」検定を行いたいのですが、
その場合、カイ二乗検定→ライアンの方法で多重比較　を行う方法は正しいでしょうか？

また、青木先生の説明およびＲコードを基に、カイ二乗検定で全体としての差が認められたことを確認後、ライアン法をひとまず行ってみたのですが、
エラーメッセージ「Error: n > 0 は TRUE ではありません 」と表示され、検定を行うことができませんでした。
質問(2)ライアン法は0が含まれていると行えないのでしょうか？

(1), (2) についてご解答のほどよろしくお願いいたします。
Re: クラスタリングされた個体の割合の比較とライアンの方法について
　　投稿者：青木繁伸 2019/06/10(Mon) 09:47 No. 22740
ライアン法は，あなたのデータには使えません（使用例参照）

全ての2群の組み合わせ（A:B, A:C, ..., E:F の15組〕について 2 x 9 の行列で chisq.test を行い，p 値をBonferroni法で解釈する（0.05/15)

手作業でやるのは大変なので，以下のようにすればよい
x = matrix(c(
  6, 6, 4, 3, 4, 9, 3, 8, 6,
  2, 1, 8, 4, 0, 1, 1, 2, 3,
  9, 2, 1, 1, 1, 0, 0, 4, 1,
  8, 2, 3, 2, 0, 3, 4, 2, 3,
  1, 0, 1, 0, 1, 0, 0, 0, 0,
  3, 3, 4, 9, 6, 4, 4, 7, 10),
  nrow = 6, byrow=TRUE)  

g = function(ij) {
  i = ij[1]
  j = ij[2]
  y= rbind(x[i,], x[j,])
  z = colSums(y)
  y2 = y[, z != 0]
  p = chisq.test(y2)$p.value
  cat(i, j, p, "\n")
  return(c(i, j, p))
}

gn = nrow(x)
cn = ncol(x)
p0 = 0.05 / choose(gn, 2)
res = data.frame(t(combn(6, 2, g)))
colnames(res) = c("i", "j", "p.value")
res$sig = ifelse(res$p.value > p0, "n.s.", "*")
print(res)
結果は以下のようになる
   i j     p.value  sig
1  1 2 0.060309971 n.s.
2  1 3 0.095874639 n.s.
3  1 4 0.365146072 n.s.
4  1 5 0.558144611 n.s.
5  1 6 0.384082770 n.s.
6  2 3 0.034933166 n.s.
7  2 4 0.234596686 n.s.
8  2 5 0.246627056 n.s.
9  2 6 0.177712096 n.s.
10 3 4 0.267205617 n.s.
11 3 5 0.438290867 n.s.
12 3 6 0.005751997 n.s.
13 4 5 0.160816574 n.s.
14 4 6 0.077873765 n.s.
15 5 6 0.441102843 n.s.
Re: クラスタリングされた個体の割合の比較とライアンの方法について
　　投稿者：PV 2019/06/10(Mon) 12:05 No. 22741
青木先生、

早速のご解答、ご丁寧な解説を頂き誠にありがとうございました。
Ｒで回して無事回答を得られることができました。

引き続き質問となり大変恐縮ですが、ライアン法について
もう少しお伺いさせていただいてもよろしいでしょうか？

今回私のデータではライアン法は適していないとのことでしたが、
私のデータは、郡ごとに各クラスターにあてはめられた個体数の頻度をみており、
青木先生のライアン法およびチューキー法の説明にある
「いくつかの群の比率に全体として差が見られたときに，どの群間に差があるかを検定する手法。（比率の差の多重比較（対比較））」
にあてはまると（素人的に）思ってしまったのですが
私の解析は「比率の差」を見ているわけではないということでしょうか？

また、私のデータは以下のホームページにある例題（ここから青木先生のＨＰを知りました）

猫も杓子も構造化
http://nekomosyakushimo.hatenablog.com/entry/2017/08/26/222728

のカテゴリー数を増やしただけだと思ったのですが、この解釈は間違っているということでしょうか？
統計解析が素人のため愚問で大変恐縮です。

今後解析を行う上でぜひ参考にさせていただきたいためご解答いただけましたら幸いです。
どうぞよろしくお願いいたします。
Re: クラスタリングされた個体の割合の比較とライアンの方法について
　　投稿者：青木繁伸 2019/06/10(Mon) 12:15 No. 22742
ライアン法，チューキー法は，群（グループ）の数を k としたとき，k x 2 の行列データを対象にするのです。例えば，正否とか失敗・成功のような2値データが対象です。一方（もう一方）の「比率」が対象になります。

一般的に k x m 行列は，群により「分布が違うか」を対象にすることになります。

「カテゴリー数を増やしただけ」ではありません。
Re: クラスタリングされた個体の割合の比較とライアンの方法について
　　投稿者：PV 2019/06/10(Mon) 18:42 No. 22743
青木先生、

なるほど、大変勉強になりました。
お忙しいところ大変丁寧にご説明頂き誠にありがとうございました。


R textでのイタリック表示
　　投稿者：コロン 2019/05/26(Sun) 17:39 No. 22733
お世話になっております。

R plotで描画をしておりますが，グラフ内に p<0.01と表示させたいと思っております。

text(1.5, 70, labels="p <- 0.01)

pだけをイタリックで表示させたいのですが，どのようにすればよろしいでしょうか。

ご教示頂けますでしょうか。
Re: R textでのイタリック表示
　　投稿者：青木繁伸 2019/05/26(Sun) 21:08 No. 22734
実行可能な最小限のテストプログラムを添付してください
Re: R textでのイタリック表示
　　投稿者：コロン 2019/05/27(Mon) 08:38 No. 22735
青木先生

申し訳ございませんでした。簡単なコードではございますが，よろしくお願いいたします。

m1 <- 30
m2 <- 40
x <- c(1,2)
allm <- c(m1, m2)
plot(x, allm, ylim=c(0, 100))
text(1.5, 60, labels="p < 0.05")
Re: R textでのイタリック表示
　　投稿者：青木繁伸 2019/05/27(Mon) 18:13 No. 22736
できばえはちょっといまいちなのですが（字が薄い），以下のようにすれば可能
m1 <- 30
m2 <- 40
x <- c(1,2)
allm <- c(m1, m2)
plot(x, allm, ylim=c(0, 100))
text(1.5, 60, labels="p", vfont = c("sans serif", "italic"))
x2 = 1.58 # この値を locator で決めるか，試行錯誤で決める
text(x2, 60, labels=" < 0.05", vfont = c("sans serif", "plain"))
または，書き始めの位置指定を文字列の左端として（pos=4），2 番目の文字列の最初を p のための空白を一つ確保し，
text(1.5, 60, labels="p", vfont = c("sans serif", "italic"), pos=4)
text(1.5, 60, labels="  < 0.05", vfont = c("sans serif", "plain"), pos=4)
とする（この方が試行錯誤不要なのでよいかも）。
"sans serif" の代わりに "serif" も試して，良いと思う方を使ってください。
Re: R textでのイタリック表示
　　投稿者：コロン 2019/05/27(Mon) 20:52 No. 22737
青木先生

お忙しい中，ありがとうございました。今作成しておりますグラフにも使用することができました。


Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/07(Tue) 18:23 No. 22724
https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
でも紹介されているが，R の
lawstat パッケージの brunner.munzel.test() も
brunnermunzel パッケージの brunnermunzel.test() も
alternative="greater", "less" の場合の信頼区間がおかしい（"two.sided" の場合とまったく同じになっている）
nparcomp パッケージの npar.t.test で method="t.app" を指定して得られるのが正しい信頼区間のようだ。
Re: Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/09(Thu) 11:38 No. 22725
> nparcomp パッケージの npar.t.test で method="t.app" を指定して得られるのが正しい信頼区間のようだ。

と書いたが，そもそも npar.t.test は片側検定で "greater" と "less" を取り違えているようだ。p 値が t.test と比べて逆になっている。これは困ったことになる。
df = data.frame(
  d = c(1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1, 3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4),
  g = rep(1:2, c(14, 11)))
library(nparcomp)
cl = 0.95
alt = "greater"
a = npar.t.test(d ~ g, data=df, round=5, info=FALSE, method="t.app", alternative=alt, conf.level=cl)
a$Analysis
b = t.test(d ~ g, data=df, alternative=alt, conf.level=cl)
b
を実行すると，npar.test の結果は
  Effect Estimator  Lower Upper       T p.Value
1 p(1,2)   0.78896 0.6291     1 3.13747 0.00289
だが，t.test の結果は
t = -2.9486, df = 15.916, p-value = 0.9953
p 値が 0.00289 と 0.9953 と 180 度異なっている。
alt="less" では 0.99711 と 0.004739 である。
ちなみに，wilcox.test の p 値は t.test と同一方向だ。
Re: Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/10(Fri) 18:00 No. 22727
引き続き

並べ替え検定については，奥村先生が
https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
でも独自のプログラム例を示しているが，
brunnermunzel パッケージには brunnermunzel.permutation.test がある。
しかし，これは場合によって正しい答えを出さない。というか，奥村先生のプログラムと結果が異なることがある。
そのようなデータの例は，x = c(36, 45, 51, 49, 57)，y = c(48, 42, 49, 39, 45) の場合。

奥村先生のプログラムでは 0.365079365079365 になるのだが，brunnermunzel.permutation.test では 0.341269841269841 になる。

前者は 92/252，後者は 86/252 だ。92 の内訳は，データを並べ替えたデータにより計算される検定統計量の絶対値が，観察値における検定統計量の絶対値と同じかより極端な場合の数 45*2 と，検定統計量が +Inf, -Inf になる 2 通りの計 92 である。

観察値における検定統計量の絶対値と同じになるのは 4*2 = 8 通りあるが brunnermunzel.permutation.test では，計算誤差のためか 8 個全てが同じとなっていないのであろう。つまり，数え落としている。

また，Inf, -Inf となる場合の扱いもチェックする必要があるだろう。
前に書いて消したが，検定統計量が Inf になるのは R の場合だけで，Python の場合には 0 による割り算エラーで計算が中止される。

いろいろ問題ありで，それぞれのプログラムどれも信頼して使うわけにはいかな状態である。

自分で書くしかないか。
Re: Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/10(Fri) 21:15 No. 22728
Brunner-Munzel検定 最強！！

なんて，書かれているページが多くて，ですね。ちょっと，立ち止まって検証してみませんか？という趣旨で書きました。

ほかの人の意見も聞きたいので，「拡散希望」（^_^;）
Re: Brunner-Munzel検定
　　投稿者：it may be that 2019/05/19(Sun) 19:00 No. 22729
npar.t.test の greater/less の扱いですが、取り違えているわけではないと思いますよ。
summary(a)
に
#---------------------------Interpretation---------------------------------------------#
 p(a,b) > 1/2 : b tends to be larger than a
と出ます。
したがって、群1と2 (g1とg2) の観測量をX1, X2として、対立仮説は
'greater'はP(X1 < X2) > 1/2, 
'less'はP(X1 < X2) < 1/2
です。
Estimate が 0.78896 > 1/2 であることから、P(1,2)はP(X2 - X1 > 0)であると思われます。
よって、書き直すと
'greater'はP(X2 - X1 > 0) > 1/2,
'less'はP(X2 - X1 > 0) < 1/2 あるいは P(X2 - X1 < 0) > 1/2
となります。

一方、g1とg2の位置母数をそれぞれm1とm2とすると、t.testの対立仮説は
'greater' :m1 - m2 > 0
'less' : m1 - m2 < 0
です。
つまり、npar.t.testの対立仮説は
m2 - m1 > 0
m2 - m1 < 0
に相当し、基準が逆になっています。
確かに紛らわしいですが、lm, glmなどでcontr.treatmentに馴染んでしまった身としては、むしろt.testやwilcox.testのfactorの扱いに違和感を覚えます。
普段contr.SASを使っていればt.test,wilcox.testが自然でしょう。
Re: Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/21(Tue) 11:01 No. 22730
「npar.t.test の greater/less の扱い」についてのコメント，ありがとうございます。
確かに，説明の全体を見ればその通りですね。
Re: Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/21(Tue) 11:04 No. 22731
追加の考察

http://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

では，p = Pr(x＜y)+0.5*Pr(x=y) ではなく，t 分布による漸近近似を使っている。

n1, n2 が小さいときに行う並べ替え検定で，n1, n2 が大きいことを仮定する漸近近似を行うという，相矛盾する状況である。

p そのままを使って並べ替え検定を行うと，
foo2 = function(X) {
  brunner.munzel.test(xandy[X], xandy[-X])$estimate
}
z2 = combn(1:N, n1, foo2)
bm2 = brunner.munzel.test(x, y)$estimate
mean(abs(z2-0.5) >= abs(bm2-0.5))
0.006690223 が得られる。

なお，この値の正当性としては，マン・ホイットニーの U 検定統計量と Brunner-Munzel の p との間に p=1-U/n1/n2 という単純な関係があるので，マン・ホイットニーの U 検定統計量を計算して使う並べ替え検定と同じ値になることを示せば十分。

または，coin パッケージの wilcox_test で正確な p 値を求めると，同じ値 0.006690223 が得られる。結局は wilcox_test が並べ替え検定をやっているということだ。
> df = data.frame(d=c(x, y), g=factor(rep(1:2, c(length(x), length(y)))))
> wilcox_test(d ~ g, data=df, distribution="exact")

    Exact Wilcoxon-Mann-Whitney Test

data:  d by g (1, 2)
Z = -2.6934, p-value = 0.00669
alternative hypothesis: true mu is not equal to 0
ということは，マン・ホイットニーの U 検定と Brunner-Munzel 検定って，前者が正規分布，後者が t 分布（など）で漸近近似を行うだけで，本質は同じということだ...
> npar.t.test(d ~ g, data=df, info=FALSE, method="perm", nperm=400000, round=5)$Analysis
       Estimator Statistic   Lower   Upper p.value
id       0.78896   3.13747 0.59708 0.99057 0.00708
logit    0.78896   2.38394 0.57650 0.91561 0.00694
probit   0.78896   2.51949 0.58109 0.92446 0.00694
となれば，奥村先生のページの冒頭の説明，

> 二つの確率変数 X1，X2 が同じ分布に従うという帰無仮説を検定するには，有名なWMW検定（Wilcoxon-Mann-Whitney test）が使えます。しかし，分布が異なる（例えば分散が異なる）場合には，これでは正確に検定できません。

> そこで，分布が同じことは仮定せず，両群から一つずつ値を取り出したとき，どちらが大きい確率も等しいという帰無仮説を検定することにします。

という，Brunner-Munzel 検定の優位性というのもなくなるのか。（前提が異なるといえども，検定統計量が同じなんだからねぇ）
Re: Brunner-Munzel検定
　　投稿者：青木繁伸 2019/05/24(Fri) 18:57 No. 22732
検定統計量ではなく，その近似（t分布に近似する）のときに「ベーレンス・フィッシャー問題を考慮した」ということですね。


クロンバックαの係数
　　投稿者：wtcs 2019/05/05(Sun) 08:24 No. 22721
青木先生

はじめまして。お忙しい中 またお休み中に失礼いたします。

昨年尺度開発を行い、現在はその開発尺度を用いて縦断調査にて影響要因を探索しています。
尺度開発の論文を執筆し投稿しましたところ、査読者から信頼性の検討のところでご指摘をいただきました。
尺度は全15項目で、下位因子は5項目、4項目、6項目の3因子です。
クロンバックαが0.89,0.77,0.70で全体では0.78でした。
それぞれの下位因子のクロンバックαが全体のクロンバックαより高いのはおかしいとご指摘をいただいたのですが、どこを調べてもわかりませんでした。文献も同じようなデータのものもあります。
確かに項目数が多い全体の方が高くなりやすいとは思いますが、必ずしも全体の方が高い値でないとその尺度の信頼性は低いのでしょうか？
尺度開発は初めてのため、初歩的な質問だと思いますがどうぞよろしくお願いいたします。
Re: クロンバックαの係数
　　投稿者：青木繁伸 2019/05/05(Sun) 15:00 No. 22722
ある数学的（統計学的）命題が誤りであることを示すには，反例を一つあげれば十分ですね。

妥当なデータを，信頼のおけるプログラムで分析した結果がその命題に反するものであった場合，その命題は誤っているということでしょう。

今の場合，命題は「下位尺度のクロンバックのアルファは全体のクロンバックのα信頼性係数より小さい」ということ。

信頼のおけるプログラムは，たとえば R とか SPSS でいいでしょう。
妥当なデータというのはちょっとやっかいで，確かにそのようなデータがあるということでは不十分なのです。その統計手法が想定している条件を満たしているデータでなければいけません。

クロンバックのα信頼性係数は各質問項目への回答の不偏分散の和と各質問項目への回答の合計値（尺度値）の不偏分散で決まります。合計値の計算の元になる各質問項目間の相関の符号は一致していますか？もし，符号が一致しない項目があると合計を取るときに項目間で多少の値の打ち消しが起こり合計値の分散は小さくなり，必然的に各項目の分散の合計との割合は大きくなります。そして，クロンバックのαは小さくなります。もし，下位尺度にこのような符号の不一致がない場合はその尺度のクロンバックのα信頼性係数は小さくはならないでしょう。そうすれば，全体のクロンバックのα信頼性係数のほうが小さくなってしまうということもありえるでしょう。

以下のようなn=15, 9変数（下位因子数3個ずつ）のデータで
x1  x2  x3  y1  y2  y3  z1  z2  z3
39  37  50  72  50  52  47  63  65
42  48  39  46  47  56  48  48  40
48  54  50  53  43  55  52  55  58
32  30  29  43  26  48  25  50  40
54  39  39  46  57  55  47  31  48
42  50  44  34  38  49  59  29  55
56  62  61  52  50  54  60  49  75
57  51  55  61  56  39  68  51  51
47  65  69  58  59  51  61  61  53
48  50  45  36  51  24  43  43  41
73  66  54  40  56  49  46  60  44
49  48  47  45  43  39  45  54  45
46  50  56  58  67  50  55  49  45
61  56  51  54  59  60  48  57  50
57  44  61  53  47  67  45  51  40
全体のクロンバックのα信頼性係数は 0.8022528 であり，下位尺度のクロンバックのα信頼性係数は 0.8212634, 0.5153157, 0.4307981 となっていることが，エクセルでちょっとした計算をするとわかります。
原因は，相関係数行列をみるとわかりますが，4つの項目間の相関係数がほんのちいさな負の相関を持っているからです。
通常は，尺度を作る際にこのような項目は選ばれないと思いますが，こんなこともあるよということを示すために掲示しておきます。
Re: クロンバックαの係数
　　投稿者：wtcs 2019/05/05(Sun) 17:21 No. 22723
青木先生

早々のご回答をありがとうございました。
データの因子行列を確認しました。符号の一致しない（負の相関）因子構造はありませんでした。
分散の大きさなどはあまり意識していませんでした。とても勉強になりました。
青木先生のデータ例から、査読者へも回答できそうですし、尺度の信頼性を確認でき安心しました。

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


勉強中
　　投稿者：T-test 2019/05/01(Wed) 15:52 No. 22716
青木先生

お世話になります。お忙しい中大変申し訳ないのですが、アドバイスを頂ければと思います。

現在、同じ人物が回答したpooled データの解析を行っており、pooled T-test ができないか方法を探しております。

例：

血液型 健康良い群　　　　悪い群
A　　　　　　30　　　　　　　　45　　　
B　　　　　　29　　　　　　　　39　　
O　　　　　　34　　　　　　　　50
AB　　　　　35　　　　　　　　　42

30などの数字は、医療費です。

この場合、各血液毎、例えばA型の列の30と45で、対応のないT-testは、統計としてダメでしょうか？

サンプルは、同じ人物が3か月ごとに回答した結果が3回分あり、その合計です。
Pooled データのため、どのように有意差を解析していいかわかりません。
お手数をお掛けし大変申し訳ございませんが、どうぞよろしくお願い致します。
Re: 勉強中
　　投稿者：青木繁伸 2019/05/02(Thu) 10:57 No. 22717
実験でサイン的にはあまり芳しいものではないと思うが，繰り返しのない二元配置分散分析でしょうかね？
http://aoki2.si.gunma-u.ac.jp/lecture/TwoWayANOVA/TwoWay.html（繰り返しが1，すなわち繰り返しがない場合）
Re: 勉強中
　　投稿者：T-test 2019/05/02(Thu) 19:17 No. 22718
青木先生

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

私の書き方が悪かったのですが、同じ回答者が1回から3回、繰り返して回答しています。

データは、
ID 血液　健康状態　医療費
1　　A 良い　　 30
1 A　　　 良い　 　 28　　
1 A　　　　良い　　　29
2 B　　　　悪い　　　35　
2 B　　　　悪い　　　40
2 B　　　　良い　　　32
3 O　　　　良い　　　35
3 O　　　　良い　　　32　

ただし、3回回答している者と2回や1回の者がおります。

このようなデータでも、繰り返しのない二元配置でよろしいでしょうか？

多重比較は、シェッフェでよろしいでしょうか？

分散分析は初めてで、勉強してみたいと思います。
お手数をお掛けし大変申し訳ございませんが、どうぞよろしくお願い致します。
Re: 勉強中
　　投稿者：青木繁伸 2019/05/03(Fri) 14:08 No. 22719
> このようなデータでも、繰り返しのない二元配置でよろしいでしょうか？

よくわからないけど。無理だと思う。

複数の観察がある場合には（ランダムに）どれか一つに絞れば，対応のない二元配置分散分析にはなりますよね。

分析のことを考えずにデータを取ってしまうと，後で困るという典型例です。
実験デザイン（調査デザイン）をしっかり立ててからデータを取るべきでした。
Re: 勉強中
　　投稿者：T-test 2019/05/03(Fri) 21:25 No. 22720
青木先生

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

はい、その通りで、データの二次利用のため、非常に困っております。
なんでこういう調査を行うのかと思うのですが、やってしまった後ではどうにもなりません。

対応のない二元配置分散分析、勉強してみます。

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


lda関数の判別係数の再利用について
　　投稿者：初学者 2019/04/14(Sun) 20:16 No. 22713
青木先生

こんにちは。
いつも掲示板で勉強させていただいております。

このたびは正準判別分析の判別係数を別の研究で再利用できないかと思い、
ご質問させていただきました。
Rのlda関数で分析しております。
print()で出力される判別係数を使って、別の新しいデータの分類に利用できれば、
と考えているのですが、そのようなことは可能なのでしょうか。
R上ではpredict()等で簡単に実現できると思うのですが、
Rが使いこなせない人のために、例えばEXCEL上で判別係数から新データの分類が
簡便にできないかと考えました。

ご教示のほど、どうぞよろしくお願いいたします。
Re: lda関数の判別係数の再利用について
　　投稿者：青木繁伸 2019/04/15(Mon) 11:58 No. 22714
library(MASS)
MASS:::predict.lda
の 2 行を入力すると 89 行のソースプログラムが表示されます。
まるまる全部が必須というわけではないですが，その程度の処理が必要ということです。
それを Excel のワークシート機能だけを使って R を知らない人のために用意するというのはあまりお勧めできません。たとえ用意できたとして，使い方もある程度は複雑になるでしょう。データを用意してスタートボタンを押すだけというようにすることもできるでしょうが，かなりたいへんでしょう。それくらいなら predict(obj, newdata) と入力するのが遙かに簡単では？
Re: lda関数の判別係数の再利用について
　　投稿者：初学者 2019/04/15(Mon) 12:28 No. 22715
青木先生

早速のご回答、誠にありがとうございました。
たしかに、おっしゃる通りだと感じました。

他の人に使っていただくときには、
・判別分析から学習したobjectをRDSファイル等にセーブ
・RDSファイル、実行するRのscript、簡単な解説、などを渡す
・predict(obj,newdata)してもらう
等の方法で検討してみたいと思いました。

大変勉強になりました。
ありがとうございました。


マンホイットニーのＵ検定〜同順位の補正〜
　　投稿者：確率は難しい 2019/04/03(Wed) 14:29 No. 22708
青木先生　どうかよろしくお願いします。

マンホイットニーのＵ検定について、先生の解説
http://aoki2.si.gunma-u.ac.jp/lecture/Average/U-test.html
を見ると、同順位の場合に補正する数式が掲載されています。
この数式通りやってみたのですが、同順位が全体の中で多いような分布の場合に確率がうまく出ない場合があることに気がつきました。

例えば、満月の日の翌日は雨が降りやすいという仮説に対して、帰無仮説を立てるとして、
１６日間観測して、満月の翌日だけに１００ミリの雨が降り、あとは雨が降らなかったとします。
第１群 満月の翌日ではない日→雨量０のデータが１５個
第２群 満月の翌日→雨量１００のデータが１個
これを上記の同順位の数式も含めて計算してみると、確率が0.0001ということになりました。

しかし、満月の日の翌日は平均気温が低くなるということで考えると、小数点以下３桁ぐらいまで平均気温を計算して同順位がないとすると、
第１群 気温９.２５３、１０．８６５度・・・・など、９度〜１１度の日が１５個
第２群 気温８．０２５度のデータが１個（最も低い）
これで計算してみると、確率は0.1ということになりました。

感覚的には、満月の日だけに雨が降ったとしても、確率が0.0001というのはおかしいと思いまして、後者の確率（０．１）の方が感覚に合っているように思います。

この、同順位の式を数学的に理解することは私には難しすぎると思いますが、それはさておき、

これを苦肉の策で解決するために、データに適当に小さい乱数（他のデータと逆転しないような十分に小さい乱数）を足して無理に順位をつけるような数字に置き換えて計算しました。例えば、０の日の雨量を、 0.143、0.054、0.658・・・・というふうに置き換えた計算です。そうすると、同順位がなくなり、気温の例と同じ確率になります。

このように、乱数をプラスして同順位をなくして解決する方が、同点の場合の数式を使うよりも実態に合っているように思うのですが、致命的な問題はありますでしょうか。
また、他に良い解決策はないでしょうか？
よろしくお願いします。
Re: マンホイットニーのＵ検定〜同順位の補正〜
　　投稿者：青木繁伸 2019/04/03(Wed) 18:03 No. 22709
stats の wilcox.test で，同順位があるために「タイがあるため、正確な p 値を計算することができません」となる場合は，coin パッケージの wilcox_test を使ってみましょう。
library(coin)
x <- c(rep(0, 15), 100)
g <- factor(c(rep("a", 15), "b"))
wilcox_test(x ~ g, distribution="exact")
以下が表示されます。
	Exact Wilcoxon-Mann-Whitney Test

data:  x by g (a, b)
Z = -3.873, p-value = 0.0625
alternative hypothesis: true mu is not equal to 0
Z は，同順位補正をしたときの正規近似の Z 値で，表示されている p 値の算出には使われていません。

同時に表示されている p-value = 0.0625 は正確な p 値です。

coin パッケージの出現前には exactRankTests パッケージの wilcox.exact が使われていました。
library(exactRankTests)
wilcox.exact(x ~ g)
以下が表示されます。
	Exact Wilcoxon rank sum test

data:  x by g
W = 0, p-value = 0.0625
alternative hypothesis: true mu is not equal to 0
W は wilcoxon 統計量です（どちらを1群にするかで，もう一方の極端な側の W = 15 が返されることもあります）。表示された p 値は，coin:::wilcox_test と同じ 0.0625 です。

もう少し他の道具でも試してみるには，
http://aoki2.si.gunma-u.ac.jp/R/index.html
の中の，
XV. exact 検定とモンテカルロ法による近似検定
	3. マン・ホイットニーの U 検定(exact test)
に記載した R プログラム exact.mw を使ってみると以下のようになります。
exact.mw(x=rep(0, 15), y=1)
U = 0, Z = 3.87298, P 値 = 0.000107511
正確な P 値 = 0.0625
査察した分割表の個数は 2 個
示されている U = 0, Z = 3.87298 は，coin:::wilcox_test の Z と同じですね。

次の行に正確な p 値が示されています。0.0625 で，coin:::wilcox_test，exactRankTests:::wilcox.exact の p 値と同じです。

exact.mw は exact.mw(matrix(c(15, 0, 0, 1), 2)) のようにも書けますが，
http://aoki2.si.gunma-u.ac.jp/exact/utest/getpar.html
にオンライン版を用意してあります。
列数を 2 として，次の画面で
15 0
0 1
の2×2分割表を入力すれば， 0.0625 が表示されます。

ということで，正確な p 値は 0.0625 でよいようです。

====================

最後に，「乱数をプラスして同順位をなくして解決する方が、同点の場合の数式を使うよりも実態に合っているように思う」ということですが，乱数を加えて検定を行うという流儀は存在します。
しかし，今回の場合は
stats:::wilcox.test(1:15, 16)

	Wilcoxon rank sum test

data:  1:15 and 16
W = 0, p-value = 0.125
alternative hypothesis: true location shift is not equal to 0
とすることと同じなので，わざわざ乱数を加えて同順位をなくすということがどの程度適切か疑問ですね。
Re: マンホイットニーのＵ検定〜同順位の補正〜
　　投稿者：確率は難しい 2019/04/04(Thu) 10:16 No. 22710
青木先生

アドバイスありがとうございました。
結局のところ、ｚ値を用いて正規分布を用いるのは近似に過ぎないということに尽きるということだと理解しましたが、そういう理解で合ってますでしょうか？

実際に適用して統計分析をしており、同順位が少なく、データ数が多い場合には比較的正確な値に近く、有効に適用できていると思っています。たくさんのケースをエクセルで分析して、機械的にＺ値を算出し、正規分布だと仮定してｐ値を出しているため、
ひとつひとつ統計ソフトに入力するというのは、出来れば避けたいと考えています。

この正確なｐ値の計算は、エクセルの表計算で単純に行えるようなものではないのでしょうか？
Re: マンホイットニーのＵ検定〜同順位の補正〜
　　投稿者：青木繁伸 2019/04/04(Thu) 12:19 No. 22711
それぞれのデータを wilcox_test(x ~ g, distribution="exact") で分析する数行のプログラムを書けばよいだけでしょう。

エクセルの表計算のほうが手間が掛かるのでは？
Re: マンホイットニーのＵ検定〜同順位の補正〜
　　投稿者：確率は難しい 2019/04/04(Thu) 15:28 No. 22712
先生　ありがとうございます。

タイムリーな回答をいただきまして、非常に助かりました。
本当にありがとうございました。

おそらく、「同順位が多くて、近似的なＰ値と正確なＰ値が大きく乖離するケース」はそれほど多くないと思いますので、その部分だけプログラムを書くように勉強してみます。


二項ロジスティック回帰分析
　　投稿者：春 2019/03/30(Sat) 02:00 No. 22701
青木先生

初めて質問をさせていただきます。
spssを用いて二項ロジスティック分析を行っていますが、どうにもわからなくて困っております。ご教示いただけないでしょうか。よろしくお願いいたします。

仮説に基づいて、独立変数群を5モデル作り、それぞれの影響をみるため順番に強制投入しています。以下の２つの方法を試みましたが、投入したモデルは同じなのに結果が違っています（B,標準偏差、wald,有意差、EXP（B)）。なぜ違っているのか、また、どちらのやり方が適しているのかを教えていただけないでしょうか。

�‘販�変数投入の際にブロックごとに次々と指定していった場合
�▲屮蹈奪�ごとの指定せずに「総当たり法」を繰り返し、�，里箸�と同じブロックごとに変数を増やした場合

以上です。

お忙しいところ、申し訳ありませんが、どうか、よろしくお願いいたします。
Re: 二項ロジスティック回帰分析
　　投稿者：青木繁伸 2019/03/30(Sat) 09:16 No. 22702
（2）の実際の分析過程があいまいなのではっきりしたことが言えませんが，一般的には，「分析に用いる変数セットにより，結果的に選択される変数セットが異なることはよくあること」

（1）の「独立変数投入の際にブロックごとに次々と指定」する場合は，投入順序により結果が異なることがある。増減法であろうと同じことが起きる。
増加法と減少法，増減法と減増法で結果が違うのもよくあること。

「ブロックごとの指定せずに「総当たり法」を繰り返し、（1）のときと同じブロックごとに変数を増やした場合」も同じ。ブロックごとに総当たり法で増やすべき変数を決めても，結局は（1）と同じように，投入順序で結果が違うこともある。

「分析に用いる変数セット」に依存しない方法は，全ての変数を対象に総当たり法でモデルを選ぶ。もっとも，理論的に解釈可能なモデルであるかはチェックする必要がある。

=======

以下のような思考実験をするとよく分かる。

以上のようにして「最良のモデル」が得られた後，ある独立変数を分析に加えたらどうなるだろうか？

その変数が橋にも棒にもかからないヘボ変数で，結果は変わらない。
その変数が，新たにモデルに加えられる。
その変数が，モデルに加えられる一方，今までにモデルに含まれていた変数が不適切とわかり除去される（変数1個とは限らない）
その変数が，モデルに加えられる一方，今までにモデルに含まれていた変数が不適切とわかり除去される（1個とは限らない），そして今までモデルに含まれていなかった変数がモデルに加えられる（1個とは限らない）。
その際，同じブロック内の変数に入れ替えが生じる。
また，ブロック内に適当な変数がなくなる。
さらに，ブロック内の複数の変数がモデルに含まれる。
などなど，全ての状況が考えられる。

だって，それが多変量解析だもの
Re: 二項ロジスティック回帰分析
　　投稿者：春 2019/03/30(Sat) 12:09 No. 22703
青木先生

早々に明快なご説明をいただきまして、誠に、ありがとうございます。

�，皚△睚竸瑤療蠧�順序は同じで、結果的にも変数は除去されなかったのですが、値（B,標準偏差、wald,有意差、EXP（B)）が異なるのは、なぜだろうかと不思議に思っていました。
しかし、先生が明示してくださった「全ての変数を対象に総当たり法でモデルを選ぶ」方法で解析を進めようと思います。ありがとうございました。

今後も、どうぞよろしくお願い致します。
Re: 二項ロジスティック回帰分析
　　投稿者：青木繁伸 2019/03/30(Sat) 16:58 No. 22704
文意がとりにくかったのですが，「最終的に採用された変数は同じ」ということですか？

データに欠損値があったりしますか？
2つの方法で，実際に使用されたサンプルサイズがいくつかわかりますか？
変数が同じでも，実際に分析に使われるサンプルサイズが異なれば，B,標準偏差、wald,有意差、EXP（B)）が異なるのは当たり前です。
Re: 二項ロジスティック回帰分析
　　投稿者：春 2019/03/31(Sun) 23:27 No. 22705
青木先生

ご連絡ありがとうございます。
本当に、このように親身になって頂いてありがたいです。

ご質問して頂いてように「最終的に採用された変数は同じ」です。
」
そして、「実際に使用されたサンプルサイズ」」を確認すべく、もう一度、spssの結果のlog�，鉢△鮓�比べてみました。すると、青木先生の仰るようにサンプルサイズが違っていました！
�，諒�法だと投入した全ての変数の欠損値を含むので、�△離汽鵐廛襯汽ぅ困茲蠅眈�さくなっていました。
これが原因だったのですね。よくわかりました。

ちなみに、�，任呂匹離皀妊襪�Nの値が同じになりますが、�△世肇皀妊襪瓦箸�Nの値が違ってくるということになりますね。論文にモデル1、モデル2、モデル3、モデル4、モデル5と階層的に示したいときには、Nを揃えた�，諒�法が適切なのでしょうか？
何度も質問させていただき、恐縮ですが、ご教示いただけないでしょうか。

どうぞ、よろしくお願い致します。
Re: 二項ロジスティック回帰分析
　　投稿者：青木繁伸 2019/04/01(Mon) 08:29 No. 22706
モデルの評価を AIC などで行う場合は，サンプルサイズが同じでないと評価不能ですね。

ただし，対象者の属性により欠損値が生じやすくなるような状況がある場合は，結果の評価に注意が必要になるかもしれません。

なお，変数選択の結果を見て，モデルに残った変数だけを指定して再度分析を行うべきです。
以下のような状況を考察すれば分かるでしょう。
変数が a ~ e，対象者が 9 人。o は測定値がある（欠損値ではない），x は欠損値。
       a b c d e
case1  o o o o x
case2  o o o x x
case3  o o x o o
case4  o o o o o
case5  o o o o o
case6  o o o o o
case7  o o o o o
case8  o o o o o
case9  o o o o o
5変数全部を使って分析を始める場合，c, d, e が欠損値の case1, case2, case3 は分析対象から外されます。分析対象は 6 人です。
その後の変数選択で a, b, d が残ったとします（c, e は採用されなかった）。
これで終わってはいけません。
a, b, d を使って分析することにすると，欠損値を持つとして分析から外されるのは case2 のみとなります。分析対象は 8 人です。 case1 も case3 も分析対象となり，2 人増えます。
Re: 二項ロジスティック回帰分析
　　投稿者：春 2019/04/01(Mon) 10:56 No. 22707
青木先生

何度も丁寧にご教示いただきまして誠にありがとうございました。

ようやく分かってまいりました。
Nを揃えるときの注意点に気をつけながら、モデルに残った変数でやり直しをしてみます。

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


Split plot + RCBD with repeated measures
　　投稿者：森 2019/03/25(Mon) 18:32 No. 22700
青木先生

自分で調べてみましたが、解決できなかったので質問させていただきます。
本当に申し訳ありませんが、よろしくお願いいたします。

10つの森林を対象に、それぞれの森林に種Aの樹木1個体と樹種Bの樹木1個体を植栽し、それらの個体の高さ（樹高）を2ヶ月おきに観察しています。
私が明らかにしたいのは、樹高に対する種の影響と経時的な変化（測定時期の影響）、それらの交互作用です。
森林をブロックとして捉え、その中に、処理区（樹種）を設置し、さらに、その経時変化を調べていますが、それぞれのブロックの処理区に個体が1つずつしか存在しません。
経時変化の測定はそれぞれの個体を繰り返し測定します（repeated measures）。
これを統計ソフトRで解析したいのですが、一般化線形混合モデルを用いる場合、変量効果をどのようにモデルに組み込めばいいのかが分かりません。
glmer(樹高~樹種＋(1|森林/樹種/測定時期)＋測定時期＋樹種×測定時期)
などの形を考えていますが、これでいいのかも分かりません。
よろしければご教示いただけないでしょうか。

また、2つ目の質問になってしまいますが（不躾で本当に申し訳ありません）、上記の実験デザインの中で、それぞれのブロックの処理区に2つの個体を配置した場合はどのようなモデルになりますでしょうか。
こちらにつきましても、
glmer(樹高~樹種＋(1|森林/樹種/個体/測定時期)＋測定時期＋樹種×測定時期)
などの形を考えていますが、よくわかっていません。
よろしければお教えいただけないでしょうか。

私の勉強不足なところをお尋ねしてしまい、大変申し訳ありません。
お手数をおかけいたしますが、どうぞよろしくお願いいたします。


最小選択肢にほとんど丸のつかない可能性のある尺度開発
　　投稿者：風 2019/03/11(Mon) 08:13 No. 22695
青木先生
初めて質問させていただきます。
初学者ですよろしくお願い致します。

現在、自己評価尺度の開発を行っています。尺度原案の信頼性妥当性検証を行っている
ところです。
尺度原案のリッカートを6件法「1．全く当てはまらない」、「2．当てはまらない」、「3．どちらかというと当てはまらない」、「4．どちらかというと当てはまる」、「5．当てはまる」、「6．非常に当てはまる」として調査を行いました。
返送されてきた調査票は、「4.どちらかというと当てはまる」「5．当てはまる」に回答が集まる傾向がありました。

このように最小選択肢にほとんど〇のつかない可能性のある尺度はどうなのでしょうか…？

（リッカートを変更する予定ですが…）
Re: 最小選択肢にほとんど丸のつかない可能性のある尺度開発
　　投稿者：青木繁伸 2019/03/11(Mon) 09:20 No. 22696
もっと極端な場合，つまり一つの選択肢にほとんど全てが回答するという場合は，そのような質問は不適切ということになりますね。
そのような場合と比べてどうすべきかということでしょう。明確な基準というものはないと思います。他の研究者とのコンセンサスが得られるかどうかも判断基準になるでしょう。
Re: 最小選択肢にほとんど丸のつかない可能性のある尺度開発
　　投稿者：風 2019/03/11(Mon) 12:20 No. 22697
青木先生

早速お返事をいただきありがとうございました。
極端な場合は、質問項目として除外されることになると考えます。

他の研究者の意見も伺って検討していきたいと思います。ありがとうございました。


重回帰分析の基準変数
　　投稿者：北海道 2019/03/08(Fri) 08:59 No. 22692
久しぶりにメールさせていただきます。よろしくお願いします。
重回帰分析の基準変数と言っていいのかわからないのですが、「不安を持ちやすい」という変数をリッカートで大変当てはまる、当てはまる、どちらかというと当てはまる、どちらかというと当てはまらない、当てはまらない、全く当てはまらないの6つで聞きました。最近は6件だとダミー変数を作らない方法でも良いようなので、他の変数と一緒にそのまま投入しspssで解析しました。全く当てはまらないを基準としてそれ以外のベーターが出てきたので、これを表にしました。
論文を投稿したところ、査読者が全く当てはまらないのベーターがない、といってきました。査読者にどのように説明したらいいかご教授願えますでしょうか
Re: 重回帰分析の基準変数
　　投稿者：青木繁伸 2019/03/08(Fri) 09:35 No. 22693
「全くあてはまらない」のβは 0 なのです。
それが不都合だと思う人はその人のかって（統計知らずなだけ）なのですが。
あなた自身も，「カテゴリー変数を重回帰分析に使う」ということ，その際に「ある変数を基準とすること」について十分理解されていないようです（あなたはダミー変数を使ったつもりはなくても，解析プログラムはダミー変数を使っています）。それらについて書いてある本とか WEB ページで勉強すれば，あなたの言葉で説明できるようになるでしょう。

あまりよくないページが多いのですが，
http://jojoshin.hatenablog.com/entry/2016/05/08/005747
などは，結果の解釈法についても簡潔に書かれていると思います。
Re: 重回帰分析の基準変数
　　投稿者：北海道 2019/03/08(Fri) 09:42 No. 22694
青木先生
早速のお返事ありがとうございました。
よく分かりました。
まだまだ未熟ですが、これからも研究を続けていきます。


Python
　　投稿者：波音 2019/03/01(Fri) 17:31 No. 22691
質問でも意見でもなく恐縮ですが、、

何気に今気づいたのですが、青木先生はPythonのスクリプトを大量にまとめられていたのですね。

http://aoki2.si.gunma-u.ac.jp/Python/index.html

少し勉強しようと意気込みましたが、Rでできるからといって学ぶのを呆けておりました。
（先生には感服いたします）


簡単なサンプルサイズの求め方について
　　投稿者：清水せつ子 2019/02/26(Tue) 17:59 No. 22686
ある特定の母集団の血液中の化学物質の濃度を求めるために実験を委託したいと考えており，
血液中の化学物質の濃度の90％タイル値を信頼水準95％で予測するためのサンプル数を計算しています．
計算はできる限りシンプルにしたいと思っています．

私は，過去のデータを整理して，血液中の化学物質の濃度を正規分布に従い，母分散が既知と仮定することができるから，母平均の検定に使われる検出力の式で求められるのではないかと考えています．

検出力＝P（｜z｜＞Zalpha/2）

ここで，Zは，標本平均を期待値と分散で標準化した値．alphaは有意水準．

そんな中，某研究仲間から，二項分布に従うと仮定すれば，一つのサンプルをとり，それが90％タイル以上となる確率は0.9，n個のサンプルのすべてが90％タイル以上となる確率は0.9＾nだから，

1-0.9＾n= (1-信頼水準)　　　　(1)

という数式でサンプルサイズnが求まると言われました．

・その研究者に「血液中濃度を二項分布に仮定することはできないのではないか」と問いましたところ，不良品の計算と同じだから，0.9％タイル値を超えるか超えないかといった二値でみればいいのだと言われました．また，血液濃度の分布は，ノンパラメットリックで見ているという回答を得ました．

母集団の分布が変われば，90％タイル値を95％の信頼水準で予測するためのサンプルサイズは異なるはずですが，(1)式では母集団がどのような分布であっても，サンプルサイズは同じとなるのでおかしいと思います．

すなわち，この回答に様々な理論的な間違いが含まれていると思いますが，どこに間違いがあるか今のところ答えが見つかりません．

・さらに，二項分布で仮定して，90％タイル値を超えるか超えないかといった二値でみれるならば，
(1）式は　1-p(n)=(1-信頼水準），ここでp(n)はBin(n,p)ではないかと思います．

以上，不適切な質問かもしれませんが，以上の2点についてご助言いただけませんでしょうか．

なお，私がこのFAQで問いたいのは，上記の間違いを理論的に説明したく，そのご助言をいただきたいということです．参考のため，私が以前に作成した血液濃度のサンプルサイズに関するファイルを添付しようと思いましたが，ファイルサイズが大きく添付できませんでした．この件でまた悩みがでましたら，改めて投稿させていただきます．
Re: 簡単なサンプルサイズの求め方について
　　投稿者：青木繁伸 2019/02/26(Tue) 18:50 No. 22687
利用上の注意に書いてあるように，

画像などのバイナリーファイルをアップロードすることが可能です。

添付可能ファイル : PNG, ZIP
最大投稿データ量 : 150 KB

ということです。


> それが90％タイル以上となる確率は0.9
　　？ 0.1 では？

> 0.9％タイル値を超えるか超えないか
　　？ 90%タイルの書き間違い？

検出力と信頼率（信頼度）を混同している？
Re: 簡単なサンプルサイズの求め方について
　　投稿者：清水せつ子 2019/02/27(Wed) 11:17 No. 22688
青木先生

早速の回答ありがとうございます．
添付ファイルについては承知いたしました．
ご丁寧な情報をありがとうございます．

> それが90％タイル以上となる確率は0.9
　　？ 0.1 では？

> 0.9％タイル値を超えるか超えないか
　　？ 90%タイルの書き間違い？

はい，1-0.9，90％タイルが正しいです．

また，次の文章にも訂正があります．

>n個のサンプルのすべてが90％タイル以上となる確率は0.9＾nだから，

0.9＾n→(1-0.9＾n)です．

どれも私の書き間違いです．
申し訳ありません．

・(1)式は検出力と信頼度を混合しているということですね．
そのとおりだと思いました．
まとめます．

検出力というのは，ある検定方法を用いて,ある値thetaに対する帰無仮説を棄却する確率です．
つまり，標本平均が母平均と同等という帰無仮説を立てた場合，標本平均が母平均と異なる対立仮説が示せる確率，ある標本の不良率が，母不良率と同等であるという帰無仮説を立てた場合は，それが異なる対立仮説が示せる確率を表すものです．

一方，信頼度は，帰無仮説が正しいとされる確率です．ここでは，95％信頼水準が当てはまります．

検定を行うことによる間違いが2種類あって，帰無仮説が正しいのに間違いだとする過誤を示すのが信頼度，対立仮説が正しいのにそれを間違いだとする確率を示すのが（1-検出力）です．

(1）式では，明らかに，検出力と信頼度が区別されていません．

・また，信頼度か検出力を求めようとしている式が，1-0.9＾nというのも間違い．

もし混合していて，(1）式の左辺が検出力，右辺が1-信頼水準としているのであれば，

左辺は不良率と同じように検出力を求めることができ，1-0.9＾nではなく，

1-P(n)，

ここで，p(n)は二項分布(n,p),nは標本数，pは不良確率

となります．

以上，これが私の理解です．
間違いがあればご指摘ください．よろしくお願いいたします．
Re: 簡単なサンプルサイズの求め方について
　　投稿者：青木繁伸 2019/02/27(Wed) 12:10 No. 22689
> 帰無仮説が正しいのに間違いだとする過誤を示すのが信頼度

違います。「帰無仮説が正しいのに間違いだとする過誤」は第一種の過誤（その過誤を犯す確率がα）。信頼度は 1-α で表される。同じようにしてデータを収集し，信頼限界を計算したとき，その信頼限界が母数（帰無仮説で仮定した統計量）を含む確率。信頼度をを大きくするにはサンプルサイズを大きくする必要がある。

> 検出力というのは，ある検定方法を用いて,ある値thetaに対する帰無仮説を棄却する確率です．

微妙だけど正確ではない。検出力とは，帰無仮説が間違っている場合に，帰無仮説を棄却できる確率（1-β）。βは，帰無仮説が間違っているのに帰無仮説を棄却できない確率。サンプルサイズを大きくすれば，検出力は大きくなる。

つまり，あなたの場合は，「血液中の化学物質の濃度の90％タイル値」を「誤差範囲内」で推定する際に，推定値「90%タイル値±誤差」が真値である「90%タイル値」を信頼度95%で含むようにするためにはサンプルサイズが幾つ必要かということでしょう。しかも，調査は1回しか行われないので，その調査が成功する為には検出力も高くする必要があるので，サンプルサイズはもっと必要であると。

母平均を求めるのならば http://aoki2.si.gunma-u.ac.jp/lecture/SampleSize/muconf.html ，特にその後半の方法が使えるのだが，90%タイル値を求めるのはこれではだめでしょうね。

解を求める直接的な方法はないかもしれない。モンテカルロ法で求めることはできるかも。
Re: 簡単なサンプルサイズの求め方について
　　投稿者：清水せつ子 2019/02/27(Wed) 12:33 No. 22690
青木先生

早速のお返事ありがとうございます．

訂正ありがとうございます．
また，モンテカルロ法に対する助言についてもありがとうございます．
大変参考になります．
平均ではなく，90％タイル値の求め方については，少しやってみて，後日，躓きましたら再度投稿いたします．

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


独立変数の水準サイズについて
　　投稿者：畠山 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乗検定を遺伝子型の各組み合わせでやって、組み合わせ数で割ることにします。

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

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

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