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

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


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

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


重回帰での切片有り無しモデルについて
  投稿者:さくらとお城 2016/09/29(Thu) 17:24 No. 22153
いつもお世話になっております。

下記のデータで切片有無の重回帰を行った結果の解釈に悩んでいます。
R G B N
1 99.0 128.0 27.0 1.61
2 91.5 99.5 40.5 2.16
3 90.0 113.0 30.0 2.32
4 72.5 72.0 45.0 3.07
5 80.0 77.5 46.0 2.97
6 110.0 142.0 24.0 1.58
7 96.5 114.5 31.5 1.87
8 98.0 119.5 30.0 2.02
9 112.5 139.5 21.5 1.61
10 88.0 104.0 36.5 2.02
11 81.5 74.0 44.0 2.87
12 90.0 102.5 34.5 1.99
13 84.5 92.5 39.0 2.17
14 99.5 120.5 34.5 2.49
15 81.0 79.0 44.0 2.99
16 91.0 103.0 38.0 2.25
17 95.5 111.5 32.5 2.07
18 85.5 88.0 42.5 2.31
19 83.0 79.5 45.5 2.59
20 81.0 89.0 38.5 2.36

data=read.table("data.txt",header=TRUE)

lm1=lm(N~R+G+B+0,data=data)
lm2=lm(N~R+G+B,data=data)

summary(lm1)
AIC(lm1)

summary(lm2)
AIC(lm2)

以下はRでの結果の一部だけです。

> summary(lm1)
lm(formula = N ~ R + G + B + 0, data = data)

Residual standard error: 0.2448 on 17 degrees of freedom
Multiple R-squared: 0.9904, Adjusted R-squared: 0.9888
F-statistic: 587.5 on 3 and 17 DF, p-value: < 2.2e-16

> AIC(lm1)
[1] 5.21581

> summary(lm2)
lm(formula = N ~ R + G + B, data = data)

Residual standard error: 0.2407 on 16 degrees of freedom
Multiple R-squared: 0.765, Adjusted R-squared: 0.7209
F-statistic: 17.36 on 3 and 16 DF, p-value: 2.764e-05

> AIC(lm2)
[1] 5.327225

切片無しモデルはAdjusted R-squared: 0.9888で,良すぎるようにも見えて,大丈夫でしょうか?

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

Re: 重回帰での切片有り無しモデルについて
  投稿者:青木繁伸 2016/09/30(Fri) 07:06 No. 22154
基本的には,「理論的に考えて」独立変数が 0 のときに従属変数は 0 になるだろうかということに尽きるでしょう。
データだけから見る場合,R, G を独立変数とした単回帰では原点を通らないので,重回帰でも R, G, B が 0 のとき原点を通るというモデルはどうなんでしょうか?
切片なしのモデルと切片ありのモデルでは決定係数の計算方法が異なるので,両者を比較することはできません。

Re: 重回帰での切片有り無しモデルについて
  投稿者:さくらとお城 2016/09/30(Fri) 08:29 No. 22155
青木先生,早速のご返事ありがとうございます。

>基本的には,「理論的に考えて」独立変数が 0 のときに従属変数は 0 になるだろうかということに尽きるでしょう。

今回のデータは,RGBでNを予測する目的で実験で得たものです。理論的にRGBが0なら,Nも0になることについて,確証もないし,確かめることも難しい。

>切片なしのモデルと切片ありのモデルでは決定係数の計算方法が異なるので,両者を比較することはできません。

これに関連した図書や参考論文があれば,ぜひ教えてください。

Re: 重回帰での切片有り無しモデルについて
  投稿者:青木繁伸 2016/09/30(Fri) 17:00 No. 22156
以前はこの件のページがあったのですが,Excel の馬鹿馬鹿しさ,および Excel をひいき・寵愛・擁護するひとからのクレームによりそのページ(と Excel に関する全ページ)は削除しました。

当たり障りのないところを再現すると

======= 引用開始
スネデッカー・コクラン著(畑村,奥野,津村訳)「統計的方法 原書第6版」岩波書店の6章18節「原点をとおる直線のあてはめ」というところに詳しく書いてあります。

SPSS はどんな風にしているかなと思ってやってみて,以下のようなことであるらしいと思います。
  x  y  x^2  x・y  予測値     y^2  予測値^2
1 1 1 1 0.40779221 1 0.16629448
2 2 4 4 0.81558442 4 0.66517794
3 2 9 6 1.22337662 4 1.49665036
4 4 16 16 1.63116883 16 2.66071176
5 3 25 15 2.03896104 9 4.15736212
6 3 36 18 2.44675325 9 5.98660145
7 2 49 14 2.85454545 4 8.14842975
8 2 64 16 3.26233766 4 10.6428470
9 3 81 27 3.67012987 9 13.4698533
10 4 100 40 4.07792208 16 16.6294485
合計 385 157 76 64.0233766

a = 0.4077922 = 157 / 385
R^2 = 0.8424129 = 64.0233766 / 76 =予測値の平方和 / yの平方和
R = 0.9178305

分散分析表
  平方和  自由度 平均平方    F値
回帰 64.023377   1   64.0233766 48.1112557
残差 11.976623   9  1.33073593
全体 76     10
回帰の平方和は「予測値の平方和」,全体の平方和は「yの平方和」

プログラムの出力に下手でさっぱりわけの分からない日本語訳の注が付いています(SPSS Release 6.1.2)。
曰く,

Note # 10572
原点(切片なしのモデル)の場合,R二乗は,回帰によって証明(ママ;「説明」のミスタイプだと思う)される原点について,yにおける変動性の割合を測定します。これは,切片を含むモデルのR二乗値とは,比較できません。

========== 引用終了

まあ,定義が違うのだから比較できるわけはないのですけどね。

R の summary.lm がやっていることを抜粋すると以下のようになります。
> func = function(z) { # R.squared を求める箇所
+ r = z$residuals
+ f = z$fitted.values
+ if (attr(z$terms, "intercept")) { # 切片のあり/なし判定
+ # 切片ありのモデル 予測値と予測値の平均値からの偏差の二乗和
+ mss = sum((f - mean(f))^2)
+ } else {
+ # 切片なしのモデル 予測値の二乗和(予測値と0からの偏差の二乗和)
+ mss = sum(f^2)
+ }
+ rss = sum(r^2) # 残差の二乗和
+ r.squared = mss/(mss + rss) # 決定係数
+ c(r.squared = r.squared, mss = mss, rss = rss)
+ }
>
> lm1 = lm(Sepal.Length ~ ., iris) # 切片あり
> summary(lm1)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.17127 0.27979 7.760 1.43e-12
Sepal.Width 0.49589 0.08607 5.761 4.87e-08
Petal.Length 0.82924 0.06853 12.101 < 2e-16
Petal.Width -0.31516 0.15120 -2.084 0.03889
Speciesversicolor -0.72356 0.24017 -3.013 0.00306
Speciesvirginica -1.02350 0.33373 -3.067 0.00258

Residual standard error: 0.3068 on 144 degrees of freedom
Multiple R-squared: 0.8673, Adjusted R-squared: 0.8627
F-statistic: 188.3 on 5 and 144 DF, p-value: < 2.2e-16

> func(lm1)
r.squared mss rss
0.8673123 88.6118483 13.5564851

> lm0 = lm(Sepal.Length ~ . + 0, iris) # 切片なし
> summary(lm0)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
Sepal.Width 0.49589 0.08607 5.761 4.87e-08
Petal.Length 0.82924 0.06853 12.101 < 2e-16
Petal.Width -0.31516 0.15120 -2.084 0.03889
Speciessetosa 2.17127 0.27979 7.760 1.43e-12
Speciesversicolor 1.44770 0.28149 5.143 8.68e-07
Speciesvirginica 1.14777 0.35356 3.246 0.00145

Residual standard error: 0.3068 on 144 degrees of freedom
Multiple R-squared: 0.9974, Adjusted R-squared: 0.9973
F-statistic: 9224 on 6 and 144 DF, p-value: < 2.2e-16

> func(lm0)
r.squared mss rss
0.9974049 5210.2935149 13.5564851



Rを使った成績などの判定
  投稿者:初学者 2016/09/26(Mon) 20:19 No. 22150
お忙しいところ失礼致します。
Rについて、ご教示いただきたいと思いました。
データを集計していたのですが、例えば図のように数値と成績が対応している場合、各対象の数値から成績判定するには、どのようにコードを組めばよいのでしょうか。
初歩的なことで申し訳ございません。ご教示いただけますと幸いです。

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


Re: Rを使った成績などの判定
  投稿者:青木繁伸 2016/09/26(Mon) 21:24 No. 22151
やり方はいろいろありますが,区切りが不等間隔なこともあるので,一般的な方法を示しましょう。findInterval 関数を使います。
x <- seq(6.0, 10.5, by=0.1)
y <- 11-findInterval(x, c(0, 6.7, 6.9, 7.1, 7.3, 7.6, 8.0, 8.5, 9.1, 9.8))
data.frame(x, y)
findInerval の第2引数 c(0, 6.7, 6.9, 7.1, 7.3, 7.6, 8.0, 8.5, 9.1, 9.8)
0 以上 6.7未満を 1,6.7 以上 6.9 未満を 2,...,9.8 以上を 10 にするための区切りです。
11 からその値を引けば,欲しい数値になりますね。

Re: Rを使った成績などの判定
  投稿者:初学者 2016/09/26(Mon) 22:45 No. 22152
早速のお返事をいただき、誠にありがとうございます。
ご教示いただいた通り、findinterval関数を使って一度やってみます。

ありがとうございました。
今後とも何卒よろしくお願いいたします。


R 分割表から元データを再現
  投稿者:明石 2016/09/22(Thu) 19:56 No. 22147
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。

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

−−−

青木先生の以下のサイトについて、ご質問がございます。

 http://aoki2.si.gunma-u.ac.jp/R/tenkai.html
 分割表から元データを再現 

使用例として、以下のように記載されています。

f <- matrix(c(2,3,1,4,3,5,4,3,2,1,2,4), nrow=3, byrow=TRUE)

result <- tenkai(f) # 結果はリストで返される
x <- result$x
y <- result$y

