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

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


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

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


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

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

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

ご助言いただけますでしょうか。
よろしくお願いいたします。

Re: 介入研究の因子分析
  投稿者: 2019/06/17(Mon) 15:24 No. 22758
青木先生

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

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

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


【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
青木先生
はじめまして。

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

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

多忙のところ申し訳ありませんが、適切な基準の設定があればご教授お願い致します。

Re: 低値・高値の基準について
  投稿者:青木繁伸 2019/06/14(Fri) 21:04 No. 22751
> 他の研究でもホルモン剤内服後のホルモン値が妊娠率に影響するかについての報告がありますが

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

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

データはなるべくもとのまま使うのが原則でしょう。何らかの加工(第1四分位数以下、第3四分位数以上で低値・高値を分けコントロール群と比較し)などを行うと,必ず恣意的ではないか?といわれるのは避けがたいことです。では,といって,別の分類法も使って,別の結果が出たら,どっちが正しいのかといわれるし,いい結果が出たものだけを発表すると,結局は「恣意的だろう」といわれる。

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)「クラスターに含まれる個体数の傾向は群間で異なるか」検定を行いたいのですが、
その場合、カイ二乗検定→ライアンの方法で多重比較 を行う方法は正しいでしょうか?

また、青木先生の説明およびRコードを基に、カイ二乗検定で全体としての差が認められたことを確認後、ライアン法をひとまず行ってみたのですが、
エラーメッセージ「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
青木先生、

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

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

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

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

猫も杓子も構造化
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)してもらう
等の方法で検討してみたいと思いました。

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


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

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

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

しかし、満月の日の翌日は平均気温が低くなるということで考えると、小数点以下3桁ぐらいまで平均気温を計算して同順位がないとすると、
第1群 気温9.253、10.865度・・・・など、9度〜11度の日が15個
第2群 気温8.025度のデータが1個(最も低い)
これで計算してみると、確率は0.1ということになりました。

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

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

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

このように、乱数をプラスして同順位をなくして解決する方が、同点の場合の数式を使うよりも実態に合っているように思うのですが、致命的な問題はありますでしょうか。
また、他に良い解決策はないでしょうか?
よろしくお願いします。

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:青木繁伸 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: マンホイットニーのU検定〜同順位の補正〜
  投稿者:確率は難しい 2019/04/04(Thu) 10:16 No. 22710
青木先生

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

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

この正確なp値の計算は、エクセルの表計算で単純に行えるようなものではないのでしょうか?

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:青木繁伸 2019/04/04(Thu) 12:19 No. 22711
それぞれのデータを wilcox_test(x ~ g, distribution="exact") で分析する数行のプログラムを書けばよいだけでしょう。

エクセルの表計算のほうが手間が掛かるのでは?

Re: マンホイットニーのU検定〜同順位の補正〜
  投稿者:確率は難しい 2019/04/04(Thu) 15:28 No. 22712
先生 ありがとうございます。

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

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


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

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

仮説に基づいて、独立変数群を5モデル作り、それぞれの影響をみるため順番に強制投入しています。以下の2つの方法を試みましたが、投入したモデルは同じなのに結果が違っています(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乗検定を遺伝子型の各組み合わせでやって、組み合わせ数で割ることにします。

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


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

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

 図示した場合、信頼区間の破線が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

> ##############x軸の値を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)
>
>
> ##############x軸の値を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
オフセット項を含む単回帰モデルの切片の標準誤差の推定方法について,ご教示いただければ幸いです.
以下のようなRのスクリプトを実行したとき,切片の推定値と標準誤差を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の活用など,ご指摘ありがとうございます.
丁寧にソースコードを追って,どのような計算が行われているのかを確認したいと思います.
ご助言いただけましたこと,また,このような質問の場を維持していただいていることに感謝申し上げます.


【R】 文字列の部分的マッチング 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

この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: 【R】 文字列の部分的マッチング 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: 【R】 文字列の部分的マッチング 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本のベクトル)をリストで得る
Rプログラムをご教示いただければ大変に勉強になります。

新年早々、お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

Re: 【R】 文字列の部分的マッチング 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: 【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/04(Fri) 08:19 No. 22666
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
新年早々に、大変に良い勉強をさせていただきました。
この例題を通じて、今まで不得手であった sapply、lapplyも理解できました。
心から御礼を申し上げます。

Re: 【R】 文字列の部分的マッチング 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: 【R】 文字列の部分的マッチング 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: 【R】 文字列の部分的マッチング charmatch
  投稿者:明石 2019/01/04(Fri) 18:58 No. 22669
青木先生 様;

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

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

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

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


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

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

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

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

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

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

・評価基準(価格、品質、外見)
・代替案(商品A、商品B、商品C、商品D)
・一対比較表
 添付の画像ファイルにお示しをします。

# 各評価基準の重要度

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: 【R】青木先生のAHPプログラムの使用方法
  投稿者:青木繁伸 2018/12/21(Fri) 10:54 No. 22660
それでよいと思いますが

Re: 【R】青木先生の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)
で、行と列のサイズの与え方に不安がありましたので、ご確認をさせていただきました。

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

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

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

意思決定分野のAHPとDEA(包絡分析法)について、利用できるRプログラムを調べております。

AHPについては、CRANのRパッケージがありますが、データ作成方法が煩雑で、苦慮していました。

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

DEAについても、CRANのRパッケージがありますが、マニュアルの理解で苦慮しています。

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

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

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


確証的因子分析のモデル適合について
  投稿者:さとう 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が出る仕組み自体がまずいのでに、それを正の値しか出ない方法に変えたいのです。

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

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

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