x
[1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3

y
[1] 1 1 2 2 2 3 4 4 4 4 1 1 1 2 2 2 2 2 3 3 3 3 4 4 4 1 1 2 3 3 4 4 4 4

この例題に、行名と列名を、以下のように与えたとします。
rownames(f) <- c("r1", "r2", "r3")
colnames(f) <- c("c1", "c2", "c3", "c4")

上記のx、yの値を、数字ではなく、行名と列名で返したいと思います。

そのRプログラムをご教示いただければ、大変に助かります。

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

Re: R 分割表から元データを再現
  投稿者:青木繁伸 2016/09/22(Thu) 21:15 No. 22148
せっかく与えた行名と列名を使うなら,
rownames(f)[x]
colnames(f)[y]
ただ,tenkai が何をやっているのかを見て,単に r, c で始まる名前でよいならば,行名,列名を与えてやる必要はなくて,
paste("r", rep(row(f), f), sep="")
paste("c", rep(col(f), f), sep="")
で十分。これを関数にしてやると面倒くささがなくなる。
paste0 が使えるならば,以下のように。
paste0("r", rep(row(f), f))
paste0("c", rep(col(f), f))

【御礼】 Re: R 分割表から元データを再現
  投稿者:明石 2016/09/23(Fri) 08:45 No. 22149
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

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

私も自分で、以下のような方法を考えてみました(途中で断念)。
・x,y からデータフレームを作成し、
・x,y が返す数字を、因子の水準番号と見なして、
 因子の水準に対応する名前(ここでが、行と列の名前)を与えて、
・最後に、数値を文字列表示に変換する

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


R ベクトル要素の完全一致検索
  投稿者:明石 2016/09/19(Mon) 15:17 No. 22144
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。

初歩的な質問で大変に恐縮ですが、何卒どうぞ、よろしくお願いいたします。

−−−

Rに組み込みの文字列ベクトル state.name(アメリカ50州の名前)があります。

state.name

"Alabama" "Alaska" "Arizona" "Arkansas" "California" "Colorado" "Connecticut" "Delaware" "Florida" "Georgia" "Hawaii" "Idaho" "Illinois" "Indiana" "Iowa" "Kansas" "Kentucky" "Louisiana" "Maine" "Maryland" "Massachusetts" "Michigan" "Minnesota" "Mississippi" "Missouri" "Montana" "Nebraska" "Nevada" "New Hampshire" "New Jersey" "New Mexico" "New York" "North Carolina" "North Dakota" "Ohio" "Oklahoma" "Oregon" "Pennsylvania" "Rhode Island" "South Carolina" "South Dakota" "Tennessee" "Texas" "Utah" "Vermont" "Virginia" "Washington" "West Virginia" "Wisconsin" "Wyoming"

state.name の50要素から、
文字列の長さが最大、最小となる要素(州の名前)を
州の名前、要素の番号を対にして、すべて取り出すことを考えます。

nagasa <- nchar(state.name)

文字列の長さが最大となる要素(州の名前)は、複数あり、
> state.name[ nagasa == max(nagasa) ]
[1] "North Carolina" "South Carolina"

文字列の長さが最小となる要素(州の名前)は、複数あり、
> state.name[ nagasa == min(nagasa) ]
[1] "Iowa" "Ohio" "Utah"

のように、州の名前はすべて取り出すことはできますが、
その要素の番号までは、とりだすことができません。

初歩的な質問で大変に恐縮ですが、何卒どうぞ、よろしくお願いいたします。

Re: R ベクトル要素の完全一致検索
  投稿者:青木繁伸 2016/09/19(Mon) 22:50 No. 22145
which を使えば良いでしょう

> which(nagasa==max(nagasa))
[1] 33 40
> state.name[which(nagasa==max(nagasa))]
[1] "North Carolina" "South Carolina"
> which(nagasa==min(nagasa))
[1] 15 35 44
> state.name[which(nagasa==min(nagasa))]
[1] "Iowa" "Ohio" "Utah"
> which(nagasa==12) # 任意の長さの州名を取り出す
[1] 34 38 39 41
> state.name[which(nagasa==12)]
[1] "North Dakota" "Pennsylvania" "Rhode Island" "South Dakota"

蛇足ではあるが,which.max, which.min では仕様にあわないが,それはそれで使い道がある
> which.max(nagasa)
[1] 33
> state.name[which.max(nagasa)]
[1] "North Carolina"
> which.min(nagasa)
[1] 15
> state.name[which.min(nagasa)]
[1] "Iowa"

【御礼】 Re: R ベクトル要素の完全一致検索
  投稿者:明石 2016/09/20(Tue) 08:32 No. 22146
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

which.max()関数、which.min()関数を使うことを考えましたが、
最初に検出した1つに要素のみを返すことから、
重複の場合も考慮した、which.max()関数、which.min()関数のような関数を
探しましたが、見つけることができませんでした。

青木先生からご教示いただきました内容で、
which.max()関数、which.min()関数の活用の幅が広がりました。

よく分かりました。

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


母集団数の調整
  投稿者:蒲郡 2016/09/11(Sun) 19:55 No. 22139
お世話になります。

いくつかのイベントの発生率を地域間で比較するときに、母集団数の差を調整して計算しようと思います。

経験ベイズ法で計算した結果、すべての地域の発生率が同値になったイベントがありました。

この場合は「母集団数の差を調整する必要がない」と判断して、このイベントのみを非調整値として、他のイベント発生率とともに論じることは妥当なのでしょうか。
あるいは、すべてのイベントを非調整値で論じることが妥当なのでしょうか。

そもそも、先生方は、母集団数の差を調整する/しないは、どのように決めていますか。
ご教示いただけますと幸いです。

Re: 母集団数の調整
  投稿者:青木繁伸 2016/09/12(Mon) 17:05 No. 22140
母集団数とは,「母集団のサイズ」つまり「母集団の構成要素の個数」ですか?(そんな用語は統計学にはないと思いますが)。
例えば,47都道府県を47個の母集団とし,母集団のサイズとは各都道府県の人口みたいなことですか?

このような場合は,標本調査ではなく全数調査なので,母数(パラメータ)が観察されるので,検定とか推定は無意味でしょう?
つまり,ある県の例えば「○○の悪性新生物による死亡率」が人口10万対 12.345 ならそれは母比率ですということ。別の県で人口10万対 9.876 だとすれば,差は -2.469 で,それ以上でもそれ以下でもない。人口で調整する意味もない(性・年齢別人口で調整する意味はある--年齢調整死亡率)。

Re: 母集団数の調整
  投稿者:蒲郡 2016/09/12(Mon) 19:51 No. 22141
青木先生侍史

お忙しい中、ご教示をいただきありがとうございます。
用語を誤っており、お手数をおかけして失礼しました。
先生が推測して下さった通りの内容をお尋ねしたかった次第です。

特定の年齢別の人口で調整した場合に、年齢調整死亡率がすべての都道府県で同じ値になった場合は、「調整した結果、差はなかった」と考えて良いのでしょうか?

私の最初の質問を読みなおしてみますと「このイベントのみを非調整値として、他のイベント発生率とともに論じる」のは真実と異なっていると考えました。

Re: 母集団数の調整
  投稿者:青木繁伸 2016/09/13(Tue) 00:48 No. 22142
> 特定の年齢別の人口で調整した場合に、年齢調整死亡率がすべての都道府県で同じ値になった場合は、「調整した結果、差はなかった」と考えて良いのでしょうか?

例では都道府県としましたが,47 の年齢調整死亡率がすべて同じになった?そんなことあるのかなあ?

まあ,同じ値になったのなら「全く差はない」ということでいいでしょう。

> 「母集団数の差を調整する必要がない」

調整したから差がないことがわかったのでしょう?だったら,調整する必要はあるのでしょう。違いますか?

Re: 母集団数の調整
  投稿者:蒲郡 2016/09/14(Wed) 01:16 No. 22143
青木先生侍史

ご教示をいただきありがとうございます。
同じになることがあるのか??という疑問点については、ソフトウェアをRに替えて計算し直したところ、わずかですが差がありました。
ありがとうございました。とても勉強になりました。


比の検定
  投稿者:瑯崎 2016/09/07(Wed) 20:28 No. 22136
皆様平素よりお世話になります。下記のことを検定したい場合、どのような検定手法になるかご教示頂けないでしょうか。
ある同じ品種の植物を5個体育てて、地上部(茎や葉等)と地下部(根)の質量を測定した。地上部と地下部の質量比が一定であるを帰無仮説にして危険率5%で検定したいのです。
個体No 地上部質量 地下部質量
1 20.3 5.6
2 30.7 7.8
3 16.1 3.1
4 16.8 6.9
5 25 7.3

よろしくお願いします

Re: 比の検定
  投稿者:鈴木康弘 2016/09/09(Fri) 07:44 No. 22137
 思うに、方法のひとつは、質量比の値5つについて、その平均の信頼区間を示すこと。

http://aoki2.si.gunma-u.ac.jp/lecture/Average/Mean2.html

 もう一つの方法は、地上部と地下部の相関係数を出して、それが0でないことの検定。

http://aoki2.si.gunma-u.ac.jp/lecture/Corr/corr.html

 かな?

Re: 比の検定
  投稿者:瑯崎 2016/09/09(Fri) 08:37 No. 22138
ありがとうございます


R 文字列の連結
  投稿者:明石 2016/08/28(Sun) 21:02 No. 22125
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。
何卒どうぞ、よろしくお願いいたします。

−−−

文字列ベクトルがあります。
簡単な例示です。
"1", "12", "12345", "1234", "123"

各要素の文字列の長さが不揃いですが、
最大の長さ(ここでは、要素"12345"の長さ5)にあわせて、
必要に応じて、先頭に"0"を連結して、
すべての要素の文字列の長さをあわせることを考えます。

出来上がりは、以下のようになります。

"00001", "00012", "12345", "01234", "00123"

ベクトルの長さが膨大(数十万件)なことから、効率的な方法を検討しています。

以下の2点について、ご教示をいただければ、大変に助かります。
(1) 先頭に"0"を連結する方法
(2) applyを使う方法

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

Re: R 文字列の連結
  投稿者:青木繁伸 2016/08/29(Mon) 07:53 No. 22126
やりかたはいろいろあるけど,一番簡単(?)なのは sprintf を使うもの。
a = c("1", "12", "12345", "1234", "123")
n = max(nchar(a))
unname(sapply(a, function(b) sprintf("%0*s", n, b)))
* には n の値がつかわれ,実際には sprintf("%05s", b) と同じになる。
別の方法は
unname(sapply(a, function(b) paste(paste(rep("0", n-nchar(b)), collapse=""), b, sep="")))
のようなものがあるけど,見た目でもわかるが,sprintf を使うより実行時間は若干長い。
一番外側の unname は,文字列要素に名前が付くのをはぎ取るためなので,場合によっては不要。

【御礼】 Re: R 文字列の連結
  投稿者:明石 2016/08/29(Mon) 08:51 No. 22127
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

今回も助けていただきました。
心から御礼を申し上げます。

例えば、"00000"という文字列を作成するには、
paste( rep( "0", 5) , sep="", collapse="" )
ということまでは自分でも分かりましたので、
個々の要素の文字列の長さから、先頭につける文字列の長さの計算はできました。

ただし、数十万件をループで回している方法は何とも回避したいと思っておりました。

大変に良い勉強をさせていただきました。
ありがとうございました。

Re: R 文字列の連結
  投稿者:青木繁伸 2016/08/29(Mon) 11:27 No. 22128
sapply する必要はなかったですね。これだと unname も不要です。
a = c("1", "12", "12345", "1234", "123")
n = max(nchar(a))
sprintf("%0*s", n, a)

Re: R 文字列の連結
  投稿者: 2016/08/29(Mon) 14:57 No. 22129
青木先生

横から失礼します。
R-3.2.5(Linux)では"0"の代わりに半角スペースで埋められてしまいました。

?sprintfで調べると
## Platform-dependent bad example from qdapTools 1.0.0:
## may pad with spaces or zeroes.
sprintf("%09s", month.name)
のような例が載っていました。

unname(sapply(a, function(b) paste(paste(rep("0", n-nchar(b)), collapse=""), b, sep="")))
では0で埋まりました。

参考までに。

Special Thanks(Re: R 文字列の連結)
  投稿者:明石 2016/08/29(Mon) 20:07 No. 22130
青木先生、荒様

いつもお世話になり、ありがとうございます、明石と申します。

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

2つ、ご報告がございます。
・私は、Windows7、R-3.3.1の環境ですが、荒様と同じ結果となりました。
・100万件のデータに適用しましたが、爆走でした。

勇気を振り絞って、投稿させていただいて本当に良かった、と思いました。
ありがとうございました。

Re: R 文字列の連結
  投稿者:青木繁伸 2016/08/29(Mon) 21:34 No. 22131
Platform dependent ですか。

もとの文字列に空白が含まれないならば,
gsub(" ", "0", sprintf("%0*s", n, a))
の方が速いかもしれません。

Re: R 文字列の連結
  投稿者:明石 2016/08/30(Tue) 19:21 No. 22135
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

お示しをしてくださいました
gsub(" ", "0", sprintf("%0*s", n, a))

恐ろしいほど、早いです。

Rの文字列の連結は、paste()関数を使用することしか知りませんでしたので、
今回、sprintf()関数を教えていただきましたことを契機に、調べ物をしています。

大きな武器を得ました。
心から御礼を申し上げます。
ありがとうございました。


対応のある2要因の分散分析における交互作用についてのpost-hoc test
  投稿者:ポスドクR 2016/08/29(Mon) 22:28 No. 22132
青木先生、皆様
お世話になります。

現在、あるリハビリ介入に関する論文を投稿しているのですが、標記のことに関して査読者からの指摘に関して皆様のお知恵をお借りいたしたく投稿いたします。

今回の研究では、健常被験者20名に、Aという介入とBという介入を行い、その前後でのパラメータの変化(即時効果)を調べています。それぞれの被験者は、AまたはBの介入をランダムな順番で、十分な間隔(数週間以上)開けて行っています。
得られたデータとしては、各被験者でのA介入前、A介入後、B介入前、B介入後の4つです。
私たちは、介入の前後 と 介入の種類 の2要因での対応のある分散分析と考えて解析を起こったところ、それぞれの主効果は有意差なく、交互作用のみが有意な結果となりました。
我々はpost-hoc testとして介入前後、介入の種類のそれぞれの要因について、単純主効果を検討し、記載しました。(それぞれ前/後、A/Bしかないので多重比較は行わず、各条件におけるpaired t-testを行いました。)
すると、査読者から、「2要因の分散分析で有意な交互作用が得られた場合、すべてのパターンについての一元配置分散分析(A前、A後、B前、B後の4つについての分散分析)を行い、多重比較補正を用いて有意な組み合わせを検出することが必要である。」
との指摘を受けました。

いろいろな文献や資料を参考に調べてみましたが、査読者の指摘している内容が正しいというう確証が得られず、さりとて、それが誤っているとのはっきりした証拠も得られず、大変困っています。

post-hoc testに関する査読者の主張の妥当性について、ご意見をいただけましたら幸いに存じます。
よろしくおお願い致します。

Re: 対応のある2要因の分散分析における交互作用についてのpost-hoc test
  投稿者:青木繁伸 2016/08/29(Mon) 23:27 No. 22133
「2要因の分散分析で有意な交互作用が得られた場合」という長いキーワードでネット検索するといくつか出てきますが,そのような場合には「単純主効果の検定を行う」らしく,SPSS での分析法について説明してあるページもあるようですね。

「すべてのパターンについての一元配置分散分析(A前、A後、B前、B後の4つについての分散分析)」というのは,サンプルサイズの水増し(4倍,あなたのデータだとサンプルサイズが80!!)になるし対応のあるデータであるということについても考慮していないので,疑問だと思います。

Re: 対応のある2要因の分散分析における交互作用についてのpost-hoc test
  投稿者:ポスドクR 2016/08/30(Tue) 13:05 No. 22134
青木先生

早速のコメントありがとうございます。
やはり対応のあるデータで無理やり一元配置に持っていくのはちょっと問題ということですね。
心強いお返事ありがとうございました。

何とか差読者のご機嫌を損ねないようにやんわり返事を書いてみます。

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


R merge()関数でのエラーメッセージへの対処方法
  投稿者:明石 2016/08/24(Wed) 22:04 No. 22122
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。
何卒どうぞ、よろしくお願いいたします。

−−−

merge()関数で、現在の私の知識、経験では理解できない事象に遭遇し、
途方にくれています。

2つの csvファイルを読み込み、データフレームに代入します。
変数名は、dat1、dat2 とします。

dat1 を左側、dat2 を右側として、左外部結合でmergeします。

結合キー"注文番号"のデータ型は、integer です。

dat3 <- merge(dat1, dat2, by.x="注文番号", by.y="注文番号", all.x = TRUE, all.y = FALSE)

これら2つの csvファイルを、EXCELのxlsxファイルとして保存し、
xlsx形式で読み込むと、なんと、上記のマージ処理で、エラーが発生します。

Error in fix.by(by.x, x) : 'by' must specify a uniquely valid column

キー"注文番号"はユニークであることは、unique()関数で確認済みです。
重複はございません。

csvファイルならばよいが、xlsxファイルならばエラーになるという事象に途方にくれています。

EXCELのxlsxファイルの読み込みは、library(openxlsx) を使いましたが、
他のEXCELファイル読み込みのパッケージ(XLConnect、readxl)でも同様でした。

このエラーメッセージへの対処として、どこを確認すればよいのでしょうか?

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

失礼いたします。

Re: R merge()関数でのエラーメッセージへの対処方法
  投稿者:青木繁伸 2016/08/25(Thu) 07:45 No. 22123
あなたのやったことと環境を再現できそうにないので(私のパソコンには EXcel はない),コメントはできません
CSV で読めるのなら,CSV で保存すればよいと思います
ほかの人とデータを共用するのでなければ,Excel を介在する意味はないと思います

【御礼】 Re: R merge()関数でのエラーメッセージへの対処方法
  投稿者:明石 2016/08/25(Thu) 12:14 No. 22124
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

貴重なコメントをありがとうございました。
御礼を申し上げます。

お手数をお掛けして、大変に申しわけございませんでした。
失礼をいたします。


就職に最も役立つ生物統計学
  投稿者:ひの 2016/08/23(Tue) 09:00 No. 22121
ずいぶんご無沙汰しております。質問というわけではありませんが、この記事によると、
http://forbesjapan.com/articles/detail/13256

2016年のランキングで「就職に最も役立つ修士号」に輝いたのは、生物統計学なのだそうです。日本ではどうなのかな?


分散ゼロの場合の平均値の比較
  投稿者:ヒロ 2016/08/18(Thu) 23:22 No. 22111
下記の2つの事例に対する判断の成否を教えてください。いずれものある生物のバイオマス増加量データに関するもので、A,B群(いずれも反復数は3)の平均値に有意差があるか知りたいというのがもともとの動機です。各々の判断は正しいでしょうか?

事例1
A群:20、20、25; B群:10、10、10
B群の分散がゼロであり、パラメトリックの検定はできない。ノンパラを行うにしては反復が少なすぎて有意差はそもそも出ない。よって、2群の間に有意差があるかの判断を保留し、反復数を増やして再試験を行う。

事例2
A群:20、20、25; B群:0、0、0
B群の平均値、分散がともにゼロであり、パラメトリックの検定はできない。ノンパラを行うにしては反復が少なすぎて有意差は出ない。しかし、検定はできないが、B群は明らかにバイオマス増加量がゼロなので、A群の方がバイオマスの増加が大きい、と判断する。

アドバイスをいただければ幸いです。
よろしくお願いいたします。

Re: 分散ゼロの場合の平均値の比較
  投稿者:青木繁伸 2016/08/19(Fri) 10:03 No. 22113
> B群の分散がゼロであり、パラメトリックの検定はできない。
> B群の平均値、分散がともにゼロであり、パラメトリックの検定はできない。

そんなことないでしょう?実際に計算してみましたか?
検定の計算式を見ても,どちらかの分散が0であっても,計算できる(0での割り算は発生しない)。たとえば,
http://aoki2.si.gunma-u.ac.jp/lecture/Average/t-test.html
にある式をみてみましょう。
R での計算結果も示しておきましょうか。
> t.test(c(20, 20, 25), c(10, 10, 10))

Welch Two Sample t-test # ウェルチの方法

data: c(20, 20, 25) and c(10, 10, 10)
t = 7, df = 2, p-value = 0.0198 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
4.495579 18.837755
sample estimates:
mean of x mean of y
21.66667 10.00000

> t.test(c(20, 20, 25), c(10, 10, 10), var.equal=TRUE)

Two Sample t-test # 等分散を仮定

data: c(20, 20, 25) and c(10, 10, 10)
t = 7, df = 4, p-value = 0.002192 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
7.039258 16.294075
sample estimates:
mean of x mean of y
21.66667 10.00000

> t.test(c(20, 20, 25), c(0, 0, 0))

Welch Two Sample t-test # ウェルチの方法

data: c(20, 20, 25) and c(0, 0, 0)
t = 13, df = 2, p-value = 0.005865 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
14.49558 28.83775
sample estimates:
mean of x mean of y
21.66667 0.00000

> t.test(c(20, 20, 25), c(0, 0, 0), var.equal=TRUE)

Two Sample t-test # 等分散を仮定

data: c(20, 20, 25) and c(0, 0, 0)
t = 13, df = 4, p-value = 0.000202 # 有意
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
17.03926 26.29408
sample estimates:
mean of x mean of y
21.66667 0.00000
しかし,いずれの場合も,サンプルサイズが少ないことにクレームが付く可能性が大。

ついでに,

> B群は明らかにバイオマス増加量がゼロ

私は,「明らかにゼロ」とは思いませんが...

Re: 分散ゼロの場合の平均値の比較
  投稿者:ヒロ 2016/08/19(Fri) 11:41 No. 22114
青木先生
お返事ありがとうございます。分散ゼロの場合は検定できないと決めてかかっていました。すみません。計算式を見て納得しました。

関連してもう1点お尋ねしたいのですが、3群以上の多重比較を行う場合で、分散が他の群と異なるものがある場合、2群ずつウェルチのt検定を行い、ボンフェローニで補正したP値(3群の場合はP=0.05→P=0.017など)に照らして有意性を判断する、というやり方は適正でしょうか?

再びですが、よろしくお願いいたします。

Re: 分散ゼロの場合の平均値の比較
  投稿者:青木繁伸 2016/08/19(Fri) 21:11 No. 22116
分散が異なろうが同じであろうが,Welch の方法を使えば良いというのは,実例でも裏付けのあることです。

その上で,多重比較を行えばよいというのも,よく知られていることです。

Re: 分散ゼロの場合の平均値の比較
  投稿者: 2016/08/20(Sat) 09:19 No. 22119
3群以上の場合にはGames-Howell法を試してはどうでしょうか。

- http://aoki2.si.gunma-u.ac.jp/lecture/Average/Tukey1.html
- http://aoki2.si.gunma-u.ac.jp/R/tukey.html

Re: 分散ゼロの場合の平均値の比較
  投稿者:ヒロ 2016/08/20(Sat) 23:04 No. 22120
青木先生、荒様、
ご返信、ありがとうございます。

荒様ご紹介の方法は知りませんでした。試してみます。

今後もご指導、よろしくお願いいたします。


交互作用の検定について
  投稿者:立山 2016/08/16(Tue) 23:23 No. 22105
青木先生 始めまして。cox回帰分析のモデルでのinteractionについて質問させてください。

認知症の発症に関して解析したいのですが、要因A〜Dがあるとして、要因B〜Dで調整した要因Aによるハザード比を出すモデルは、説明変数に要因A〜要因Dを投入して行うことができました。

そして要因A*要因Bの交互作用があるのかも知りたいのですが、モデルに「要因A、要因B、要因A*要因B」の項を投入して要因A*要因Bが有意であるか調べればよいのでしょうか?それとも「要因A〜要因D」と「要因A*要因B」を投入して、要因C要因Dも調整したうえでの交互作用をみるのがよいのでしょうか。
手持ちの本など調べたのですが記載がなく困っております。

お忙しいところ恐縮ですがよろしくおねがいいたします。

Re: 交互作用の検定について
  投稿者:鈴木康弘 2016/08/19(Fri) 07:21 No. 22112
 青木先生からのコメントがないところをみると、お好きなように、ということかもしれませんね。

 主効果に「要因A〜要因D」の4つ入れると、「要因A*要因B」以外の交互作用はどうなんだ?と査読者に言われることはあるかもしれない。

Re: 交互作用の検定について
  投稿者:青木繁伸 2016/08/19(Fri) 21:19 No. 22117
いつも私の意見を期待されるというのも困った状況ですが...

これも,多変量解析は分析に関与するすべての変数に基づいて分析すべしという大原則からすれば,A,B,C,D の主効果,AB, AC, AD, BC, BD, CD, ABC, ABD, ACD, BCD, ABCD の交互作用からなるフルモデルをまず検討すべきでは?もっとも,普通は高次の交互作用の影響度は低いので
A,B,C,D,AB, AC, AD, BC, BD, CD
A,B,C,D,AB, AC, AD, BC, BD, CD, ABC, ABD, ACD, BCD
A,B,C,D,AB, AC, AD, BC, BD, CD, ABC, ABD, BCD, BCD, ABCD
の順でよいと思うのではあるけれども。

Re: 交互作用の検定について
  投稿者:立山 2016/08/20(Sat) 08:43 No. 22118
鈴木様 
アドバイスありがとうございます。
特にA×Bについて興味があったので調べてみようと考えたのですが、それ以外の交互作用についても検討をしたいと思います。

青木先生
アドバイスありがとうございます。
アドバイスを参考にさせていただき、優先度の高い順で検討したいと思います。


単変量を積み重ねても多変量にならない件
  投稿者:青木繁伸 2016/08/18(Thu) 18:53 No. 22108
以下で用いるデータセット zip 形式で圧縮,エンコーディングは utf-8

Download:22108.zip 22108.zip


Re: 単変量を積み重ねても多変量にならない件
  投稿者:青木繁伸 2016/08/18(Thu) 19:12 No. 22109
上のデータを使って分析した結果を吟味してみよう。
> d = read.table("glm.dat", header=TRUE, fileEncoding="euc-jp")
相関係数を見ると,y と x1, X2, x4 はかなりの相関があるが,x3 とはほとんど無相関(それぞれ散布図を描いてみるとよい)。
> cor(d)
x1 x2 x3 x4 y
x1 1.00000000 0.6130828 -0.08943241 0.4913103 0.5505327
x2 0.61308282 1.0000000 0.30814005 0.8306170 0.6977745
x3 -0.08943241 0.3081401 1.00000000 0.2534457 -0.0259393
x4 0.49131031 0.8306170 0.25344571 1.0000000 0.6663762
y 0.55053270 0.6977745 -0.02593930 0.6663762 1.0000000
x1 〜 x4 を用いて単変量ロジスティック回帰を行う(出力の一部を省略)
x1, x2, x4 が有意
> summary(glm(y ~ x1, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -11.33739 2.34670 -4.831 1.36e-06
x1 0.19621 0.04275 4.590 4.43e-06
Residual deviance: 80.772 on 98 degrees of freedom
> summary(glm(y ~ x2, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -19.6555 4.0133 -4.898 9.70e-07
x2 0.3542 0.0746 4.748 2.06e-06
Residual deviance: 51.258 on 98 degrees of freedom
> summary(glm(y ~ x3, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.65418 1.13867 -0.575 0.566
x3 -0.00582 0.02245 -0.259 0.795 有意ではない
Residual deviance: 118.52 on 98 degrees of freedom
> summary(glm(y ~ x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -14.13175 2.85497 -4.950 7.43e-07
x4 0.25013 0.05271 4.745 2.09e-06
Residual deviance: 65.143 on 98 degrees of freedom
x1,x2,x4 を使って分析すると
> summary(glm(y ~ x1+x2+x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -36.07419 9.16181 -3.937 8.24e-05
x1 0.22755 0.08223 2.767 0.00566
x2 0.27575 0.09711 2.839 0.00452
x4 0.15913 0.07817 2.036 0.04178
Residual deviance: 35.413 on 96 degrees of freedom AIC: 43.413
そらみろ,ちゃんと全部有意になったじゃないかと,思いきや
総当たり法でやってみて,AIC が小さい順に挙げると,x1, x2, x4 を含むモデルは 5 番目。
上位 4 モデルは「すべて x3 が入っている!!」
x3 って,y と相関なかったよね?
> print(all.logistic(d), sort.by="AIC")

deviance AIC Formula
26.56728 36.56728 y ~ x1 + x2 + x3 + x4
29.30262 37.30262 y ~ x1 + x2 + x3
33.52754 41.52754 y ~ x2 + x3 + x4
35.80734 41.80734 y ~ x2 + x3
35.41286 43.41286 y ~ x1 + x2 + x4
40.34653 46.34653 y ~ x1 + x2
43.94204 51.94204 y ~ x1 + x3 + x4
45.99639 51.99639 y ~ x1 + x4
48.08899 54.08899 y ~ x2 + x4
51.25798 55.25798 y ~ x2
58.07721 64.07721 y ~ x3 + x4
65.14314 69.14314 y ~ x4
80.77157 84.77157 y ~ x1
80.65916 86.65916 y ~ x1 + x3
118.52330 122.52330 y ~ x3
全部使うモデルはいくらなんでも
それに x4 は有意じゃないよ
> summary(glm(y ~ x1+x2+x3+x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -36.24623 10.78017 -3.362 0.000773
x1 0.22241 0.10053 2.212 0.026950
x2 0.49536 0.17206 2.879 0.003989
x3 -0.21480 0.09236 -2.326 0.020037
x4 0.14700 0.09586 1.533 0.125160
Residual deviance: 26.567 on 95 degrees of freedom AIC: 36.567
じゃあ,x4 なんか除いてしまえ。
いいんじゃないですか。
> summary(glm(y ~ x1+x2+x3, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -31.04391 9.05752 -3.427 0.000609
x1 0.20348 0.09389 2.167 0.030219
x2 0.56413 0.15877 3.553 0.000381
x3 -0.21993 0.08598 -2.558 0.010529
Residual deviance: 29.303 on 96 degrees of freedom AIC: 37.303
その他のモデルもやってみよう。
> summary(glm(y ~ x2+x3+x4, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -22.31070 5.92023 -3.769 0.000164
x2 0.49410 0.14808 3.337 0.000848
x3 -0.21307 0.07859 -2.711 0.006705
x4 0.11349 0.07827 1.450 0.147039
Residual deviance: 33.528 on 96 degrees of freedom AIC: 41.528
> summary(glm(y ~ x2+x3, data=d, family=binomial))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -19.33143 4.71465 -4.100 4.13e-05
x2 0.54806 0.13307 4.119 3.81e-05
x3 -0.21356 0.07522 -2.839 0.00452
Residual deviance: 35.807 on 97 degrees of freedom AIC: 41.807
じつは,このデータセット,他に x5 というのがあるのだけどこれが強力で,x5 一つだけで十分予測ができてしまうというシロモノ。分析に使った(使えた)変数だけを弄っているだけでは,そんな強力な変数があるなんて,お釈迦様でも気がつかないわけです。

Re: 単変量を積み重ねても多変量にならない件
  投稿者:青木繁伸 2016/08/18(Thu) 19:21 No. 22110
2群の判別分析の場合に,2群で平均値の全く同じ変数を使うと判別率が上がるというのがわかりやすく図示している
http://aoki2.si.gunma-u.ac.jp/lecture/Discriminant/ans2.html
分析例
http://aoki2.si.gunma-u.ac.jp/BlackBox/how2use/dataset2.html

つまり,2群の平均値の差が有意な変数を使って判別分析をするのは NG ということ

重回帰分析の場合も同様

Re: 単変量を積み重ねても多変量にならない件
  投稿者:さくらとお城 2016/08/19(Fri) 16:07 No. 22115
青木先生,いつもありがとうございます。大変わかりやすい例示で,よくわかりました。


単変量で有意差なし、多変量で有意
  投稿者:内科医8 2016/08/10(Wed) 10:03 No. 22096
青木繁伸先生、みなさま
上記に関して質問させて頂きます。
AはZに対して正の関連をするのではないかとの仮説に基づいた横断研究において、AとZとの関連を見るときに、単変量解析では有意な相関関係はありませんでしたが、多変量解析ではAはZに対して独立した正の関連因子でした。
 Aと同時に因子Bに関しても解析をしており、AとBは正の相関を認めるのに対し、BとZは負の相関を認めておりました。先の多変量解析にはBも含まれおり、BはZに対して独立した負の関連因子でした。Bの影響により、AとZの関連が単変量解析ではマスクされていたと考えました。
このような考えは統計学において、認められる考えなのでしょうか。また、このような現象を示す統計学用語などありますでしょうか。ご教授のほど、よろしくお願い申し上げます。

Re: 単変量で有意差なし、多変量で有意
  投稿者:鈴木康弘 2016/08/15(Mon) 07:02 No. 22100
 これまで青木先生が何度も書かれているように、単変量解析では有意な相関はないが、多変量解析では有意な相関がみられるというのは、何の不思議もありません。

 このような現象を示す統計学用語は?..さあ、わかりません。

Re: 単変量で有意差なし、多変量で有意
  投稿者:さくらとお城 2016/08/17(Wed) 11:22 No. 22106
横からすみません。

>これまで青木先生が何度も書かれているように、単変量解析では有意な相関はないが、多変量解析では有意な相関がみられるというのは、何の不思議もありません。

確かに青木先生が書かれているコメントは,私も何回か読んだことがあります。ですが,”どうして”と突っ込まれると,返事に困っている状況です。何かわかりやすい説明とか,具体的な例示がありませんか?

特に単回帰で相関の高い変数を選んでも,重回帰で必ずしも良くならない例示などです。どうぞよろしくお願いします。

Re: 単変量で有意差なし、多変量で有意
  投稿者:青木繁伸 2016/08/18(Thu) 18:46 No. 22107
別記事で


統計方法について、ご教示ください。
  投稿者:初心者 2016/08/15(Mon) 07:43 No. 22101
青木先生、いつもお世話になります。

初歩的な質問で、申し訳ありません。

私は、認知症発症の危険因子を調べるために、地域における65歳以上の高齢者を対象とした疫学調査を行いました。初回の疫学調査では、65歳以上の高齢者住民1000人のうち、650人を健常高齢者(認知症ではない高齢者)と診断し、高血圧や糖尿病、脂質異常症などの内科的疾患、うつ病などの既往について聴取しました。そして、その7年後に第2回目の疫学調査を行ったところ、初回調査で健常高齢者と診断した650人のうち、35人が認知症を発症していました。

そこで、ご質問があります。

―蕾鵑猟敢困之鮠鏐睥霄圓反巴任気譟2回目の調査では認知症と診断された人のうち、初回の調査時に高血圧や糖尿病、脂質異常症、うつ病などの既往がなくても、2回目の調査時までの7年間の間に、これらの疾患を発症された方が存在します(2回目の調査時に聴取しました)。これらの疾患が、認知症発症の危険因子となったかどうか調べるための統計方法として、初回調査時における疾患の存在の有無と2回目調査時における疾患の存在の有無について、おのおのにロジスティック回帰分析を行うという方法でよろしいのでしょうか?

認知症を発症した人が35人とサンプルサイズが小さいために、本当は危険因子となる疾患であるにもかかわらず、ロジスティック回帰分析で有意差がでない可能性も考えられると思います。その場合、effect sizeを求めたらいいのでしょうか?

初歩的な質問で申し訳ありませんが、何卒よろしくお願い致します。

Re: 統計方法について、ご教示ください。
  投稿者:青木繁伸 2016/08/15(Mon) 22:30 No. 22102
初回調査で健常高齢者と診断した650人を対象としてロジスティック回帰を行うということでしょう。それ以外ありますか?
有意差が出ないから effect size をもとめるということではないでしょう。排反じゃないでしょう。ロジスティック回帰分析をする,effect size も示すということでしょう?

Re: 統計方法について、ご教示ください。
  投稿者:初心者 2016/08/16(Tue) 13:16 No. 22104
青木先生、御教示ありがとうございました。
今一度、自分なりに勉強してみます。
これからもよろしくお願い申し上げます。


R 文字列の分割 strsplit()関数
  投稿者:明石 2016/08/10(Wed) 22:00 No. 22097
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。
どうぞ、よろしくお願いいたします。

−−−

文字列"_"で連結した文字列ベクトル dat があります(膨大な長さ)。

(例)dat
"123_aaa"
"4567_bbbb"
"787_ccc"
"12345_ddddd"
"56787_boeing"
,,,,,,

やりたいことは、以下の通りです。

(1) まず、文字列"_"の左右で分割して、2つのベクトルを作成します。

  左側のベクトル名を lhs とします。データ型は整数
  (例) 123,4567,787,12345,56787,,,,,

  右側のベクトル名を rhs とします。データ型は文字列
  (例) "aaa","bbbb","ccc","ddddd","boeing",,,,,

(2) 次に、lhs、rhsの2列からなるデータフレーム(変数名は df )を作成したいと思います。

lhs,rhs
123,"aaa"
4567,"bbbb"
787,"ccc"
12345,"ddddd"
56787,"boeing"
,,,,,,

以下、トライしてみましたが、listに負慣れなために出来ないで苦慮しております。

dat <- c("123_aaa","4567_bbbb","787_ccc","12345_ddddd","56787_boeing")

ds <- strsplit(dat, "_")

ds が listになることは分かりましたが、listは負慣れなために、要素の取り出しができないでおります。

要素の取り出しができ、lhsベクトル、rhsベクトルを生成できれば、データフレーム df を作成することができます。

上記の手順に沿って、データフレーム df までを作成するRプログラムを、ご教示いただければ、
大変に助かります。

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

Re: R 文字列の分割 strsplit()関数
  投稿者:青木繁伸 2016/08/11(Thu) 08:26 No. 22098
リストを取り出すのは [[ ]] です
> ds[1]
[[1]]
[1] "123" "aaa"

> ds[2]
[[1]]
[1] "4567" "bbbb"

> ds[[2]]
[1] "4567" "bbbb"
> ds[[2]][1]
[1] "4567"
> ds[[2]][2]
[1] "bbbb"
リストのまま処理することもできますけど,strsplit の戻り値を unlist すれば 1 つのベクトルになります。
> unlist(ds)
[1] "123" "aaa" "4567" "bbbb" "787" "ccc" "12345" "ddddd" "56787"
[10] "boeing"
それ以降はいろいろなやり方があるでしょう。1列目が numeric,2列目は character にしたいというのでちょっと面倒くさいですが。
> dat <- c("123_aaa","4567_bbbb","787_ccc","12345_ddddd","56787_boeing")
> ds <- strsplit(dat, "_")
> x <- unlist(ds)
> f <- rep(c(TRUE, FALSE), length(dat))
> lhs <- as.numeric(x[f])
> rhs <- x[!f]
> df <- data.frame(lhs, rhs, stringsAsFactors=FALSE)
> df
lhs rhs
1 123 aaa
2 4567 bbbb
3 787 ccc
4 12345 ddddd
5 56787 boeing
> sapply(df, class)
lhs rhs
"numeric" "character"
楽しんでください。

【御礼】Re: R 文字列の分割 strsplit()関数
  投稿者:明石 2016/08/11(Thu) 09:39 No. 22099
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

大変に有難いご教示をいただき、誠にありがとうございました。

今回も、助けていただきました。

listの扱いかたが見えてきましたので、
この夏休み、楽しみながらlistの勉強をさせていただきます。

心から御礼を申し上げます。

Re: R 文字列の分割 strsplit()関数
  投稿者: 2016/08/16(Tue) 12:20 No. 22103
少し遅くなりましたが、
以下のようにしてもデータフレーム化できました。
(簡単なものしか回答できていませんね)

> dat <- c("123_aaa","4567_bbbb","787_ccc","12345_ddddd","56787_boeing")
> ds <- strsplit(dat, "_")

> library(plyr)
> df <- ldply(ds)
> df
V1 V2
1 123 aaa
2 4567 bbbb
3 787 ccc
4 12345 ddddd
5 56787 boeing

> colnames(df) <- c("lhs","rhs")
> df$lhs <- as.numeric(df$lhs)


評価法の作成方法
  投稿者:山里 2016/08/06(Sat) 03:29 No. 22089
評価法の作成方法についてご質問させて頂きます。
ある症状の重症度の予測を行える評価法を作成したいと考えています。
評価法には10種類以上の検査結果を統合し、100点満点などで採点できるようにしたいと考えています。
しかし、それぞれの検査結果が予測に与える影響度が違うため、単純にすべての項目を10点満点で採点するのではなく、ある項目は5点満点、ある項目は20点満点と重みを変えるべきと思っています。
統計を用いてこのような評価法を作成することは可能でしょうか?
ご教授の程、宜しくお願い致します。

Re: 評価法の作成方法
  投稿者:青木繁伸 2016/08/06(Sat) 06:23 No. 22090
項目群ごとに重回帰分析を行って項目の選定と重み付けをするということを繰り返し,更に統合された数値で重回帰分析をするということになるでしょう。

しかし,多変量解析的にいえば,変数セットを変えて重回帰分析を繰り返すというのはさけるべきでしょう。

この掲示板にも繰り返し出てきている,「単変量ロジスティック回帰で有意な変数を探索し,最後にそれら変数を使って多変量ロジスティック回帰を行うというのは,多変量解析的に誤った方法である」というのと,根本的に同じやり方だからです。

Re: 評価法の作成方法
  投稿者:山里 2016/08/06(Sat) 06:39 No. 22091
青木先生
早速のご回答ありがとうございます。
「項目群ごとに」とのことですが、どの項目とどの項目を合わせて重回帰分析を行うかはどのように判断すればよろしいでしょうか?
また、「統合された数値で」とはどのような数値のことなのでしょうか?
初歩的な質問で大変申し訳ございません。
宜しくお願い致します。

Re: 評価法の作成方法
  投稿者:青木繁伸 2016/08/06(Sat) 21:44 No. 22092
深読みしたようです

要するに,すべての関連する項目を持ちいて,変数選択によって重回帰分析を行えば良いというだけのことですね。

Re: 評価法の作成方法
  投稿者:山里 2016/08/07(Sun) 07:30 No. 22093
ご回答ありがとうございます。
可能であれば、重回帰式ではなく以下のような評価シートを作成したいのですが難しいでしょうか?

________20点___10点____0点___
歩行時間__5秒以下_10秒以下__10秒以上__
歩行距離___300m__200m____100m___
意欲_____あり__________なし___
_・______・___・______・___
_・______・___・______・___
→合計点数80点以上 歩行生活に支障なし
50点以上 一部介助が必要
40点以下 歩行困難

かなり見にくいですし、例えなのでざっくりですがこのようなチェックしていくことで点数化される評価表を作成したいです。そのために、上記の例で説明しますと、歩行時間を何秒以下なら何点と採点させるかなどはどのように決めるべきでしょうか。
ややこしい質問で申し訳ありませんが、宜しくお願い致します。

Re: 評価法の作成方法
  投稿者:青木繁伸 2016/08/07(Sun) 10:06 No. 22094
重回帰分析の評価を簡単にしたいということであれば,ノモグラムにするとそれほど精度を落とすこともないのでは?
ノモグラムの例は
http://www.dynacom.co.jp/tech/tech_column/tech003.html
などを見ると作り方がわかるでしょう。そのページはロジスティック回帰についてですが,そのページの Linear Predictor が,重回帰分析では予測値になります。各項目の測定値に対応するポイントの求め方,ポイントの合計値から予測値を求める方法は,説明の通りです。

Re: 評価法の作成方法
  投稿者:山里 2016/08/08(Mon) 17:14 No. 22095
青木先生

ご回答ありがとうございます。
ノモグラフというのは知らなかったため、非常に勉強になりました。
しかしながら、ノモグラフというのが私の領域では馴染みがないため、例えば参考ページでいいますとage10〜19=8点、20〜29=20点…というように中央値をポイント化し、質問No.22093のような評価表を作成するのは無理がありますでしょうか?
評価表を作成したいばかりに質問ばかりして大変申し訳ありません。


残差分析について
  投稿者:コロン 2016/07/28(Thu) 17:27 No. 22085
いつもお世話になっております。

基本的なことかと思いますが、教えてください。

2×2のカイ二乗検定においてp値が0.05未満の場合、どのセルが影響を与えているかをチェックする残差分析を行う必要はありますでしょうか。

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

Re: 残差分析について
  投稿者:青木繁伸 2016/07/28(Thu) 21:59 No. 22086
必要なら行う,そうでなければ行わなくてよい...のでは?

Re: 残差分析について
  投稿者:コロン 2016/07/29(Fri) 08:53 No. 22088
青木先生

いつもつまらない質問で申し訳ございません。また最初の投稿の内容が不十分で申し訳ございませんでした。

実は査読者から「2×2の場合は,度数の偏りは見たらわかるので残差分析は不要」とあり,このコメントに疑問を感じお尋ねした次第です。

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


重回帰分析の従属変数の正規性について
  投稿者:睡眠専門医 2016/07/27(Wed) 23:03 No. 22083
いつもお世話になっております。
お尋ねしたいのですが、重回帰分析の従属変数として用いたい変数が分布に偏りがあります。そこで、変数変換をしたところ、3乗根変換にて正規分布(Kolmogorov-Smirnov の正規性の検定による)に近い形となりました。
重回帰分析の従属変数は、正規分布である必要があるのでしょうか。
また、独立変数についても同様でしょうか。
以上につきましてご教示いただけましたら幸いです。

Re: 重回帰分析の従属変数の正規性について
  投稿者:青木繁伸 2016/07/28(Thu) 09:21 No. 22084
従属変数は正規分布する必要はありません。たとえば,ロジスティック回帰の場合,従属変数は二項分布(0/1の二値データ)です。このように,従属変数が二値データの場合には普通の重回帰分析は適用できません。また,y = ax^b のような場合は,たとえ x が正規分布でも,y は正規分布しませんね。両辺の対数をとって,log(y)=log(a)+b*log(x) のようにすれば,直線回帰(重回帰)に持ち込めます(もっとも,これも,非線形回帰を行う方がよいですけど)。

独立変数も正規分布する必要はないです。たとえば,ダミー変数(これも 0/1 の二値データ)など,正規分布するわけでもないし,データから検量線を引くような場合(単回帰...重回帰に発展)にも独立変数は任意に選ばれた(きれいな整数値)ですしね。

変数変換をするときは,闇雲にやるのではなく理論的根拠に従うべきです。

従属変数と独立変数の間に理論式(のようなもの)があるなら,非線形回帰を試みるのも一法でしょう。

Re: 重回帰分析の従属変数の正規性について
  投稿者:睡眠専門医 2016/07/29(Fri) 03:11 No. 22087
ご教示いただきありがとうございました。
変数変換はやみくもに行うべきではないのですね。
非線形回帰についても勉強させていただきます。


同一サンプルによる統計処理について
  投稿者:アルビー 2016/07/20(Wed) 10:05 No. 22080
青木先生、みなさま

 アルビーです。現在、同一サンプルを用いたオムニバス科目の授業評価について分析を進めており、この関係で質問をしたく書き込ませていただきました。

 オムニバス講義を行った4人の教員に対する授業評価をx個の評価項目を用いてn人の学生が個別に行いました。個々の教員の授業に対する授業評価ではなく、授業科目に対する授業評価(満足度)を検討したかったため、収集された結果を、x個×4n人のデータ行列を作成して、授業評価に影響する授業要素(教材や教室環境)の関係をロジスティック回帰分析により分析しました。
 得られた結果を先行研究と比較したところ、妥当な予測が見出されました。
 しかしながら、刺激対象(教員)は独立していますが、同一サンプルの反応ということからすると、何らかの不都合が起こるとも考えられます。例えば、当該サンプルの特徴が過度に誇張されると推察されます。
 また、サンプル構造そのものに問題があるとしたら、どのようなことなのかお教えいただけませんでしょうか?

Re: 同一サンプルによる統計処理について
  投稿者:青木繁伸 2016/07/20(Wed) 21:03 No. 22081
完全に学生が把握されているなら,「x個×4n人のデータ行列」ではなく「4x個×n人のデータ行列」で扱うべきでしょう。なんといっても,サンプルサイズは4nではなくnなのですから(サンプルサイズの水増しは厳禁。検定はサンプルサイズにも影響されるんだから。そもそも,データは独立でないといけない)。「個々の教員の授業に対する授業評価ではなく」というのなら,xについて4人の教員への評価点の合計(または平均)として「x個×n人のデータ行列」で扱うのが妥当ではないですか?

Re: 同一サンプルによる統計処理について
  投稿者:アルビー 2016/07/21(Thu) 19:41 No. 22082
青木先生

お返事感謝申し上げます。
「4x個×n人のデータ行列」で扱うべき」とのご指摘、改めて、反省いたしました。
学生数がどうしても少ない授業でしたので、苦肉の策として考えた方法ではありますが、再度練り直したいと思います。4名の教員の授業の評価得点を他の手法で作成してみたいと思います。
本当に、ありがとうございました。


R ファイルの突合せ
  投稿者:明石 2016/07/18(Mon) 08:55 No. 22077
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rプログラムについて、ご教示いただきたいことが出てきました。
どうぞ、よろしくお願いいたします、

------
添付の画像をご覧ください。

お示しをしました例題は、グラフネットワークのファイルです。

da1 は、ノードリストです。
id,label,score の3列からなります。

da2 は、エッジリストです。
node1,node2,weight の3列からなります。

やりたいことは、
dat2 の node1,node2 を、
dat1 と突合せをして、label を id に置き換えた dat3 を生成出力することです。

動くプログラムはできましたが、
ノードリスト、エッジリストのサイズが膨大なので、
処理効率、R言語の良さ(ベクトル化処理、論理ベクトルなど)の観点から、
もっと良いプログラムを勉強したいと思っています。

良いお手本をお示しいただけると、大変に助かります。
どうぞ、よろしくお願いいたします。

【プログラム例】

dat1 <- read.table("dat1.csv",header=T, sep="," ,stringsAsFactors=F)

dat2 <- read.table("dat2.csv",header=T, sep="," ,stringsAsFactors=F)

v <- dat1$label

dat3 <- dat2

for(i in 1:nrow(dat2)) {

n1 <- which(is.element(v, dat2$node1[i]))
dat3$node1[i] <- dat1$id[n1]

n2 <- which(is.element(v, dat2$node2[i]))
dat3$node2[i] <- dat1$id[n2]

}

dat3

> dat3
node1 node2 weight
1 1 2 0.70
2 1 3 0.64
3 1 6 0.55
4 1 9 0.44
5 2 3 0.88
6 2 5 0.50
7 2 8 0.41
8 3 4 0.66
9 3 7 0.60
10 3 8 0.55
11 3 9 0.39
12 4 8 0.38
13 4 9 0.33


Re: R ファイルの突合せ
  投稿者:青木繁伸 2016/07/18(Mon) 13:56 No. 22078
もとのプログラムも十分早いと思いますが,
data.frame(node1=apply(outer(dat1$label, dat2$node1, "=="), 2, which),
node2=apply(outer(dat1$label, dat2$node2, "=="), 2, which),
weight=dat2$weight)

data.frame(node1=outer(dat2$node1, dat1$label, "==") %*% seq_along(dat1$label),
node2=outer(dat2$node2, dat1$label, "==") %*% seq_along(dat1$label),
weight=dat2$weight)

data.frame(node1=charmatch(dat2$node1, dat1$label),
node2=charmatch(dat2$node2, dat1$label),
weight=dat2$weight)
のようなものも。
どれがどれくらい速いか試してみてください。

【御礼】Re: R ファイルの突合せ
  投稿者:明石 2016/07/18(Mon) 18:36 No. 22079
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

良いお手本を3つもお示しくださり、大変に助かりました。

apply関数を使う方法を期待していましたので、とても勉強になります。

まだ使ったことのないouter関数があり、良い勉強になります。

お休みのところ、大変に有難いご教示をいただき、誠にありがとうございました。


多重クロス表(参加者内計画をふくむ)について
  投稿者:ふくしま 2016/07/12(Tue) 20:07 No. 22069 HomePage
初めて投稿させていただきます。
参加者内計画を含む多重クロス表の分析方法について質問させていただきます。
2011年にも似たような質問をされていた方がいらっしゃいましたが(お家のマークより,URLをご参照ください),添付させていただいた図のような多重クロス表分析には,階層対数線形分析を使えば良いのではないかと指導教員の指摘を受けましたので,階層対数線形分析で分析ができないか,もしご存知の方がいらっしゃればご教授願えますでしょうか。それとも,以前の質問に対するご回答にあったように,やはり分散分析を用いたほうが良いのでしょうか?

実験デザインは,手がかりの有無(手がかりあり vs. 手がかりなし:参加者間計画)× 刺激呈示方法の種類(A vs.B:参加者間計画)× 試行の繰り返し(1回/2回:参加者内計画)の3要因混合計画です。数値はそれぞれの試行で「誤答」,「正棄却」,「わからない」反応をした人数です。

すべての独立変数が参加者間計画である場合の分析方法はわかったのですが,対数線形分析では参加者間計画を前提とした分析手法しか見当たらず,参加者内要因を含む分析手法について参考にできるものがありませんでしたので,こちらに質問させていただきました。

拙い説明で申し訳有りませんが,不十分な点は捕捉させていただきますので,どうぞ宜しくお願い致します。


Re: 多重クロス表(参加者内計画をふくむ)について
  投稿者:鈴木康弘 2016/07/14(Thu) 07:16 No. 22073
 頼りにならない意見ですが..

 分散分析はできないんじゃないでしょうか。
 引用された発言では「高い」「低い」の2値を0,1と見なして分散分析ができたわけですが、ふくしまさんの反応は「誤答」「正棄却」「わからない」の3つあるので。

 対数線形モデルのことはなおさらわからないのですが..

 階層対数線形混合モデルというのがあればいいんでしょうけれど。
 参加者内計画を参加者間計画と強引にみなして、階層対数線形分析をして有意であれば
参加者内計画でも有意であると言っていいんじゃないでしょうか。

Re: 多重クロス表(参加者内計画をふくむ)について
  投稿者:ふくしま 2016/07/14(Thu) 13:06 No. 22074
鈴木様

ご意見をありがとうございます。
なるほど,参加者内計画を参加者間計画とみなして分析してみるというのは,試みたことがなかったなのでやってみたいと思います。

もう一つ追加でご意見をいただけるとありがたいのですが,今回の参加者内要因である「試行の繰り返し(1回/2回)」の効果は表を見る限り,ないだろうと推測しています。
そこで,すべての条件(手がかりあり×呈示方法A,手がかりなし×呈示方法A.手がかりあり×呈示方法B,手がかりなし×呈示方法B.の4条件)で「試行の繰り返しによる効果はない→判断を試行間で変える人は,変えない人よりも少ない」ということを何らかの方法で示したうえで(判断を2試行間で変えた人数と,変えなかった人数に対して比率の差の検定をしてみる・・とかでしょうか),「手がかりの有無」× 「刺激呈示方法の種類」(2×2)の参加者間要因のみの分析にしてしまうというのは,やはり分析上問題でしょうか?乱暴な方法だとは承知しているのですが・・

追加で申し訳有りませんが。ご意見をいただけると幸甚です。
どうぞ宜しくお願い致します。

Re: 多重クロス表(参加者内計画をふくむ)について
  投稿者:鈴木康弘 2016/07/15(Fri) 07:11 No. 22075
 データを1回目のだけにするなら、それでよいと思います。
同じ人に2回やってデータを増やしました、というのは反則でしょうね。

Re: 多重クロス表(参加者内計画をふくむ)について
  投稿者:ふくしま 2016/07/15(Fri) 13:12 No. 22076
鈴木様

早速のご回答ありがとうございました。
階層対数線形分析ですべての変数を参加者間としてみなして分析してみたところ,やはり繰り返しの効果は最終モデルに含まれませんでした。

繰り返しの効果を潰して,「手がかりの有無」× 「刺激呈示方法の種類」(2×2) の分析を1回目,2回目それぞれで行うかどうか,それとも1回目だけ(あるいは2回目だけ)のデータに行うか,分析方法の妥当性も含めて,再度指導教員と相談してみたいと思います。

匙を投げられることが多かったので,聞いていただき大変助かりました。
ありがとうございました。


多重ロジスティックの解析について
  投稿者:安達 2016/07/11(Mon) 22:28 No. 22067
多重ロジスティック回帰分析を行うにあたり、投入する独立変数の選択方法についてご質問させて頂きます。
独立変数を選択するにあたり、『2群間比較を行い、有意差を認めた変数間において多重共線性を考慮して独立変数とするか検討すべき。』という方と、『ロジスティック回帰分析は前述の内容を踏まなくても自動的に処理されるから関係ありそうな変数全てをいきなり独立変数としても良い。』という方がいます。
二人とも指導教官であり、どちらが正解なのか分からなくなりました。
根拠と共に正解をご教授頂けますと幸いです。
宜しくお願い致します。

Re: 多重ロジスティックの解析について
  投稿者:青木繁伸 2016/07/12(Tue) 09:00 No. 22068
この掲示板にも,何回となく繰り返し出てきています。

医学系学会誌に投稿されている論文でお手本のように採用されている「単変量ロジスティック回帰分析で有意な変数を選んで,それらを使って多重ロジスティック回帰分析をする」とうやりかたは,大間違いです。
単変量ロジスティック回帰分析の代わりに,「2群間比較を行い、有意差を認めた変数」はもっと悪い(分析手法が一貫していない)。

根拠?
多変量解析は,分析に用いた変数セットの中での結果である。単変量解析の結果を積み上げても,多変量解析の結果にはならない(だって,単変量解析で済むなら,多変量解析は不要でしょ)

たとえば,「2 群で平均値が全く同じである変数は判別に役立つだろうか。」の答えは,
http://aoki2.si.gunma-u.ac.jp/lecture/Discriminant/ans2.html

> 二人とも指導教官であり

困ったものです。片方の人は,指導教員の資格なしですね。もう一方も,自分の方が正しいと相手を説得できないのかな?

Re: 多重ロジスティックの解析について
  投稿者:安達 2016/07/12(Tue) 20:28 No. 22070
早速のご回答ありがとうございます。
自分なりにも調べ、多重ロジスティックを行うにあたり、いきなり関係しそうな因子を独立変数として投入しても問題ないと確信できました。
あと一点ご質問です。多重共線性に関しましては、多重ロジスティクの中で自動で処理されているため、改めて何か確認であったりする事は必要ないのでしょうか?医学論文では当たり前のようにこちらで除外などの処理をしており、行う事が正解なのか分からなくなりまして…
ご教授宜しくお願い致します。

Re: 多重ロジスティックの解析について
  投稿者:青木繁伸 2016/07/13(Wed) 06:09 No. 22071
> 多重ロジスティクの中で自動で処理されているため

適切な分析プログラムでは「完全な多重共線性」の場合は,問題のある変数を除去してくれますが,そこまでに至らない場合には自動処理はされません。
結果を見て,正しい判断を下す必要があるでしょう。

Re: 多重ロジスティックの解析について
  投稿者:安達 2016/07/13(Wed) 06:34 No. 22072
ご回答頂き、ありがとうございました。
先生のおかげですっきりする事ができました。
今後とも宜しくお願い致します。


重回帰分析の独立変数について
  投稿者:睡眠専門医 2016/07/08(Fri) 23:28 No. 22058
お世話になっております。
お尋ねですが、重回帰モデルに、互いに相関が比較的強い(r=0.5)2つの変数AとBを投入したところ、Bの標準化回帰係数が負となり、単相関(r=0.5)とは符号が異なりました。
ちなみにBは有意な独立変数とはなりませんでした。
符号が逆転したことは、AとBの多重共線性の問題の由来するものなのでしょうか。
VIFでの共線性診断では、AとBともに1.4前後と問題はないようでした。
さらに、Bの部分寄与率を計算したいのですが、標準化回帰係数×単相関係数で計算すると負となります。これは部分寄与率としておかしいようにみえるのですがいかがでしょうか。
以上につきましてご教示いただけましたら幸いです。

Re: 重回帰分析の独立変数について
  投稿者:青木繁伸 2016/07/09(Sat) 07:23 No. 22059
A,B と従属変数の,相関係数行列(三行三列,小数点以下五桁)とサンプルサイズ(n)を教えてください。

Re: 重回帰分析の独立変数について
  投稿者:睡眠専門医 2016/07/09(Sat) 13:48 No. 22062
お世話になります。
重回帰モデルにはA、Bの2つ以外に、C,Dという独立変数が投入されております。それも含めた相関行列をお示しします。(相関係数は小数点以下3桁となっています)
サンプルサイズは109です。
よろしくお願い申し上げます。
		      A     B   C    D	従属変数
A Pearson相関係数 1 .480** .041 .174 .378**
有意確率 (両側) .000 .670 .070 .000
N 109 109 109 109 109

B Pearson相関係数 .480** 1 .176 .251** .235*
有意確率 (両側) .000 .067 .008 .014
N 109 109 109 109 109

C Pearson相関係数 .041 .176 1 .134 .224*
有意確率 (両側) .670 .067 .165 .019
N 109 109 109 109 109

D Pearson相関係数 .174 .251** .134 1 .395**
有意確率 (両側) .070 .008 .165 .000
N 109 109 109 109 109

従属変数Pearson相関係数 .378** .235* .224* .395** 1
有意確率 (両側) .000 .014 .019 .000
N 109 109 109 109 109
**. P<0.01
*. P<0.05


Re: 重回帰分析の独立変数について
  投稿者:睡眠専門医 2016/07/09(Sat) 13:58 No. 22063
数字がずれてしまい申し訳ございません。
A対B  r=0.480
A対従属変数 r=0.378
B対従属変数 r=0.235
でいずれも有意な相関を認めました。

ご教示のほどよろしくお願い申し上げます。

Re: 重回帰分析の独立変数について
  投稿者:青木繁伸 2016/07/10(Sun) 06:12 No. 22064
A,B,C,D を使った分析であると理解しましたが。

確かに,A と B の偏回帰係数の符号が違いますが,
               t 値       P 値 標準化偏回帰係数 トレランス
A 3.4955e+00 0.00069678 0.331929 0.76408
B -3.6786e-01 0.71372401 -0.035928 0.72227
C 2.0443e+00 0.04344687 0.173429 0.95733
D 3.7416e+00 0.00029981 0.323023 0.92436
Y と B の偏相関係数を見ると -0.03604836 となっておりますね。
C と D を一緒にした場合に,B は抑制変数となっているという状況です。
            A           B          C         D         Y
A 0.56236117 0.48000000 0.04100000 0.1740000 0.3780000 上三角行列:単相関係数
B 0.44676154 0.52788829 0.17600000 0.2510000 0.2350000 下三角行列:偏相関係数
C -0.11627986 0.16153367 0.28223074 0.1340000 0.2240000 対角成分:重相関係数
D -0.05115057 0.17671263 0.02228465 0.4304721 0.3950000
Y 0.32424705 -0.03604836 0.19655316 0.3444393 0.5324173
かといって,A, B だけを使った重回帰分析では,A,B の符号は一致するものの,B の偏回帰係数は有意ではなくなりますね。

A があれば B はいらないようなので,A,C,D の3変数を使えば良いように思いますが。
             t 値       P 値 標準化偏回帰係数 トレランス
A 3.7621e+00 0.00027771 0.31585 0.96941
C 2.0198e+00 0.04594799 0.16851 0.98172
D 3.7503e+00 0.00028947 0.31746 0.95360

Re: 重回帰分析の独立変数について
  投稿者:睡眠専門医 2016/07/11(Mon) 00:23 No. 22066
詳細にご教示いただきありがとうございました。
Bについては、その抑制変数としての意味をよく考えた上で対処したいと思います。


偏寄与率
  投稿者:心療内科医 2016/07/04(Mon) 19:42 No. 22053
SPSS で重回帰分析の各説明変数の偏寄与率 partial R2 を出力するにはどうすればよいでしょうか?また各説明変数の目的変数に対する影響を比較するにはこの偏寄与率 partial R2 の大きさを比較するのが適当でしょうか?ご教授よろしくお願い申し上げます。

Re: 偏寄与率
  投稿者:青木繁伸 2016/07/05(Tue) 06:49 No. 22054
偏寄与率?知りませんでしたが,検索すると色々出てきますが,二つめにあった
http://www.niph.go.jp/soshiki/jinzai/download/etc/hotetsu2009.pdf
の 3 ページ左下のスライドに「偏相関係数...重回帰分析の編寄与率の平方根に符号を付けたもの」という記述があるので,偏相関係数を二乗すれば偏寄与率ですね...偏相関は SPSS でも求まったはず。

ただし,一つめに出てきた
http://www.snap-tck.com/room04/c01/stat/stat07/stat0702.html
の (2) 各説明変数の寄与率 という項で,「説明変数全体としての影響力は評価できても、個々の説明変数が目的変数にどれ程の影響を与えているのかということは正確には評価できない」とも書かれています(長いけど,このページは全部読んでおくと吉)

Re: 偏寄与率
  投稿者:scdent 2016/07/05(Tue) 10:12 No. 22055
重回帰分析における要因の影響を比較したいときは、各変数を標準化(平均を0、標準偏差を1に変数変換 → (x-μ)/σ)して、その値を用いて重回帰分析を行った経験があります。

そうして求まった各偏回帰係数は、各変数の1標準偏差あたりの影響度を表すので、要因間の比較が可能となります。

詳しくは、杉本先生の下記「我楽多頓陳館」を参照ください。
http://www.snap-tck.com/room04/c01/stat/stat07/stat0702.html

Re: 偏寄与率
  投稿者:青木繁伸 2016/07/05(Tue) 13:23 No. 22056
> 重回帰分析における要因の影響を比較したいときは、各変数を標準化(平均を0、標準偏差を1に変数変換 → (x-μ)/σ)して、その値を用いて重回帰分析を行った経験があります。

これによって求まるのは「標準化偏回帰係数」です。
各変数を標準化して重回帰分析をやりなおさなくても,もとの変数を使って得られた偏回帰係数と各変数の標準偏差を使って,簡単な計算で標準化偏回帰係数を求めることができます。
http://aoki2.si.gunma-u.ac.jp/lecture/Regression/mreg/mreg2.html
独立変数iの標準化偏回帰係数=独立変数iの偏回帰係数×独立変数iの標準偏差÷従属変数の標準偏差
http://www.snap-tck.com/room04/c01/stat/stat07/stat0702.html
の 注1 にも書いてありますが

Re: 偏寄与率
  投稿者:心療内科医 2016/07/10(Sun) 08:52 No. 22065
お礼が遅くなりましてすみません。各要因の影響度の比較には結局「標準化偏回帰係数」が一般的ですかね。


複数回答とカテゴリー数
  投稿者:コロン 2016/07/07(Thu) 13:38 No. 22057
いつも大変お世話になっております。

「Q あなたは○○の状況で,下のどの言葉を使いますか?」という質問があり,選択肢が5つあったとします。

15年前にも同じ調査がありますので,15年前と今の使用の仕方を比較したいと思います。

ただ15年前のデータは「複数回答」OKとなっているため,被験者数よりもケースの合計は多くなっています。

質問は2点あります。

(1)複数回答は気にせず(つまり集計されている数をそのまま利用して),15年前と今の比較を5つのカテゴリーで,2×5のχ二乗検定で分析して問題ありませんか?

(2)カテゴリーの(選択肢)数ですが,理論上はいくつくらいまで用意することは可能でしょうか?先行研究のある問題には選択肢が13個あり,2×13の検定となるのかなと思うのですが,このような表を分析する際に気を付かなければならないことがあれば合わせて教えていただけますでしょうか。

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

Re: 複数回答とカテゴリー数
  投稿者:鈴木康弘 2016/07/09(Sat) 07:47 No. 22060
頼りにならない意見ですが..

(1)だめじゃないのかな?そんなことができるなら青木先生がすぐ答えてると思うし。
(2)いくつでもOKでしょうが、増やせば増やすほど統計的に有意な結果は出なくなる。

Re: 複数回答とカテゴリー数
  投稿者:コロン 2016/07/09(Sat) 07:51 No. 22061
鈴木先生

回答がつかないので,私の聞き方がまずかったのか,あるいはとんちんかんな質問をしているのか悩んでおりました。先行研究はその先行研究としてさらに遡った研究データを使用して,複数回答があるにも関わらず,χ二乗検定をしていましたので,「あれ?いいのかな」と考えた次第です。

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


消去された質問
  投稿者:青木繁伸 2016/06/28(Tue) 19:12 No. 22052
用が済んだらあとは証拠隠滅か...なんだかな

質疑応答の概要

質問:句読点の使用頻度など,文章の特性が文章作成指導法によって変化するかどうか,被検者にたいして,指導前と指導後に14000字程度の文章を書かせ,特徴をデータ化。例えば,
    読点使用有り 読点使用なし
指導前 回数     回数
指導後 回数     回数
のようになっている。どのような検定を使えば良いか。

回答1:フィッシャーの正確確率検定を使えば良いのではないか。

再質問:フィッシャーの正確確率検定は,データが少ないときと書いてあるが,多くても使って良いのか。
査読者には,マクネマー検定を使うべきではないかといわれているが,そうなのか?

回答2:データが多いときでも使える(どんなときにも正確なんだから)。マクネマー検定は対応のあるデータの場合なので,今のデータは対応がないんだから,マクネマー検定なんかつかえない。
ところで,被検者は一人なのか?

それに対する回答:一人だ。

回答3:被検者一人じゃ統計学で扱うには不適切。サンプルサイズを決めることもできない。検定ではなく,何が何パーセント増えた(減った)で十分なんじゃないか?

それに対する回答:一人だけど,何らかの方向性をいうことはできるだろう。
でも,回答ありがとう。これを参考に査読者に回答,再投稿する。

そして,きょうの夕方,発端の投稿が消去されたので,それに端を発する質疑応答が,全部消えた。

このデータにマクネマー検定を提案するようなヘボ査読者は,この実験デザインの問題点には全く気づかないということなんだろうか。なんという学会だろうね?


偏F値によるステップワイズ変数選択の参考資料をご教示願います。
  投稿者:米谷 2016/06/23(Thu) 17:12 No. 22040
お世話になっております。
以下のウェブサイトにある偏F値によるステップワイズ変数選択について質問があり、投稿をさせていただきました。
http://aoki2.si.gunma-u.ac.jp/R/sreg.html

こちらにある関数のうち、fmax、fmin、step.outの各関数内で使われている理論式を説明がある参考資料がございましたらご教示いただきたく存じます。

特にfmax、fminのsweep処理や、step.out関数のF値や偏回帰係数の標準偏差等について、参考資料がございましたら、そちらから内容を理解したく、ご教示いただければ幸いです。どうぞよろしくお願いいたします。

米谷

Re: 偏F値によるステップワイズ変数選択の参考資料をご教示願います。
  投稿者:青木繁伸 2016/06/23(Thu) 18:13 No. 22042
古い本になりますが,奥野忠一他「正・続 多変量解析法」日科技連,に書かれていたと思います。

Re: 偏F値によるステップワイズ変数選択の参考資料をご教示願います。
  投稿者:米谷 2016/06/23(Thu) 19:19 No. 22044
青木先生

お忙しいところご教示いただき誠にありがとうござました。

偏F値を使った変数選択の理論的な流れを紹介する資料がなかなか見つからず困っておりました。その本で確認させていただきます。

米谷


Shirley-Williams検定
  投稿者:赤岳 2016/06/23(Thu) 16:38 No. 22039
いつもご指導ありがとうございます。
さて、Shirley-Williams検定をRで実行するスクリプトを青木先生のサイトで見つけました。
source("http://aoki2.si.gunma-u.ac.jp/R/src/Shirley-Williams.R", encoding="euc-jp")
出力結果をみると
> Shirley.Williams(data, group, method="up")
4 3 2
2.189899 1.091791 1.243468
となっていますが、ここに表示されている2.189899等は対照群を1としたときの大きさだと思われますが、統計量やP値は、出力されないのでしょうか。
どうかP値の出し方を教えてください。よろしくお願いします。

Re: Shirley-Williams検定
  投稿者:青木繁伸 2016/06/23(Thu) 18:11 No. 22041
ページに書いてあるように
永田靖・吉田道弘(1997)「統計的多重比較法の基礎」サイエンティスト社
を参考にプログラムを書いたのですが,P 値の求め方は,書いてなかったと思います。
付表を参照だったかな?
出力される数値が何を表すかを含め,前述の本を見てください。

Re: Shirley-Williams検定
  投稿者:赤岳 2016/06/23(Thu) 18:33 No. 22043
青木先生
早速お返事いただきありがとうございました。
参考本を確認してみます。


非正規分布の検定について
  投稿者:下山 2016/06/18(Sat) 18:43 No. 22030
度々お世話になります。

ある要素Aの有無に対してスコアB(%)、スコアC(%)というN=150程度の表が有ります。
A B C
有 10% 17%
無 15% 18%
有 20% 28%
・・・

A-B、A-Cについてそれぞれ要素Aの有無とスコアの間に差が有るか、有るとすればどちらが大きいかを検定したいと考えています。B、Cの分布は正規性検定ではどちらも非正規という結果でしたので、U検定をしようかと思いました。U検定ではBもCもP<0.01でした。分散分析を行ったところA-Bは有無の2群間で分散は等しくなく、A-Cは等分散という結果になりました。

A-Cは中央値で評価しようかと思いましたが、A-Bについては有群と無群どちらがスコアが高いかをどう評価しようか、と壁に当たりました。
どのような方法がありますでしょうか?中央値を用いても良いのでしょうか?

もし正規性の検定以降の手順についても、間違っているところがありましたらご指摘いただけましたら大変ありがたく存じます。どうぞよろしくお願いいたします。

Re: 非正規分布の検定について
  投稿者:青木繁伸 2016/06/18(Sat) 20:32 No. 22031
そもそも,

> A B C
> 有 10% 17%
> 無 15% 18%
> 有 20% 28%
> ・・・・

というデータの見方がわかりません。

それと,検定に必要なのは,% ではなく,実数です。
本当のデータ(実数)が明らかになって困るとは思えないので,それぞれが何を表すかは言う必要はないのだから,ちゃんと説明する方がよいと思いますけど?
じゃないと,返答のしようがない。

Re: 非正規分布の検定について
  投稿者:下山 2016/06/19(Sun) 08:38 No. 22032
青木先生、失礼いたしました。
BとCはそれぞれ異なる計算式から算出した予測死亡率で、どちらもいくつかの調査項目から計算しています。ですので%単位になっています。

取り急ぎお返事させていただきました。

Re: 非正規分布の検定について
  投稿者:青木繁伸 2016/06/19(Sun) 17:49 No. 22033
検定できないでしょう

Re: 非正規分布の検定について
  投稿者:下山 2016/06/19(Sun) 20:39 No. 22034
青木先生 ご回答ありがとうございます。
このデータに関しては検定できない、ということで分析から外す方向で考えます。

後学のために以下のような場合にどう評価していくかも教えていただけましたらありがたく存じます。

N=150
A=某塾に通っているかどうか
B=期末試験の点数(通塾の有無で分けた場合等分散ではない)
C=統一模試の点数(通塾の有無で分けた場合等分散である)

として、点数がどちらも正規分布していない場合、通塾の有無で期末試験の点数に差があるか、あるとすればどちらが点数が高い傾向があるのか、および通塾の有無で統一模試の点数に差があるか、あるとすればどちらが点数が高い傾向があるのかを評価したい。この場合であればデータは実数となるかと思いますが、ご教授いただけましたらありがたく存じます。
大変恐縮でございますが、どうぞよろしくお願いいたします。

Re: 非正規分布の検定について
  投稿者:青木繁伸 2016/06/19(Sun) 21:04 No. 22035
この例と,前のと,どのような関係があるのですか?
本当は関係がないのに,別の問題だとして回答したことを(つまり,新たに提示されたデータについてはこれこれの検定が適切でしょうということを),あなたが勝手に(状況は全く違う)この問題にも適用するかもしれないという懸念があるのに,回答できないでしょう?

抽象的なことばかり言っていないで,データの性質も含めて記述しないと,なんの回答もできません(私のせいにされるのは心外で,そんなことこわくて,できません)。

Re: 非正規分布の検定について
  投稿者:下山 2016/06/20(Mon) 14:06 No. 22036
青木先生

検定を根本的に理解しておらず多々考え違いをしていることを痛感いたしました。ご気分を害してしまいましたことお詫びいたします。

理解が乏しい中で自分なりに解釈すると、正規分布が仮定できる場合にしろ分割表にしろ非正規分布の場合にしろ、実際のデータを吟味しないといけなくて、初心者向けの統計本に記載されているような、パラメトリックなら対応の有無を考慮してt検定やANOVAで、分割表ならχ二乗で、ノンパラメトリックなら・・・で、というように馬鹿の一つ覚えで1対1に決めてできるものではないということかと考えました。

必要十分なデータの性質についての情報を記載できる自信はありませんが(データの性質といわれても自分には分布や分散、外れ値くらいしか思いつかない程度の理解で)、添付の画像のようなデータです。Logitはそれぞれ異なる論文に記載されていた回帰式より計算しております。

はじめに先生がご回答されたように、こういうデータは検定できないのでしょうか(あるいは検定を行う意味・意義がない、検定すべきではない)。


Re: 非正規分布の検定について
  投稿者:鈴木康弘 2016/06/22(Wed) 07:15 No. 22037
 せっかく具体的なデータを出していただいたのに、一般論ですが、
 正規分布し分散も等しい2変数を利用して、2群に差があるかを検定するなら、多変量分散分析を使うはずです。

 でもこの条件に合わなければこれは使えない。
 だから2変数を一度に使うことはあきらめざるをえない(と思う)。

 正規分布もしてないし分散も異なる2群を単変量で検定するのだったら、Brunner-Munzel検定というのがあるようです。

Re: 非正規分布の検定について
  投稿者:下山 2016/06/22(Wed) 20:12 No. 22038
鈴木先生

御助言ありがとうございます。Brunner-Munzel検定は初めて聞きました。
調べてみたいと思います。

[1] [2]

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