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

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


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

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


複合アウトカム
  投稿者:あおいひげ 2018/10/17(Wed) 17:43 No. 22623
青木先生

ご回答有り難うございます。
正準相関分析・・・参考になりました。
Rを使用するのはかなりハードルが高いと思っています。
再度、研究デザインから検討し直します。

このスピードで、知りたいことに端的にお返事頂いていることに大変感銘を受けました。今後も、お世話になるかと思います。宜しくお願いします。

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/17(Wed) 22:39 No. 22624
今度は,スレッドを維持するように


複合アウトカム
  投稿者:あおいひげ 2018/10/16(Tue) 16:30 No. 22621
青木先生
回答有り難うございます。

複合アウトカムについては、死亡と再入院など、意味合いが異なるイベントを組み合わせることの意義について見直されてきているようです。
 現実的には、列挙したアウトカムの中でも1番重要視しているものをプライマリーアウトカムとし、それ以外はセカンダリーアウトカムとしてオマケぐらいの位置づけで解析をやるしかないのかもしれませんね。
 因みに、正準相関分析は私が使用しているEZRでは処理できそうにないのですが、他のソフトなら可能なものなのでしょうか。
 

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/16(Tue) 17:55 No. 22622
EZR は,R のシュガーコートなので,R の全てに対応しているわけではない
(EZR にどっぷりというのは避けた方がよい)

R には,cancor という関数がある(stats に含まれているので,library とかしないで,すぐ使える)
使用例(example)も見ることができる
> example(cancor)

cancor> ## No test:
cancor> ## signs of results are random
cancor> pop <- LifeCycleSavings[, 2:3]

cancor> oec <- LifeCycleSavings[, -(2:3)]

cancor> cancor(pop, oec)
$cor
[1] 0.8247966 0.3652762

$xcoef
[,1] [,2]
pop15 -0.009110856 -0.03622206
pop75 0.048647514 -0.26031158

$ycoef
[,1] [,2] [,3]
sr 0.0084710221 3.337936e-02 -5.157130e-03
dpi 0.0001307398 -7.588232e-05 4.543705e-06

以下略


複合アウトカム
  投稿者:あおいひげ 2018/10/15(Mon) 18:15 No. 22619
初めて投稿致します。

医学系臨床研究のデザインを検討しています。統計初心者で、いろいろ試行錯誤で苦しんでおりますのでご教授お願い致します。統計ソフトはEZRを使用しております。

集団をある値が高い群(H群)と低い群(L群)の2群に分けました。100〜200例 VS 300〜500例 (集計を開始していないので見積もり概算です) 。
 アウトカムとして、連続変数と二値変数の2種類の変数を全部で10個くらい想定しています。交絡因子は年齢、性別、既往歴、など10種類くらいを検討しています。
 多変量解析になるかと思いますが、複合アウトカムの処理方法が分かりません。アウトカムを一つずつ、重回帰やロジスティックなどを多重性に気を付けながらやるしかないのでしょうか。

曖昧な質問ですみません。
お教え下さいませ。

Re: 複合アウトカム
  投稿者:青木繁伸 2018/10/15(Mon) 19:21 No. 22620
結果が複数の変数ということならば,正準相関分析(重回帰分析の一般形)ということになるかな?実際に使った研究例は多くはないかも。


【R】 リスト構造からの読み取り
  投稿者:明石 2018/10/02(Tue) 21:46 No. 22616
青木先生 様;

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

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

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

ベイジアンネット bnlearn パッケージの出力結果のリスト構造からグラフ構造を取得したいと思います。
以下、グラフ描画までのコード例をお示しします。

library(bnlearn) # ベイジアンネットのパッケージ

df <- learning.test # 組み込みのデータセット

res <- hc(df)   # 構造学習 hill-climbing法を適用
plot(res, main = "hill-climbing法による構造学習") # 添付の画像ファイル

str(res)

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

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

所望する出力です。
A,B
A,D
B,E
C,D
F,E

--------------------------------
以下、str(res) の出力結果に、 銑┐涼躰瓩鯢嬪燭靴泙靴拭

List of 3
$ learning:List of 8
..$ whitelist: NULL
..$ blacklist: NULL
..$ test : chr "bic"
..$ ntests : num 40
..$ algo : chr "hc"
..$ args :List of 1
.. ..$ k: num 4.26
..$ optimized: logi TRUE
..$ illegal : NULL
$ nodes :List of 6
..$ A:List of 4             ⇒
.. ..$ mb : chr [1:3] "B" "C" "D"
.. ..$ nbr : chr [1:2] "B" "D"
.. ..$ parents : chr(0)
.. ..$ children: chr [1:2] "B" "D"   ⇒
..$ B:List of 4             ⇒
.. ..$ mb : chr [1:3] "A" "E" "F"
.. ..$ nbr : chr [1:2] "A" "E"
.. ..$ parents : chr "A"
.. ..$ children: chr "E"        ⇒
..$ C:List of 4              ⇒ァ      
.. ..$ mb : chr [1:2] "A" "D"
.. ..$ nbr : chr "D"
.. ..$ parents : chr(0)
.. ..$ children: chr "D"        ⇒
..$ D:List of 4
.. ..$ mb : chr [1:2] "A" "C"
.. ..$ nbr : chr [1:2] "A" "C"
.. ..$ parents : chr [1:2] "A" "C"
.. ..$ children: chr(0)
..$ E:List of 4
.. ..$ mb : chr [1:2] "B" "F"
.. ..$ nbr : chr [1:2] "B" "F"
.. ..$ parents : chr [1:2] "B" "F"
.. ..$ children: chr(0)
..$ F:List of 4              ⇒
.. ..$ mb : chr [1:2] "B" "E"
.. ..$ nbr : chr "E"
.. ..$ parents : chr(0)
.. ..$ children: chr "E"        ⇒
$ arcs : chr [1:5, 1:2] "A" "A" "C" "B" ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:2] "from" "to"
- attr(*, "class")= chr "bn"

何卒どうぞ、よろしくお願いいたします。
いつもありがとうございます。


Re: 【R】 リスト構造からの読み取り
  投稿者:青木繁伸 2018/10/02(Tue) 23:14 No. 22617
この件については,何も知りませんが,
以下のようにすれば,簡単に解が得られのでは?
> res$arcs
from to
[1,] "A" "B"
[2,] "A" "D"
[3,] "C" "D"
[4,] "B" "E"
[5,] "F" "E"
See bn-class for details.

arcs: the arcs of the Bayesian network (a two-column matrix, whose columns are labeled from and to). Undirected arcs are stored as two directed arcs with
opposite directions between the corresponding incident nodes.

Re: 【R】 リスト構造からの読み取り
  投稿者:明石 2018/10/03(Wed) 19:59 No. 22618
青木先生 様;

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

今回も大変に良い勉強をさせていただきました。
心から御礼を申し上げます。


ゼロでの割り算
  投稿者:青木繁伸 2018/09/09(Sun) 11:10 No. 22612
中澤さんが,2x2 より大きい行列で R の mcnemar.test の結果が NaN になる場合があることについて,書いている。

http://minato.sip21c.org/bulbul2/20180908.html
http://minato.sip21c.org/bulbul2/20160913.html
http://minato.sip21c.org/bulbul2/20160911.html

「対角成分の和がゼロの組合せを計算からスキップできる拡張マクネマー検定の関数」として
mct <- function(x) {
L <- NROW(x)
X <- 0
N <- 0
for (i in 1:(L - 1)) {
for (j in (i + 1):L) {
s <- (x[i, j] + x[j, i])
if (s > 0) {
N <- N + 1
X <- X + (x[i, j] - x[j, i]) ^ 2 / s
}
}
}
return(c( X2 = X, df = N, p = (1 - pchisq(X, N))))
}
を挙げている。
しかし,「このやり方が正しいという理論的根拠を見つけられないので悩ましいところだが」とあるのだが。

「N23とN32がゼロ,という状態だと,(N23+N32)が分母になる項で"Division by zero"エラーがでてしまって,カイ二乗値の計算結果自体がNaNになってしまう」そこで「分母がゼロになる組合せをスキップする」ということを書いている。

この記述は正しくない。N23=N32=0のとき,つまり(0-0)/(0+0) = 0/0 の計算結果は,NaN である("Division by zero"ではない。Python だと,ZeroDivisionError: division by zero になる)。何かに NaN を足すと NaN になり,それ以降 NaN に何を足しても NaN のまま。ということで,R の mcnemar.test は NaN を返すと,中澤さんも書いている。

N23=N32=0 のときの計算結果が NaN になった場合は足し算をスキップすればよいわけで,N23=N32=0 の場合は計算をスキップしても同じではある。またその場合には自由度のカウントはしない。

以下のような関数を書いてみた。
mct2 <- function(x) {
i <- col(x)[lower.tri(x)]
j <- row(x)[lower.tri(x)]
a <- mapply(function(i, j) (x[i, j] - x[j, i]) ^ 2 / (x[i, j] + x[j, i]), i, j)
chisq <- sum(a, na.rm = TRUE) # NaN も除かれる
df <- sum(!is.nan(a)) # !is.na(a) でもよい
return(c(X2 = chisq, df = df, p.value = pchisq(chisq, df, lower.tail = FALSE)))
}
「ゼロ割」は避けるのではなく,積極的に正しく対処するべきということか?(同じようなことではあるが)

蛇足ではあるが,自由度を減ずるのは,カイ二乗分布の再生性から明らか(自由度 1 のカイ二乗分布にしたがう確率変数 k 個の和は自由度 k のカイ二乗分布にしたがう)。

R の mcnemar.test を修正するとすれば,
STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)])

STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
PARAMETER <- sum(!is.nan(y[upper.tri(x)]^2/x[upper.tri(x)]))
にする(前の方の PARAMETER <- r * (r - 1)/2 は不要)。
プログラムにおいて,特殊な条件の発生を見逃すと命取りというところ。
しかしさすがに,sum(y[upper.tri(x)]^2/x[upper.tri(x)]) とは,にくいねえ。
> x = matrix(c(c(7, 2, 5, 3, 11, 2, 0, 3, 4, 0, 5, 9, 3, 15, 14, 10)), 4, 4)
> mct(x)
X2 df p
15.428836864 5.000000000 0.008678852
> mct2(x)
X2 df p.value
15.428836864 5.000000000 0.008678852

> mcnemar.test(x) # 元のまま

McNemar's Chi-squared test

data: x
McNemar's chi-squared = NaN, df = 6, p-value = NA

> mcnemar.test2(x) # 修正したもの

McNemar's Chi-squared test

data: x
McNemar's chi-squared = 15.429, df = 5, p-value = 0.008679
なお,「"Division by zero"を避けるために全部のセルに0.01を足す」というのは明らかに間違い(自由度が減らない)。

Re: ゼロでの割り算
  投稿者:青木繁伸 2018/09/09(Sun) 12:49 No. 22613
2 0
0 5

というのと

2 1
1 5

というのは,0-0 = 1-1 = 0 という情報を与えるという意味では同じなので,自由度は減らさなくてもよいか

STATISTIC <- sum(y[upper.tri(x)]^2/x[upper.tri(x)], na.rm=TRUE)
だけの変更で,これは実質的には NaN の場合に 0 を加えるというのと同じであるから
PARAMETER はそのまま PARAMETER <- r * (r - 1)/2

とすれば,mcnemar.test(x+0.01) もありか。ただし,加える数は他のまともなセルのことも考えるとできるだけ小さい方が良いだろうとは思う。例えば +1e-15 とか

Re: ゼロでの割り算
  投稿者:中澤 2018/09/09(Sun) 21:30 No. 22614
記事ありがとうございます。
確かにNaNが正しいですね。不正確な記述でした。
自由度の扱いは悩ましいですね。スキップでなくて分子がゼロなのでゼロとみなすということですね。勉強になりました。

Re: ゼロでの割り算
  投稿者:中澤 2018/09/19(Wed) 17:55 No. 22615
http://minato.sip21c.org/bulbul2/20180915.html
に書きましたが,カテゴリ数が3以上の場合の拡張マクネマー検定と同様の帰無仮説を検定する方法としてBhapkarの検定というのがあって,irrパッケージのbhapkar()で実行できました。いま書いている論文ではこれを使おうと思います。


パネルデータの欠損値補完について
  投稿者:羽島 2018/09/04(Tue) 09:21 No. 22609
青木先生

羽島と申します、この度はお世話になります。

Rを使って以下のようなデータを分析しており、その中で住所コードの欠損値を補完する必要が出てきました。LOCFのような手法での補完では精度に不安があるため、個人IDを参照することで欠損値を補完したいのですが、どのようにプログラムを書けば良いでしょうか?

お忙しいところ大変恐縮ですが、ご教示頂けましたら幸いです。

個人ID  年   住所コード

100   2015    810
100   2016    NA
100   2017    810
201   2015    NA
201   2016    NA
201   2017    666

Re: パネルデータの欠損値補完について
  投稿者:青木繁伸 2018/09/04(Tue) 13:03 No. 22610
あまり簡単なプログラムは書けないかな
d <- data.frame(ID = c(100, 100, 100, 201, 201, 201, 330, 330, 330),
year = c(2015, 2016, 1017, 2015, 2016, 2017, 2015, 2016, 2017),
code = c(810, NA, 810, NA, NA, 666, NA, 543, NA))
tbl <- unique(na.omit(d[c("ID", "code")]))
for (i in seq_len(nrow(d))) {
if (is.na(d$code[i])) {
d$code[i] = tbl[which(d$ID[i] == tbl[,1]), 2]
}
}
d

Re: パネルデータの欠損値補完について
  投稿者:羽島 2018/09/04(Tue) 18:15 No. 22611
青木先生

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

元のデータは3万件ほどの個人データがあり手作業での修正を覚悟しておりましたが、
先生のお力で無事補完することが出来ました。

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


【R】 二重繰り返し処理の高速化
  投稿者:明石 2018/09/01(Sat) 10:42 No. 22606
青木先生 様;

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

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

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

記事 No.22590 「繰り返し処理の高速化」では、ベクトル化を用いた、劇的な高速化のご教示をいただきました。
ありがとうございました。

今回は、二重繰り返し処理です。

例示データは、先生のサイトからお借りいたします。

http://aoki2.si.gunma-u.ac.jp/R/my-cor.html
単相関係数,偏相関係数,重相関係数を計算する

x <- matrix(c( # 5ケース,4変数のデータ行列例(ファイルから読んでも良い)
1, 5, 6, 4,
2, 14, 5, 3,
3, 3, 4, 2,
4, 2, 6, 6,
3, 4, 3, 5
), ncol=4, byrow=TRUE)

ここでは、相関係数行列を計算していますが、
今回、私がやりたいのは、ケース(行ベクトル)同士のコサイン類似度行列です。

以下、プログラムを作成しました。

x <- matrix(c( # 5ケース,4変数のデータ行列例(ファイルから読んでも良い)
1, 5, 6, 4,
2, 14, 5, 3,
3, 3, 4, 2,
4, 2, 6, 6,
3, 4, 3, 5
), ncol=4, byrow=TRUE)

nr <- nrow(x)
mx <- matrix(0,ncol=nr,nrow=nr)
diag(mx) <- 1

for(i in 1:(nr-1)) {
for(j in (i+1):nr) {
s <- (x[i,] %*% x[j,]) / (sqrt(sum(x[i,]^2)) * sqrt(sum(x[j,]^2)))
mx[i,j] <- s
mx[j,i] <- s
}
}

> mx
[,1] [,2] [,3] [,4] [,5]
[1,] 1.0000000 0.8438196 0.9183979 0.8735555 0.8992005
[2,] 0.8438196 1.0000000 0.7847512 0.5725026 0.7829858
[3,] 0.9183979 0.7847512 1.0000000 0.9132886 0.9081355
[4,] 0.8735555 0.5725026 0.9132886 1.0000000 0.9229730
[5,] 0.8992005 0.7829858 0.9081355 0.9229730 1.0000000

私が扱いたいデータでは、
ケース数は、数千
属性数は、数百
ですので、高速化処理ができれば大変にありがたいです。

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

Re: 【R】 二重繰り返し処理の高速化
  投稿者:青木繁伸 2018/09/02(Sun) 11:31 No. 22607
二重 for-loop の内側で使われる
x[i,] %*% x[j,]
は,行ベクトルの積和ですから
for-loop なしの
z <- x%*%t(x)
です。

同じく,分母は
sqrt(sum(x[i,]^2))という積和の平方根ですから x%*%t(x) の対角要素の平方根です
d <- sqrt(diag(z))
二重 for-loop で使われる
sqrt(sum(x[i,]^2)) * sqrt(sum(x[j,]^2))
は,
d2 <- outer(d, d)
です。
以上をまとめて,ループなしで答えが求まります。
z <- x%*%t(x)
d <- sqrt(diag(z))
d2 <- outer(d, d)
mx2 <- z / d2
5000x500の行列だと,私のマシンではfor-loop プログラムは約400秒,ベクトル演算を使ったプログラムは約10秒で,およそ40倍速くなります

Special Thanks !(Re: 【R】 二重繰り返し処理の高速化)
  投稿者:明石 2018/09/02(Sun) 12:24 No. 22608
青木先生 様;

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

劇的な改善方法をご教示をいただき、誠にありがとうございます。

表層的ではなく、計算の本質を見抜くことの重要性を教えていただきました。

心から、心より、御礼を申し上げます。
本当にありがとうございました。


適合度のexact testの結果の書き方について
  投稿者:黒岩 2018/08/24(Fri) 19:13 No. 22603
青木先生

いつもお世話になっております。
下記についてご教示頂けると幸いです。

gftで適合度のexact testを行うと、
カイ二乗値、自由度、P値、正確なP値、が出力されますが、
レポートなどでこの結果を書くときの
適切な記述方法を教えて頂けないでしょうか?

単純にこの4つを書き並べれば良いのでしょうか?
P値をはっきり区別するために、
例えば漸近有意確率P、正確有意確率P、などと
書き分ける必要はないのでしょうか?
あるいは、正確なP値のみではいけないのでしょうか?

Re: 適合度のexact testの結果の書き方について
  投稿者:青木繁伸 2018/08/24(Fri) 21:39 No. 22604
「適合度の正確検定による有意確率 p < 0.xxx」でよいと思いますが...

Re: 適合度のexact testの結果の書き方について
  投稿者:黒岩 2018/08/25(Sat) 04:32 No. 22605
青木先生

早速の回答、ありがとうございました!


NaNを1に置き換えたい
  投稿者:赤岳 2018/08/21(Tue) 17:30 No. 22600
いつもお世話になっております。
次の文をRで指定しましたが、logの値がマイナスになり、エラーとなります。
logの中の乱数がマイナスやゼロの場合、1を返すようにしたいのですが、うまく文をかけません。
ご教示いただけないでしょうか。
*****************
> N <-100
> P<- rnorm(N, 5, 1)
> dP <- c(1:N)*0
> EP <- c(1:N)*0
> count<- 0
> for (i in 1:N) {
+ dP[i] <- P[i] - log(rnorm (1, 1, 0.5))
+ EP[i] <- exp (dP[i])/exp(P[i])
+ if (EP[i] >0.5)
+ count <- count+1
+ }
log(rnorm(1, 1, 0.5)) で警告がありました: 計算結果が NaN になりました
if (EP[i] > 0.5) count <- count + 1 でエラー:
TRUE/FALSE が必要なところが欠損値です
> count/N
[1] 0.41

Re: NaNを1に置き換えたい
  投稿者:青木繁伸 2018/08/21(Tue) 22:23 No. 22601
本筋ではないところをコメント
dP <- c(1:N)*0 は dp <- numeric(N) のほうが適切でしょう。

rnorm(1, 1, 0.5) が 0 以下のとき,log(rnorm(1, 1, 0.5)) を 1 にするのですか?

for ループを使っている場合は,特に何の工夫も必要ないです
N <-100
P <- rnorm(N, 5, 1)
dP <- c(1:N)*0
EP <- c(1:N)*0
count<- 0
for (i in 1:N) {
x <- rnorm (1, 1, 0.5) # ここから
if (x <= 0) {
x <- 1
} else {
x <- log(x)
}
dP[i] <- P[i] - x # ここまで追加(変更)
EP[i] <- exp (dP[i])/exp(P[i])
if (EP[i] >0.5)
count <- count+1
}
count/N


for ループをなくしてみましょう。ベクトル演算します。
N <- 100
p <- rnorm(N, 5, 1)
x <- rnorm(N, 1, 0.5) # まとめて N 個生成
x <- ifelse(x <= 0, exp(1), x) # x <= 0 のとき log(x) を 1 にするなら,x = exp(1) にしておく
dP <- P - log(x) # x が正なら log(x),正でなければ log(exp(1)) = 1 になる
EP <- exp(dP) / exp(P) # EP の計算
mean(EP > 0.5) # count = sum(EP > 0.5); count/N は,mean(EP > 0.5)


注: ifelse(x <= 0, 1, log(x)) としても,一度全ての x に log(x) を適用して,x <= 0 のときに 1,そうでないときに log(x) の値を取るように動作するので,エラーが出る。
そこで,上のようにへんな(?)ことをする。

Re: NaNを1に置き換えたい
  投稿者:赤岳 2018/08/21(Tue) 22:53 No. 22602
青木先生
ご教示ありがとうございました。
> mean(EP > 0.5) # count = sum(EP > 0.5); count/N は,mean(EP > 0.5)
このmeanがこんな風に使われるのは初めて知りました。
また、確かにifelse(x <= 0, 1, log(x)) のようにしましたが、エラーメッセージはクリアーできませんでした。
解決しました。ありがとうございました。
(p <- rnorm(N, 5, 1)は、P <- rnorm(N, 5, 1)と修正しました。)


直接法と群間差の検定
  投稿者:M 2018/08/18(Sat) 21:48 No. 22593 HomePage
青木先生  皆様

お世話になります。いろいろ調べたのですが、自信がないので、質問させていただきたいと存じます。

添付のサイトのP182 表88のような検定を行いたいのですが、検定方法がわかりません。

質問1: 直接法とは何でしょうか?

表88には、注釈に、注1)年齢階級(20−29歳,30−39歳,40−49歳,50−59歳,60−69歳,70歳以上の6区分)と世帯員数(1人,2人,3人以上世帯の3区分)での調整値。割合に関する項目は直接法,平均 値に関する項目は共分散分析を用いて算出。 む

と書いてあります。

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

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

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

群間比較について,年齢調整を行い分析を行ったものについては,割合に関す る項目はCochran-Mantel-Haenszel検定,平均値に関する項目は,共分散分析を 行った。

とあります。

お手数をお掛け致しますが、どうぞよろしくお願い致します。

Re: 直接法と群間差の検定
  投稿者:青木繁伸 2018/08/19(Sun) 06:12 No. 22594
サイト情報がありません(家マークでのリンクではなく,本文中に記載してください。)

というわけで,詳細は分かりませんが,直接法については以下のごとく。

直接法とは年齢調整死亡率の計算に用いられる方法。ほかに,間接法もある。
死亡率に限らず,率の調整に使われることもある。
直接法についての説明は,
http://www.takenet.or.jp/~hayakawa/u-tanqa04.htm
http://www.kenko.pref.yamaguchi.lg.jp/img/kenko21/download/h26map/h26kaisetsu-new.pdf
など,
標準化死亡比の検定は,
https://www.niph.go.jp/soshiki/07shougai/datakatsuyou/data/h22smr-kinki.pdf
などを参照のこと。

ここで質問しても詳細な(適切な)回答が得られるとは限りません。
著者に直接お問い合わせになるのが筋でしょう。

Re: 直接法と群間差の検定
  投稿者:M 2018/08/19(Sun) 14:43 No. 22596
青木先生 

お返事、どうもありがとうございました。
また、リンクが添付できず大変失礼いたしました。

リンクは、
https://www.mhlw.go.jp/bunya/kenkou/eiyou/dl/h26-houkoku.pdf

ページ 182 表88 です。

特に、群間差の検定を行いたいのです。この表88では、運動習慣のない者(%)と表示してあり、P値で群間の検定結果が表示されています。
どうやって検定したものか、知りたいのです。

表の注ではなく、報告書の冒頭に「割合に関す る項目はCochran-Mantel-Haenszel検定」
とありますので、Cochran-Mantel-Haenszel検定かなと思うのですが、そうするとすべて2値にしなければならないでしょうか?

例えば 運動習慣 ある・ない× 性別 男・女 

層として、収入200万円未満・600万円以上、
     200-600未満・600万円以上
     200万円未満・200-600未満

の3回、Cochran-Mantel-Haenszel検定を行い、有意確率で、層に設定した収入の違いに関係なく有意に差があるかを検定すればよろしいでしょうか?

T検定では、こういう組み合わせにした検定はご法度のはずで、躊躇しております。

うまく説明できたか不安なのですが、厚労省が回答してくれるかわかりませんので、再質問させていただきました。お手数をお掛け致しますが、どうぞよろしくお願い致します。

Re: 直接法と群間差の検定
  投稿者:青木繁伸 2018/08/19(Sun) 18:32 No. 22597
厚生省は発注元?で,実際は「立研究開発法人医薬基盤・健康・栄養研究所は,厚生労働省に提出された調査票について,入力・集計・作表を行った。」わけで,研究者レベルの質問には答えてくれるかも。
主要な成果については論文として投稿される(されている)かも。この年のデータのものでなくても,以前のものにも同じ分析手法が使われていると思うので,もう少し丁寧な説明や参考文献も記載されているかも。

さて,182ページ表88ですが,

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

「割合に関する項目は Cochran-Mantel-Haenszel検定」というのは,たとえば,50ページの表8の脚注に同じ記述がありますよね。層は年齢階級なのかな(よくわからない)。

統計解析過程についての詳細な説明がないので,よくわかりません。

Re: 直接法と群間差の検定
  投稿者:青木繁伸 2018/08/19(Sun) 20:54 No. 22598
解析過程をもっとも正確に記述できるのはプログラミング言語なので,全部とはいわず代表的なデータ解析を行ったプログラムを掲載するのは優れた試みだと思いますね。
R なんかは,最適でしょう。まさか Excel でやってますということはないと思いますけど。

Re: 直接法と群間差の検定
  投稿者:M 2018/08/20(Mon) 13:02 No. 22599
青木先生

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

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

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

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


【R】 繰り返し処理の高速化
  投稿者:明石 2018/08/14(Tue) 18:06 No. 22590
青木先生 様;

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

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

---------------------------
受験者の生年月日と試験日から、試験日での年齢を計算したいと思います。

受験者データはデータフレームdatです。
生年月日の変数名はbirthです。日付型です。

試験日の変数名はtestです。日付型です。
test <- as.Date("2018-9-13")

以下の繰り返し処理で、受験者の試験日での年齢yearを計算できることを確認しました。
for (i in 1:nrow(dat)) {
dat$year[i] <- length(seq(dat$birth[i], test, by = "year")) - 1
}

受験者データが大規模ですので、繰り返し処理を高速化したいと思います。

ご教示いただければ大変に助かります。
何卒どうぞ、よろしくお願いいたします。

Re: 【R】 繰り返し処理の高速化
  投稿者:青木繁伸 2018/08/15(Wed) 02:12 No. 22591
50000 人分のテストデータを作ってやってみました
test <- as.Date("2018-9-13")
dat <- data.frame(birth=as.Date(Sys.Date()-1:50000)) # 50000 人分のテストデータ

system.time({
for (i in 1:nrow(dat)) {
dat$year[i] <- length(seq(dat$birth[i], test, by = "year")) - 1
}
})
ユーザ システム 経過
11.711 6.689 18.254

年齢の計算は,簡単なのです。試験年から生まれ年を引くだけ。
ただし,「生まれ月が試験月より後(大きい)なら」,
または,「生まれ月が試験月と同じで,生まれ日が試験日より後(大きい)なら」
1を引く。

わざわざ as.Date となっているけど,そこから年月日(スカラー y, m, d)を取り出す。
試験年月日は yyyy, mm, dd の整数で与える(as.Date から取りだしてもよいけど)

age2 <- function(x, yyyy = 2018, mm = 09, dd = 13) {
y <- as.integer(substr(x, 1, 4))
m <- as.integer(substr(x, 6, 7))
d <- as.integer(substr(x, 9, 10))
val <- yyyy - y
if (m > mm || (m == mm && d > dd)) {
val <- val-1
}
val
}
system.time({
for (i in 1:nrow(dat)) {
dat$year2[i] <- age2(dat$birth[i])
}
})

でも,これでは for を含むし,年月日の取り出しコストが掛かるのでかえって遅い。

ユーザ システム 経過
13.662 6.667 20.312

一応答えは合っているようだ。

> all(dat$year == dat$year2)
[1] TRUE

for を避けるためにベクトル演算を行うようにする。
関数には dat$birth ベクトルを渡す。
年月日の取り出し,それに基づく年齢(val)の計算も全部がベクトル演算。

age3 <- function(birth, yyyy = 2018, mm = 09, dd = 13) {
y <- as.integer(substr(birth, 1, 4))
m <- as.integer(substr(birth, 6, 7))
d <- as.integer(substr(birth, 9, 10))
val <- yyyy - y
ifelse(m > mm | (m == mm & d > dd), val-1, val)
}
system.time(dat$year3 <- age3(dat$birth))

27 倍ほど速くなった。

ユーザ システム 経過
0.424 0.023 0.448

答えも合っているようだ。

> all(dat$year == dat$year3)
[1] TRUE

なお,val を計算している 2 行を(ifelse を使うのが気持ち悪いからといって)

yyyy - y - (m > mm | (m == mm & d > dd))

にしても,ほとんど速度に違いはない。

年月日が3つのベクトルで与えられれば as.Date から取り出す必要がないのでもっと速い。
birth を "2018-08-15" とわざわざ入力しているのか,2018, 8, 15 という3つの整数を入力しているのかによる。
後者なら,その3つの整数から文字列を作って as.Date していることになる。
3つの整数値を入力してそれをそのまま使うなら,人間もコンピュータも無駄なことをしなくてもよいということになる。

御礼(Re: 【R】 繰り返し処理の高速化)
  投稿者:明石 2018/08/15(Wed) 08:04 No. 22592
青木先生 様;

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

多面的にご教示いただき、
特に、for を避けるためにベクトル演算を行うようにするところは、
大変に勉強になります。

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


Rの代入演算子について
  投稿者:波音 2018/08/09(Thu) 12:18 No. 22587
素朴な疑問ですが、ここのところ青木先生が代入演算子を

<-

ではなく、

=

と書かれるようになったのには、何か意味があるのでしょうか?

※書名は忘れましたが = ではなく <- の方がいいと書いてあった記憶があります。

最近になってソースコード表現の定石も変わってきているのでしょうか、、、

Re: Rの代入演算子について
  投稿者:青木繁伸 2018/08/09(Thu) 15:00 No. 22588
ソースコード表現の定石は変わっていないと思います
「GoogleのRスタイルガイドは、割り当てのために「=」を禁止する」ということもあり,
? "<-" で表示される R のマニュアルでも,The operator <- can be used anywhere, whereas the operator = is only allowed at the top level (e.g., in the complete expression typed at the command prompt) or as one of the subexpressions in a braced list of expressions. とありますので,<- を使うのが妥当であることに変わりありません。

= を使うようになったのは些細な理由からですが,= を使うようになっての感想は,何の不便もなく,タイプ数も(わずか1文字の違いですが)少なく,他の言語を使ってきた経験からも = の方が違和感がないなあというだけです。

公の場では <- を使うようにしたいと思います。

Re: Rの代入演算子について
  投稿者:波音 2018/08/10(Fri) 12:01 No. 22589
回答ありがとうございます。

私はRが経験年数的には一番長いので <- に慣れてしまっていますが、多言語習得ユーザーとしては(明示的な問題が生じない限り) = でも不都合はなさそうと思います。

※代入演算子とは別に、最近のR使いの方はパイプを多用した記述をされるので、コードを読むのに一苦労することが多いです。


Dataを部分的に取り出す方法
  投稿者:赤岳 2018/08/07(Tue) 22:31 No. 22584
いつも勉強させていただいています。
Rの基本操作だと思いますが、わかりません。
今、次のようなデータが1000個あるとします。
x<- rnorm(1000, 15.5, 3.3)

このデータの95%区間にあるデータを取り出したいのですが、どのように取り出せばよいのでしょうか。
たとえば、
> y<- quantile(x, c(0.025, 0.975))
とすると、yのlength=2となり、950個のデータが取り出せません。
どなたか、yが950個のベクトルになるようなコマンドを教えてもらますでしょうか。

Re: Dataを部分的に取り出す方法
  投稿者:青木繁伸 2018/08/08(Wed) 05:48 No. 22585
x <- rnorm(1000, 15.5, 3.3) のベクトル x から,10 以上, 15 以下のデータを取り出すときはどのようにしますか?

x[x >= 10 & x <= 15] ですね

では y <- quantile(x, c(0.025, 0.975)) から,下限と上限を取り出すのはどのようにしますか?

y[1], y[2] ですね

それをあわせれば

x[x >= y[1] & x <= y[2]] です

Re: Dataを部分的に取り出す方法
  投稿者:赤岳 2018/08/08(Wed) 08:26 No. 22586
青木先生、
ありがとうございました。
ここ2日くらいそのような関数というかコマンドばかり探していました。
今後ともよろしくお願いします。


ノンパラ、パラメトリックの決め方
  投稿者:コロン 2018/07/27(Fri) 07:51 No. 22581
いつもお世話になっております。

基本的なことで恐縮なのですが、ご教示いただけますか。

手元に、テスト得点データがあります。実験群、統制群の2グループで、サンプルサイズはそれぞれ、30、25です。平均点の差の検定を行いたい場合は、普通はt検定を行うのではと考えます。

これを、サンプルサイズが小さいという理由で、ノンパラメトリック検定で行うことは大丈夫でしょうか。なぜこのような質問をしたかと言いますと、ある論文の査読をしているのですが、サンプルサイズが小さいとありますが、GPowerで必要な被験者数を求めると上のサイズは小さいとは言えず、テスト得点ですので、母分布は正規分布をなしていると考えます。また、本論文では、正規分布をなしていない根拠として、手元のデータに対して正規分布をなしているかどうかの検定を行い、正規分布をなしていないとも言っています。

そもそもとしまして、ノンパラにするかどうかは手元のデータが正規分布をなすかどうかで決めるものなのか、以前、この掲示板のどこかで指摘ぐあったような記憶があります。

どういう視点で、ノンパラか否かを決めるべきか、また、決定のプロセスの中でやってはいけないこと、ぜひご教示ください。

Re: ノンパラ、パラメトリックの決め方
  投稿者:青木繁伸 2018/07/28(Sat) 20:49 No. 22582
> サンプルサイズが小さいという理由で、ノンパラメトリック検定で行う

あまり説得力はないでしょう
しかし,

> GPowerで必要な被験者数を求めると上のサイズは小さいとは言えず

とはいっても,マン・ホイットニーの U 検定の検定力は,条件を満たしているときの平均値の差の検定(t 検定)の 95% ほどもあるので,劣っている検定手法を採用したと批難するのはちょっと厳しすぎるかも。

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

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

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

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

「データを見てから検定法を選択すること」はやってはいけないことと言えるでしょう。

Re: ノンパラ、パラメトリックの決め方
  投稿者:コロン 2018/07/29(Sun) 08:08 No. 22583
青木先生

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

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

失礼いたします。


R 半角文字→全角文字への変換
  投稿者:明石 2018/07/25(Wed) 19:23 No. 22578
青木先生 様;

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

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

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

半角、全角、大文字、小文字が混在する文字列を、全角大文字に統一したいのです。

例示ですが、

moji <- "ボーイング787 ボーイング787 boeing787 Boeing787 boeing787 BOEING787"

"ボーイング787 ボーイング787 BOEING787 BOEING787 BOEING787 BOEING787"

大文字に置換は、toupper(moji)を使うえばよいことは分かりましたが、
誠に情けないのですが、半角→全角の関数を見つけることができません。

ご教示いただけましたら大変に助かります。
お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

Re: R 半角文字→全角文字への変換
  投稿者:青木繁伸 2018/07/25(Wed) 19:53 No. 22579
文字の変換には基本的には chartr 関数を使えば良いです
例題の場合には,大文字変換に toupper を使えば,少し簡単に書けます
> moji = "ボーイング787 ボーイング787 boeing787 Boeing787 boeing787 BOEING787"
> moji = toupper(moji)
> chartr("A-Z0-9", "A-Z0-9", moji)
[1] "ボーイング787 ボーイング787 BOEING787 BOEING787 BOEING787 BOEING787"

> moji = "ボーイング787 ボーイング787 boeing787 Boeing787 boeing787 BOEING787"
> chartr("a-zA-Za-z0-9", "A-ZA-ZA-Z0-9", moji)
[1] "ボーイング787 ボーイング787 BOEING787 BOEING787 BOEING787 BOEING787"


Re: R 半角文字→全角文字への変換
  投稿者:明石 2018/07/26(Thu) 09:22 No. 22580
青木先生 様;

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

テキスト分析で、語のゆらぎを吸収する際に必要な前処理であり、大変に助かりました。

つまらない質問で、ご高名な青木先生に失礼にならないか、
小心者ですので、躊躇していましたが、
勇気をふりしぼってご相談して良かったと思います。

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


対応のある2標本の母分布の違いの分析について
  投稿者:HG 2018/06/21(Thu) 12:51 No. 22574
対応のある2標本の母分布に違いがあるかを検証したいのですが、
どのような分析手法を採用すべきでしょうか。

2標本コルモゴロフ・スミルノフ検定は独立2群の比較なので
この場合は使えないと考えております。

ご教示のほど、よろしくお願いいたします。

Re: 対応のある2標本の母分布の違いの分析について
  投稿者:青木繁伸 2018/06/21(Thu) 21:07 No. 22575
「母分布に違い」というのは,位置の母数の違いなのか,散布度の違いなのか?それとも,両方含めての違いなのか?

統計手法の選択ガイド
http://aoki2.si.gunma-u.ac.jp/FlowChart/Tutorial.html
をやってみるとどうなりました?

Re: 対応のある2標本の母分布の違いの分析について
  投稿者:HG 2018/06/22(Fri) 13:22 No. 22577
青木先生
コメントありがとうございます。

「母分布に違い」とは散布度(分布の形)の違いを意図しておりました。説明不足であり、
申し訳ありません。なお、データは離散型の比率尺度データです。

統計手法の選択ガイド
http://aoki2.si.gunma-u.ac.jp/FlowChart/Tutorial.html
では、「検定→その他」で適合度の検定のページ
http://aoki2.si.gunma-u.ac.jp/lecture/tekigoudo.html
にたどり着きますが、対応のある2群の散布度(分布の形)に違いがあるかを検証する分析が
見つけられず、質問させていただきました。

また、
http://aoki2.si.gunma-u.ac.jp/FlowChart/kentei_heikinti_taiouari.2.html
のページで紹介されている検定は平均値(代表値)の比較が目的であり、
散布度(分布の形)の違いを検証するものではないと思っており、判断に迷っております。


R 複数のNA判定に基づく代入
  投稿者:明石 2018/06/16(Sat) 19:06 No. 22571
青木先生 様;

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

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

添付ファイルでお示しをします。 ご覧ください。

表1から表2を作成したいと思います。

表1のデータフレームは、以下のようになっています。
・socreは、すべてNA
・score1とscore2の組み合わせは、以下の3ケースのいずれか
  ・score1がNA,score2は数値
  ・score1が数値,score2はNA
  ・score1,score2両方ともにNA

score1、score2で、数値がセットされている側の値をscoreにセットする。

データフレームを dat として、Rコードを書きました。

jun <- which(! is.na(dat$score1))
dat$score[jun] <- dat$score1[jun]

jun <- which(! is.na(dat$score2))
dat$score[jun] <- dat$score2[jun]

score1、score2を、各々、NAの判定をして、scoreに代入しています。

これらを1つにまとめる書き方があれば、ご教示をお願いいたします。

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


Re: R 複数のNA判定に基づく代入
  投稿者:青木繁伸 2018/06/16(Sat) 21:35 No. 22572
いろいろやり方はあるでしょうが,もっとも真っ正直に書くならば
> dat = data.frame(
+ score=NA,
+ score1=c(NA, NA, 2.174, 6.042, NA, 2.656, NA, NA, NA, NA, 3.153, 4.344),
+ score2=c(1.481, 4.113, NA, NA, 2.674, NA, 2.500, 2.276, 4.162, NA, NA, NA))

> dat$score = ifelse(is.na(dat$score1) & !is.na(dat$score2), dat$score2, ifelse(!is.na(dat$score1) & is.na(dat$score2), dat$score1, NA))
> dat
score score1 score2
1 1.481 NA 1.481
2 4.113 NA 4.113
3 2.174 2.174 NA
4 6.042 6.042 NA
5 2.674 NA 2.674
6 2.656 2.656 NA
7 2.500 NA 2.500
8 2.276 NA 2.276
9 4.162 NA 4.162
10 NA NA NA
11 3.153 3.153 NA
12 4.344 4.344 NA
トリッキーな書き方ならば,出現する可能性のある数値よりも小さい値を dummy にわりあてて,以下のようにすることも可能
> dummy = -99999 # 出現する可能性のある数値が正ならばこのような値を仮定する
> dat[is.na(dat)] = dummy
> dat$score = ifelse(dat$score1 > dat$score2, dat$score1, dat$score2)
> dat[dat == dummy] = NA
> dat
score score1 score2
1 1.481 NA 1.481
2 4.113 NA 4.113
3 2.174 2.174 NA
4 6.042 6.042 NA
5 2.674 NA 2.674
6 2.656 2.656 NA
7 2.500 NA 2.500
8 2.276 NA 2.276
9 4.162 NA 4.162
10 NA NA NA
11 3.153 3.153 NA
12 4.344 4.344 NA
どちらが効率がよいか?

Re: R 複数のNA判定に基づく代入
  投稿者:明石 2018/06/16(Sat) 21:54 No. 22573
青木先生 様;

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

今回も、ご丁寧なご教示をいただき、大変に良い勉強をさせていただきました。

2つの方法とも理解できましたが、ifelse()関数を使わせていただきます。

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


サンプリング誤差について
  投稿者:統計迷い人 2018/06/07(Thu) 11:20 No. 22563
サンプリング誤差について質問です。

サンプリング誤差について、視聴率での以下のサイトでの考え方や計算方法は分かったのですが、回答選択肢が3以上ある場合、例えば、非常に好き、やや好き、どちらでもない、やや嫌い、非常に嫌い、などでも、誤差について同様に考えて良いのでしょうか?

https://toukeigaku-jouhou.info/2015/08/25/audience-rate/


よろしくご教示のほど、お願いいたします。

Re: サンプリング誤差について
  投稿者:鈴木康弘 2018/06/10(Sun) 09:53 No. 22564
 サンプリング誤差って何かな?と思ったら、母集団の比率の95%信頼区間の求め方のことですか。
 選択肢がたくさんあっても同じだと、僕は思います..?

Re: サンプリング誤差について
  投稿者:統計迷い人 2018/06/10(Sun) 11:55 No. 22565 HomePage
鈴木康弘 先生

ご回答、ありがとうございます。

視聴率や知名度、購入経験のように、yes no いずれかから一つ選ぶ場合は二項分布であって、それが正規分布に近似することから、95%の信頼区間推定が出来るとテキストにあったのですが、

非常に好き、やや好き、普通、やや嫌い、非常に嫌いとかの、
多項分布の場合、誤差をどう考えたらよいかご教示下さい。

例えば、サンプル数600、非常に好き10%、やや好き15%、ふつう30%
やや嫌い20%、非常に嫌い25%の場合、
95%信頼区間で

非常に好き、非常に好き+やや好きは

どうなるのでしょうか?

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

Re: サンプリング誤差について
  投稿者:青木繁伸 2018/06/10(Sun) 15:45 No. 22566
各カテゴリーは排他なので,
「非常に好き」と「それ以外」で,sqrt(0.1*0.9/600)
「非常に好き+やや好き」と「それ以外」で,sqrt(0.25*0.75/600)
というように考えるということです。

Re: サンプリング誤差について
  投稿者:統計迷い人 2018/06/12(Tue) 10:16 No. 22570 HomePage
青木教授、アドバイス、大変、ありがとうございました。


信頼性係数 KR-20
  投稿者:やまと 2018/06/10(Sun) 18:02 No. 22567
よろしくお願いします。

RでKuder RichardsonのKR20を計算するためのスクリプトは公開されているのでしょうか。教えていただけますでしょうか。

Re: 信頼性係数 KR-20
  投稿者:青木繁伸 2018/06/10(Sun) 21:13 No. 22568
検索しましたか?

https://rdrr.io/github/cddesja/validateR/man/kr20.html
などもあります。

ここに書かれていることが難しければ,ともかく以下のようにすればよいでしょう。

devtools がまだインストールされていないなら,以下の 1 行を一度だけ
install.packages("devtools")

インストルされていた(インストールした)なら,以下の 2 行を一度だけ
library(devtools)
install_github("cddesja/validateR")
そして,kr20 を使うときには R を立ち上げるたびに,以下の 1 行を 一度だけ
library(validateR)
さて,準備が整ったら,データを以下のように準備する
例では 7 人の被検者に 6 項目の変数を 0/1 で準備
x = matrix(c(
1,1,0,1,1,0,
1,0,0,1,1,1,
0,1,1,0,1,1,
0,0,0,1,1,0,
0,1,0,1,0,1,
0,0,0,1,1,0,
0,0,1,0,0,0), byrow=TRUE, ncol=6)
その後,
kr20(x)
とすれば,結果は
>  kr20(x)
[1] 0.1382488
となる。

====

なお,やっていることは簡単なので,自分で以下のような関数を書けば,面倒なことは一切なしで結果が得られる。
計算方法につては Wikipedia にさえ書かれているが,例えば,
https://www.radford.edu/.../KR21%20&%20KR20%20examples.xls
には,Excel での計算例が書かれているので,計算方法を理解すれば以下のような関数を書くのは容易である(R がそこそこできればではあるが)
KR20 = function(x) {
cm = colMeans(x)
k = ncol(x)
k/(k-1)*(1-sum(cm*(1-cm))/var(rowSums(x)))
}

> KR20(x)
[1] 0.1382488

Re: 信頼性係数 KR-20
  投稿者:やまと 2018/06/11(Mon) 07:38 No. 22569
青木先生

大変詳細な説明をいただきありがとうございます。

無事解決できました。


t検定で処理していいでしょうか
  投稿者:コロン 2018/05/27(Sun) 07:54 No. 22556
いつもお世話になっております。

2クラス(実験群と統制群)で、授業を行う前と行った後の平均値の差を見たいのですが、分散分析(被験者間×被験者内)で分析をしたところ、交互作用は非有意、主効果がそれぞれの要因で有意、となりました。

解釈に苦しむところがあり、t検定でpreの平均を調べたところ、差がない、という結果となりました(等分散性の検定もOK)。

この結果を利用して、事前の能力が同じだと仮定して、postのみをt検定で処理してもよろしいでしょうか?実際に行ったところ、有意でした。

アドバイスをお願いいたします。

Re: t検定で処理していいでしょうか
  投稿者:鈴木康弘 2018/05/30(Wed) 11:29 No. 22557
 単に授業を行った後の2クラスの比較をしたいならt検定でいいんでしょうが、
2要因分散分析ではいけないんですか?

Re: t検定で処理していいでしょうか
  投稿者:コロン 2018/05/30(Wed) 19:18 No. 22558
鈴木康弘先生

お返事、ありがとうございます。簡単に説明させていただきますと以下のようになります。

・(要因A)被験者間要因(実験群と統制群)×(要因B)被験者内要因(preとpost)の2要因分散分析(ともに2水準)
・preにはなく、postにおける平均値に差があることが実験仮説
・分散分析の結果、交互作用がなく、要因AおよびBそれぞれに主効果があった(グラフに書きますと、ともに右肩上がりの平行線)

これを、実験仮説に即してうまく解釈できず...つまり、要因Aの主効果に差があると言うことは、「実験群の方が統制群よりも平均が高い」ということですよね?これだと「postで差が出ている」と解釈できないのではないかと考えた次第です。なので、t検定で行ったところ、うまい具合に出てくれたので、上のような質問をさせていただいた次第です。

もしよろしければ再度コメントをいただけますと幸いに存じます。

Re: t検定で処理していいでしょうか
  投稿者:鈴木康弘 2018/06/03(Sun) 07:05 No. 22562
 授業前に有意差はありませんでした、授業後には有意差はありました、と言うことはできるでしょう。

 でも、分散分析での交互作用はどうなのか?と問われたら、ありませんでした、と白状するしかないのでは。


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

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

取り急ぎお礼まで


予測区間の式
  投稿者:hikon 2018/05/16(Wed) 14:29 No. 22552
予測区間の式について教えて頂けましたら幸いです。

検索したところ、下記サイトがあり、式がすべて違い、値も違います。
どれを使えばよいのかわかりません。
コメント頂けましたら幸いです。
https://qiita.com/ksksk/items/75ba95337ccdb32e7cb1
https://ameblo.jp/cabay/entry-12184838325.html
http://heartland.geocities.jp/ecodata222/ed/edj1-1-4.html
https://www.monodukuri.com/gihou/article/946

それはそれとして、今回やりたいことは、
重回帰分析をした後に予測区間をプロットすることです。
お勧めのやり方をご教示頂けましたら助かります。
Rは少し使えます。

お手数おかけいたしますがよろしくお願いいたします。

Re: 予測区間の式
  投稿者:青木繁伸 2018/05/16(Wed) 16:58 No. 22553
2 番目のページと 3 番目のページは,標本平均値の信頼区間と予測区間の話です
1 番目のページと 4 番目のページは,単回帰分析における信頼区間と予測区間の話です

1 番目と 4 番目は,基本的に同じですが,

1 番目のページの式で tn-p-1((1_0.95)/2) の部分は,(1_0.95)/2 は,下に「自由度 n-p-1 の t 分布の (1+0.95)/2(1+0.95)/2 パーセンタイル」と書かれているので誤解を生じるのでしょう。1 番目のページの結果は正しいです。

4 番目のページの式は,t(n-2, α) を「α は有意水準で 5%,(1-α) が信頼水準で,ここへ 95% が入る。」と紛らわしいことを書いているが,信頼水準を 95% にするには t(n-2, α/2) でなければならない。つまり,上側確率が 0.025 である t 値を求める必要がある。αを使ってしまうと間違った答えが出る。

> 重回帰分析をした後に予測区間をプロットする

予測区間の求め方については,
http://sayalabo.com/labo4.html
に説明があります。

ただ,プロットは難しいでしょうね。

Re: 予測区間の式
  投稿者:hikon 2018/05/18(Fri) 08:39 No. 22554
早速のご回答ありがとうございます。
自分の勉強不足を痛感いたしました。
昨日一日勉強いたしました。
まだ理解できていませんが、取り急ぎお礼とご報告をさせて頂きます。
理解が進んだ後、また質問させて頂くこともあると思います。
今後ともよろしくお願いいたします。


ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データであ...
  投稿者:tosh 2018/05/09(Wed) 11:39 No. 22537
ロジスティック回帰分析を行う際、説明変数として用いたい変数が、3種類以上の値を含む名義尺度データであった場合に、それを複数のダミー変数に置き換える方法が、いろいろあることについて悩んでいます。

インターネットで調べられた限りですが、おもに次の(1)〜(3)の方法があるようでした。
これらはどのように使い分ければよいでしょうか。

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

説明変数Aに、"x","y","z" の3種類の名義が含まれている場合、

(1)ダミー変数 { d1, d2, d3 } を作り、変数Aが、

"x" なら { 1, 0, 0 }
"y" なら { 0, 1, 0 }
"z" なら { 0, 0, 1 }

 に置き換える方法。

(2)ダミー変数 { d1, d2 } を作り、変数Aが、

"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { 0, 0 }

 に置き換える方法。

(3)ダミー変数 { d1, d2 } を作り、変数Aが、

"x" なら { 1, 0 }
"y" なら { 0, 1 }
"z" なら { -1, -1 }

に置き換える方法。

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

上に挙げたうち、どれか1つの方法を紹介しているサイトはいくつか見つかるのですが、それぞれの方法の比較や、どう使い分ければよいか、という点についての情報がなかなか見つからず、質問させていただきました。

また、この問題が体系的に解説されている書籍などで、もしおすすめのものがありましたら、ご紹介いただければ幸いです(それぞれの方法の結果の解釈、オッズ比への逆変換の方法、ステップワイズ法で変数を増減する場合にダミー変数はどう扱うのか、などについてもよく知りたいと思っています)。

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:後医は名医 2018/05/09(Wed) 12:52 No. 22538
少なくとも(1)は多重共線性でダメなのでは

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/10(Thu) 12:31 No. 22541
ご返信ありがとうございます。

多重共線性は、互いに相関が強い複数の変数が説明変数の中に含まれていると、起こる問題ですよね。
(1)の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…?


(追記)
すみません、次のように修正しました
修正後「 d1 が 1 なら d2 + d3 は必ず0に決定されてしまうので、」
修正前「 d1 が 1 なら d2 とd3 は必ず0に決定されてしまうので、」

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:青木繁伸 2018/05/10(Thu) 22:50 No. 22542
> 1)の方法の場合、たとえば d1 が1なら d2+d3は必ず0に決定されてしまうので、完全な相関となってしまい、多重共線性につながる、ということなのでしょうか…?

そういうこと

R では,変数Aが名義尺度であるということを
FA = factor(A)
としておけば,
glm(result ~ FA)
でちゃんと処理してくれます。
変数 A は 2 個のダミー変数で表現できます。
3 個以上の変数を使っても,R では 3 個目からは分析に使われません。他の統計ソフトではエラーメッセージを吐いて死んでしまうこともあるでしょう。
独立な2変数変数であれば,どんな値であっても(0 と 1 でなくても -6.23123 と 578.21234 でもなんでも),その変数に対する係数は見かけは違いますが,その係数を使って計算される予測値は全く同じになります(そうでないと困ります。よい結果が出るような変数値のセットがある!なんてことになると,客観性がなくなります)。
ではなぜ,(3)のような変数を作らない!!かというと,出てくる値が3種類になるからです。0と1だと計算も簡単だしわかりやすい。
それにしても 0,1,-1 を与えるというのは,どんなメリットがあるのかなあ?

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/11(Fri) 16:38 No. 22543
ご返信、ありがとうございます。

(3)の意義については、この方法が載っていたpdf資料の内容を後でアップいたします。
(私のパソコンのDLフォルダの整理が悪く、見つかり次第)

通常は (2) の置き換えを行い、(1)はダメ、ということですね。

さっそく、R で、名義尺度の変数を factor型 にし、glm メソッド にかけてみました。
(スクリプトは、別に投稿します)

このスクリプトは次のようなロジスティック回帰分析を行うつもりで書きました。

・データセット: birthwt 低体重出生とそのリスク因子のデータ
・目的変数: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
・説明変数: 母親の年齢、母親の体重、母親の人種、...etc

母親の人種 race {1:白人 2:黒人 3:その他 } が名義尺度なので、スクリプトの中であらかじめ次のように変換しています。

# race を factor型 に変換する
DF1$race <- factor(DF1$race)

glmメソッド の summary をみると、以下に引用したように race2, race3 という2つの変数が生成されているのを確認しましたが、これらは race { 1:白人 2:黒人 3:その他 } をどのように変換した変数なのかを確認する方法がわかりません。

---------------------------------------------------------------
# 結果表示
summary(out1.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF1)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8946 -0.8212 -0.5316 0.9818 2.2125

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.48062 1.19689 0.402 0.68801
age -0.02955 0.03703 -0.798 0.42489
lwt -0.03397 0.01524 -2.229 0.02580 *
race2 1.27226 0.52736 2.413 0.01584 *
race3 0.88050 0.44078 1.998 0.04576 *
smoke 0.93885 0.40215 2.335 0.01957 *
ptl 0.54334 0.34540 1.573 0.11571
ht 1.86330 0.69753 2.671 0.00756 **
ui 0.76765 0.45932 1.671 0.09467 .
ftv 0.06530 0.17239 0.379 0.70484
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.28

Number of Fisher Scoring iterations: 4
---------------------------------------------------------------


ここでの {race2, race3} は、たとえばですが、内部的に

race=1 なら { 1, 0 }
race=2 なら { 0, 1 }
race=3 なら { 0, 0 }

のように変換されているのだと思いますが、それをどのように確認できるのでしょうか? 

str(out1.glm)

としてみたりもしましたが…

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/11(Fri) 16:40 No. 22544
以下、スクリプトです。
----------------------

###################################################
#
# ロジスティック回帰分析例
#
# https://www.gixo.jp/blog/3200/
# より一部改変
#

library(MASS)

str(birthwt)

# *** 目的変数 ***
# low:新生児の体重が2.5kg未満か否か(2.5kg未満が1)
#
# *** 説明変数 ***
# age:母親の年齢
# lwt:母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke:母親の喫煙有無(喫煙有りが1)
# ptl:母親の早産経験の有無(経験有りが1)
# ht:母親の高血圧の有無(有りが1)
# ui:母親の子宮神経過敏の有(有りが1)
# ftv:妊娠第1期における医師の訪問回数
#
# *** 使わない ***
# bwt:新生児の体重(グラム単位)

# ■ 前処理

# bwtは使わないので除く
DF1 <- birthwt[,1:9]

# lwt をポンド単位からグラム単位に変換する
DF1$lwt <- DF1$lwt * 0.454

# race を factor型 に変換する
DF1$race <- factor(DF1$race)

str(DF1)

# ■ ロジスティック回帰

out1.glm <- glm(
low ~ . ,
family = binomial(link = "logit"),
data = DF1)

# 注: モデル式 low ~ . は、目的変数に low、説明変数に low 以外の全ての変数を使うという意味

# 結果表示
summary(out1.glm)

# 各説明変数の係数を指数変換し、オッズ比を表示
exp(out1.glm$coefficients)

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:青木繁伸 2018/05/11(Fri) 21:30 No. 22545
例としてあげられたデータ birthwt$race は整数型で 1, 2, 3 をとる変数でしたが,factor(birthwt$race) で factor型の変数にすると,1, 2, 3 というレベルを持つ factor 変数になります。
> DF1$race
[1] 2 3 1 1 1 3 1 3 1 1 3 3 3 3 1 1 2 1 3 1 3 1 1 3 3 1 1 1 2 2 2 1
途中省略
[161] 1 3 3 3 2 1 3 3 1 1 2 2 2 3 3 1 1 1 1 2 3 3 1 3 1 3 3 2 1
Levels: 1 2 3 ← ここ

このような factor 変数を多変量解析などに使用すると,一番最初の Level が基準値 0 になります。つまり,必要なダミー変数において,取る値の全てが0
一般的には,変数を factor 変数に変換すると,"辞書順の最初の値が基準値 0 になります。
詳しくは factor 関数のオンラインヘルプを参照してもらいたいですが,ちょっと見以下のような結果になりますが,訳ワカメに陥らないようにしっかり理解してください。

> sex = c("male", "female", "U.K.", "male", "male", "female")
> factor(sex, labels=c(1,2,0))
[1] 2 1 0 2 2 1
Levels: 1 2 0
> factor(sex, levels=c("male", "female", "U.K."))
[1] male female U.K. male male female
Levels: male female U.K.
> factor(sex, levels=c("male", "female", "U.K."), labels=c(1,2,0))
[1] 1 2 0 1 1 2
Levels: 1 2 0
さて,このようなことから,R が作っているダミー変数を自分で作るとすると以下のようになります
> DF2 = DF1
> DF2$race = NULL # factor である race を削除
> Race = birthwt$race # 元のデータは整数型 1, 2, 3
> DF2$Race2 = (Race == 2)+0 # Race が 2 のときのみ整数値 1,それ以外なら 0
> DF2$Race3 = (Race == 3)+0 # Race が 3 のときのみ整数値 1,それ以外なら 0
> head(DF2) # Race = 1 のときは DF2$Race2 = DF2$Race3 = 0
low age lwt smoke ptl ht ui ftv Race2 Race3
85 0 19 82.628 0 0 0 1 0 1 0
86 0 33 70.370 0 0 0 0 3 0 1
87 0 20 47.670 1 0 0 0 1 0 0
88 0 21 49.032 1 0 0 1 2 0 0
89 0 18 48.578 1 0 0 1 0 0 0
91 0 21 56.296 0 0 0 0 0 0 1
> out2.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF2)
> summary(out2.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF2)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.8946 -0.8212 -0.5316 0.9818 2.2125

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.48062 1.19689 0.402 0.68801
age -0.02955 0.03703 -0.798 0.42489
lwt -0.03397 0.01524 -2.229 0.02580
smoke 0.93885 0.40215 2.335 0.01957
ptl 0.54334 0.34540 1.573 0.11571
ht 1.86330 0.69753 2.671 0.00756
ui 0.76765 0.45932 1.671 0.09467
ftv 0.06530 0.17239 0.379 0.70484
Race2 1.27226 0.52736 2.413 0.01584
Race3 0.88050 0.44078 1.998 0.04576

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.28

Number of Fisher Scoring iterations: 4
ということで,DF1 を使ったときと同じ結果になりました。

なお,reorder という関数がありますが,この存在を知ると,きっと悪夢を見ます(^_^;)

複雑怪奇な変換規則を理解するよりは,自分で「適切な」ダミー変数を作成する方が幸せかも知れません。

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/12(Sat) 04:06 No. 22547
factor関数をテストしてみると、自動的に定まる Levels の順序が独特ですね…。

> factor(c(100,1,3,1,-5))
[1] 100 1 3 1 -5
Levels: -5 1 3 100 → 数値は昇順なのはまあ分かりますが、

> factor(c("b","a","A","c","A","b","C"))
[1] b a A c A b C
Levels: a A b c C → 文字列は辞書順だが小文字が先…?

最初のLevelがどの値になるか、いまいち確信が持てないので、ご説明にあるように、名義尺度をfactor型に変換する際は、 Levels パラメータ と Labels パラメータをちゃんと与えて、最初のLevelになる値を明示すればいい、と考え、次のように進めました。

> DF3 = DF1
> DF3$race = NULL
> DF3$race = factor(birthwt$race,levels=c(1,2,3),labels=c("白人","黒人","その他"))
> out3.glm <- glm(low ~ ., family = binomial(link = "logit"), data = DF3)

> summary(out3.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF3)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.4806 1.1969 0.40 0.6880
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race黒人 1.2723 0.5274 2.41 0.0158 *
raceその他 0.8805 0.4408 2.00 0.0458 *

この結果が、お示しいただいた out.glm2 の結果と一致するのを見て、やっと気付いたのですが、out1.glm に含まれるダミー変数 race2 , race3 の "2"や"3"は、2:黒人、3:その他 のラベルなんですね。
(No.22543を投稿した時点では、通し番号のようなものかな、と誤解していました)

ダミー変数 { race黒人, raceその他 } は、おそらく

race=1 について { 0, 0 }
race=2 について { 1, 0 }
race=3 について { 0, 1 }

という対応になっていて、

> exp(out3.glm$coefficients)
(Intercept) age lwt smoke ptl
1.62 0.97 0.97 2.56 1.72
ht ui ftv race黒人 raceその他
6.44 2.15 1.07 3.57 2.41

から、

・人種が「黒人」であることは、白人である場合と較べ、低体重児の出生リスクが 3.57倍になる
・人種が「その他」  〃  2.41倍になる

という解釈に行き着くのか、と。

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/14(Mon) 19:48 No. 22550
投稿No. 22537に挙げた(3)の方法
(ダミー変数を -1, 0, 1 などとおく方法)

についてですが、掲載されているのは、ある大学病院で開催された「EBM・統計セミナー」という研修会のPDF資料でした(現在はDLできなくなっていました)。
この資料には、次の2つの変換方法が紹介されており、

ア)基準とするカテゴリ(reference group)を仮定する場合
イ)カテゴリ全体の平均をreferenceとする場合

-------------------
ア)基準とするカテゴリ(reference group)を仮定する場合

例 名義尺度データが{ "薬剤なし", "薬剤P", "薬剤Q" } のような場合

ダミー変数 D1, D2 を

 "薬剤なし" = {0,0}
 "薬剤P"  = {1,0}
 "薬剤Q"  = {0,1}

とおく。ロジスティック回帰の結果、D1, D2 について、

 b1 : D1 に対応する回帰係数
 b2 : D2 に対応する回帰係数

が推定され、

 exp(b1) は"薬剤なし"に対する"薬剤P"のオッズ比、
 exp(b2) は"薬剤なし"に対する"薬剤Q"のオッズ比、

という解釈となる。

--------------------
イ)カテゴリ全体の平均をreferenceとする場合

例 名義尺度データが{ "地域A", "地域B", "地域C" }のような場合

ダミー変数 D1, D2 を

 "地域A" = {1,0}
 "地域B" = {0,1}
 "地域C" = {-1,-1}

とおく。
 
 b1: "地域A"の回帰係数
 b2: "地域B"の回帰係数
 b3: "地域C"の回帰係数

とし、モデル上で、

 b1 + b2 + b3 = 0

という制約を置いたうえで推定する。
計算上は、b1, b2 しか推定されないが、b3 = - b1 - b2 として計算可能である。
exp(b1)をとればオッズ比が得られる。

--------------------
引用ここまでです。

最後の一文にある、「exp(b1)をとれば」は、「exp(-b1-b2)をとれば」ではないのかな…? という疑問がありますが、原文ママです。

なお、この資料の末尾の参考文献には次が挙げられていました。

青木繁伸(2009)「Rによる統計解析」
大橋靖雄(2011)「医師のための臨床統計学基礎編」
丹後俊郎(1996)「路地すてぃく回帰分析SASを利用した統計解析の実際」
豊田秀樹(1992)「原因をさぐる統計学」
服部環・海保博之(1996)「Q&A心理データ解析」
森實敏夫(2004)「入門 医療統計学」

Re: ロジスティック回帰に用いる説明変数が、3種類(以上)の値を含む名義尺度データ...
  投稿者:tosh 2018/05/14(Mon) 20:12 No. 22551
No.22545でお示しいただいた、自分でダミー変数を作る方法を真似して、(3)の方法でダミー変数を作成する処理を追加してみました。

###################################################
#
# ロジスティック回帰分析例 2
# https://www.gixo.jp/blog/3200/
# より一部改変

library(MASS)
str(birthwt)
#
# low: 新生児の体重が2.5kg未満か否か(2.5kg未満が1)
# age: 母親の年齢
# lwt: 母親の体重(ポンド単位)
# race: 母親の人種(1=白人,2=黒人,3=その他)
# smoke: 母親の喫煙有無(喫煙有りが1)
# ptl: 母親の早産経験の有無(経験有りが1)
# ht: 母親の高血圧の有無(有りが1)
# ui: 母親の子宮神経過敏の有(有りが1)
# ftv: 妊娠第1期における医師の訪問回数
# bwt: 新生児の体重(グラム単位)

options(digits = 4)

### DF1 ダミー変数を、Rで自動生成する

# bwtは使わないので除く
DF1 <- birthwt[,1:9]

# lwtをポンド単位からg単位に変換する
DF1$lwt <- DF1$lwt * 0.454

# raceをfactor型に変換する
DF1$race <- factor(DF1$race,levels = c(1,2,3),labels = c("白人","黒人","その他"))
str(DF1)

# ロジスティック回帰
out1.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF1)
summary(out1.glm)
exp(out1.glm$coefficients)

#### DF4 ダミー変数を、白人{1,0} 黒人{0,1} その他{-1,-1} とする方法で自作する

DF4 <- DF1

# ダミー変数に変換する
race = birthwt$race
DF4$race白人 = ifelse(race==1, 1, ifelse(race==2, 0, -1))
DF4$race黒人 = ifelse(race==1, 0, ifelse(race==2, 1, -1))
DF4$race = NULL
str(DF4)

# ロジスティック回帰
out4.glm <- glm(low ~ . , family = binomial(link = "logit"), data = DF4)
summary(out4.glm)
exp(out4.glm$coefficients)

# "raceその他" の偏回帰係数、オッズ比
b3 = - coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"]
names(b3) = c("raceその他_coef")

b3_exp = exp(- coef(out4.glm)["race白人"] - coef(out4.glm)["race黒人"])
names(b3_exp) = c("raceその他_exp(coef)")

-----------------
スクリプトここまでです。
以下、結果を貼ります。
-----------------

> summary(out4.glm)

Call:
glm(formula = low ~ ., family = binomial(link = "logit"), data = DF4)

Deviance Residuals:
Min 1Q Median 3Q Max
-1.895 -0.821 -0.532 0.982 2.212

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1982 1.1509 1.04 0.2978
age -0.0295 0.0370 -0.80 0.4249
lwt -0.0340 0.0152 -2.23 0.0258 *
smoke 0.9388 0.4021 2.33 0.0196 *
ptl 0.5433 0.3454 1.57 0.1157
ht 1.8633 0.6975 2.67 0.0076 **
ui 0.7676 0.4593 1.67 0.0947 .
ftv 0.0653 0.1724 0.38 0.7048
race白人 -0.7176 0.2699 -2.66 0.0079 **
race黒人 0.5547 0.3232 1.72 0.0861 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 201.28 on 179 degrees of freedom
AIC: 221.3

Number of Fisher Scoring iterations: 4

> exp(out4.glm$coefficients)
(Intercept) age lwt smoke ptl ht ui ftv race白人
3.3142 0.9709 0.9666 2.5570 1.7217 6.4450 2.1547 1.0675 0.4879
race黒人
1.7414

> b3
raceその他_coef
0.1629
> b3_exp
raceその他_exp(coef)
1.177

------------
結果は以上です。
資料の説明からすると、

「低体重児の出生確率は、全人種の平均出生確率を基準にすると、
 ・白人の場合、0.4879倍
 ・黒人の場合、1.7414倍
 ・その他の人種の場合、1.177倍」

というような解釈になるのかな、と思いましたが…。
これ以上の詳細は、資料を作った方に直接質問するのが筋かと思いますが、具体的にやってみたことの報告のつもりで投稿しました。


条件付き確率
  投稿者:azyu 2018/05/08(Tue) 00:54 No. 22535
すいません、教えて下さい。確率変数X,Yがあった時に、条件付き確率 P(X|Y)を考えます。P(X=2|Y=1)とすれば、何らかの確率が定数として求まると思うのですが、P(X=2|Y)というものを考えると、Yは確率変数なので、このP(X=2|Y)も確率変数になるのでしょうか?

Re: 条件付き確率
  投稿者:韮澤 2018/05/10(Thu) 08:20 No. 22539
ベイズの定理に関わる式 P(X|Y)は、Yがある値に確定した後での確率P(X)ですから、先にXが確定するなら、YはPに関与しないと考えるべきでしょう。従って、P(X=2)と等価の確定値でしょう。

Re: 条件付き確率
  投稿者:azyu 2018/05/12(Sat) 12:32 No. 22548
どうもありがとうございました。
P(X=2|Y=1)もP(X=2|Y=2)も値としては同じになるのですかね。というか、そもそも、この「Y=」の部分は、関与しないわけですから、記法としての意味がないということなのでしょうか。

Re: 条件付き確率
  投稿者:韮澤 2018/05/14(Mon) 08:25 No. 22549
あくまで、Yが事前に確定してこその、P(X|Y)です。おっしゃる通り、Xが先に確定するなら、記法として成り立ちません。


繰り返しのある共分散分析
  投稿者:ならしか 2018/05/01(Tue) 09:51 No. 22521
よろしくお願いします。

指導法の効果を見るために,実験群と統制群のクラスで,それぞれテストを3回行いました。

初回のテストは実験前に,2回目は授業終了直後,3回目は遅延テストとして行いました。

データ分析の段階で,初回のテストの等分散性が確保できませんでした。この場合,2要因分散分析ではなく,共分散分析をしなければならないところまではなんとなく理解しています。

【質問1】分散分析の場合は,すべてのデータを一括で処理できるのですが,今回のデータは,

(1)初回と2回目の共分散分析
(2)初回と3回目の共分散分析

と分析しなければいけないのでしょうか?

【質問2】2回目と3回目の差についても興味があるのですが,この場合はどういう方法で見ることが可能なのでしょうか?

よろしくお願いします。

Re: 繰り返しのある共分散分析
  投稿者:鈴木康弘 2018/05/08(Tue) 11:11 No. 22536
 あまり自信はないのですが..

 等分散性が確保出来ない場合に共分散分析をしなければならないって、本当でしょうか?
 そういう場合にはいい方法はないので、2要因分散分析でごまかすしかないような...


ログランク検定におけるサンプルサイズ算出時の打ち切り分布
  投稿者:kamo 2018/05/01(Tue) 15:28 No. 22522
お世話になります。

ログランク検定のサンプルサイズを計算する際、以下の理論式(プログラム)が広く使われていると学習いたしました。理解を深めるために、この計算結果で得られる理論値(サンプルサイズ)、と検出力の関係をシミュレーションで再現したいと考えております。
例えば、以下の例ですと、受け入れR年、追跡期間T年とした場合、各群57例という結果が得られます。ここで、群ごとに57例分の生存時間を乱数で発生させて、1000回か10000回程度、検定を回したときに、検出力(有意になった回数)が実際に8割を超えることを確認したいと考えております。

生存時間は、ハザードをパラメータとした指数分布に従うことを仮定していることは理解しております。しかし、打ち切り症例が発生する確率と、打ち切り例の場合の打ち切りまでの時間に、どのような分布とパラメータを仮定すべきなのか分からず、煮詰まっております。この点について、ご教授頂けないでしょうか?

お手数をおかけして恐縮でございますが、ご返信頂けますと幸いです。また、質問に不備がございましたら、ご指摘頂けますと幸いです。
#  R      : 受け入れ期間
# Tr : 追跡期間
# TIME : TIME年生存率
# A : 100A(%)生存率(<1)
# B : 100B(%)生存率(<1)
# ALPHA : 有意水準(<1)
# BETA : 検出力(1-BETA)
# S : 例数比(Na/Nb)

R <- 3
Tr <- 5
TIME <- 5
A <- 0.2
B <- 0.05
ALPHA <- 0.05
BETA <- 0.2
S <- 1

LA <- -log(A)/TIME
LB <- -log(B)/TIME
LBER <- (LA + LB)/2
QA <- 1/(1 + S)
QB <- S/(1 + S)
ZA <- qnorm(1 - ALPHA/2)
ZB <- qnorm(1 - BETA)
P00 <- LBER * LBER
P01 <- LA * LA
P02 <- LB * LB
P10 <- LBER^2/(1 -exp(-LBER*Tr))
P11 <- LA^2/(1-exp(-LA*Tr))
P12 <- LB^2/(1-exp(-LB*Tr))
P20 <- LBER^3*Tr/(LBER*Tr-1+exp(-LBER*Tr))
P21 <- LA^3*Tr/(LA*Tr-1+exp(-LA*Tr))
P22 <- LB^3*Tr/(LB*Tr-1+exp(-LB*Tr))
P30 <- LBER^2*(1-(exp(-LBER*(Tr-R))-exp(-LBER*Tr))/(LBER*R))^(-1)
P31 <- LA^2*(1-(exp(-LA*(Tr-R))-exp(-LA*Tr))/(LA*R))^-1
P32 <- LB^2*(1-(exp(-LB*(Tr-R))-exp(-LB*Tr))/(LB*R))^-1
N <- ((ZA*sqrt(P00*(1/QB+1/QA)) + ZB*sqrt(P01/QA+P02/QB))/abs(LA-LB))^2/(1+S)
M <- ((ZA*sqrt(P10*(1/QB+1/QA)) + ZB*sqrt(P11/QA+P12/QB))/abs(LA-LB))^2/(1+S)
O <- ((ZA*sqrt(P30*(1/QB+1/QA)) + ZB*sqrt(P31/QA+P32/QB))/abs(LA-LB))^2/(1+S)
L <- ((ZA*sqrt(P20*(1/QB+1/QA)) + ZB*sqrt(P21/QA+P22/QB))/abs(LA-LB))^2/(1+S)

c(S=S, 受け入れ期間=R, 生存率=TIME, ハザードA=LA, ハザードB=LB,
ALPHA=ALPHA, ZA=ZA, BETA=BETA, ZB=ZB,
"R=0,T=Inf"=N, "R=0, T>0"=M, "0<R<T"=O, "O<R=T"=L)
参考:[新版]実用SAS生物統計ハンドブック[SAS®9.4/R3.2.0対応](臨床評価研究会)

Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布
  投稿者:青木繁伸 2018/05/02(Wed) 23:23 No. 22524
> 打ち切り症例が発生する確率と、打ち切り例の場合の打ち切りまでの時間に、どのような分布とパラメータを仮定すべきなのか

シミュレーションデータの各例は,以下の順で生成されるのでは?
(1) 受入期間内に一様分布で研究に参入する
(2) エンドポイントまでの期間は指数分布で決まる
(3) エンドポイント発生が追跡期間後になるなら「打切り例」,前なら「故障例」
(4) 「打切り例」なら,研究対象であったのは研究参入から追跡期間終了までの期間
(5) 「故障例」なら,研究対象であったのは研究参入からエンドポイントまでの期間
これで,「打ち切り症例が発生する確率と、打ち切り例の場合の打ち切りまでの時間に、どのような分布とパラメータを仮定すべきなのか」というのは無関係,意味がないことになる。

なお,あなたがいっている「打ち切り症例」が「脱落例(追跡不能例)」を含んでいるなら,そもそもそういうものは考える必要がない。例えば,2群の平均値の差のサンプルサイズの決定などでも,測定不能,測定不全,拒否などによる例数減は考えないでしょう?
受入開始         受入終了        追跡終了
+..........................+.......................+
o--------------x 受入終了前に死亡
o-------------------------------x 追跡終了前に死亡
o-------------------------------------@=======x 打切り例: 追跡終了後に死亡
- と = の和が指数分布で決定されるエンドポイントまでの期間
- が研究対象であった期間

Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布
  投稿者:kamo 2018/05/03(Thu) 21:10 No. 22528
青木 繁伸 先生

ご返信頂きましてありがとうございます。

打ち切り例は、脱落例を含むという前提で質問申し上げておりました。
二群の平均値の検定の場合は、(補完などの特別な処理をしなければ)脱落例は解析対象に含められませんが、生存時間解析においては、そういった対象も情報として利用するようですので、
サンプルサイズの計算においても実験途中の脱落を考慮していると理解しておりました。

ご教授頂いた内容を踏まえ、Rコードを書いてみました。
しかし、各群57例の場合の検出力は70~75%程度にとどまり、上の計算結果と整合しませんでした。
どこか誤りがありますでしょうか。
恐れ入りますが、改めてご指導賜りたく、何卒よろしくお願い申し上げます。
#初期化
sign <- NULL
pval <- NULL

#ハザードを計算
B <- 0.05
A <- 0.2
LB <- -log(B)/(5*365.25)
LA <- -log(A)/(5*365.25)

n <- 57*2

accrualtime <- 3*365.25#登録期間
followuptime <- 5*365.25#観察期間※登録期間を含む

for(N in 1:1000){
#初期化
x1 <- NULL#群
time <- NULL#イベントまでの時間
timec <- NULL#打ち切りまでの時間
randomized.at <- NULL
dt <- NULL
status <- NULL
i <- 0
for(cal in 1:accrualtime){#登録期間の初日を1として、
entry <- rbinom(size=1,n=1,prob=n/accrualtime)#当該時点に登録があるか(登録期間を通じて一様に集積されることを仮定する
if(entry==1){ #当該時点に登録ありの場合
randomized.at <- c(randomized.at,cal) #登録日
i <- i+1 #例数カウンター
x1[i] <- rbinom(size=1,n=1,prob=0.5) #ランダム化

if(x1[i]==0){time[i] <- rexp(n=1,rate=LB)}#対照群の場合の生存時間
else {time[i] <- rexp(n=1,rate=LA)}#治療群の場合の生存時間

timec <- followuptime-randomized.at[i]#試験終了による打ち切り

if(time[i] < timec){ #打ち切り前にイベント発生
status[i] <- 1 #イベント
}

if(timec < time[i]){ #イベント発生前に打ち切り
status[i] <- 0 #打ち切り
time[i] <- timec #打ち切り症例には打ち切りまでの時間を代入置換
}

}

else {}#登録なしの場合はレコード発生させない

if (i==n) break#予定症例数に達したら終了

}

#ログランク検定
result <- survdiff(Surv(time)~x1)
#p値を直接取り出せないので、検定統計量を格納し、そこからp値を導出
pval[N] <- pchisq(result$chisq, length(result$n)-1, lower.tail = FALSE)
if (pval[N] < 0.05){sign[N] <- 1} else{sign[N] <- 0}
}

sum(sign)/N
hist(pval)
> sum(sign)/N
[1] 0.724

Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布
  投稿者:青木繁伸 2018/05/04(Fri) 14:35 No. 22529
一番の間違いは,
survdiff(Surv(time) ~ x1)
としたところ。これは,
survdiff(Surv(time, status) ~ x1)
でしょう。
もう一つは,
x1[i] <- rbinom(size=1,n=1,prob=0.5) #ランダム化
ですが,これでは群ごとのサンプルサイズが同じになりません。
その他も,無駄を省き,分かりやすく書いて,ベクトル演算による高速化を図ったので,参考にしてみてください。
f = function(n.each, B, A, accrualtime, followuptime, REPETITION = 1000) {
library(survival)
LB = -log(B) / followuptime
LA = -log(A) / followuptime
n = n.each * 2
x1 = rep(0:1, each = n.each)
pval = numeric(REPETITION)
for (N in 1:REPETITION) {
# 研究に登録される時刻(研究開始時刻を 0 とする)
entry = runif(n, min = 0, max = accrualtime) # 一様分布に従う
# 追跡可能如何に関わらず,登録時刻からエンドポイントまでの時間
# 前半に B 群,後半に A 群をまとめて生成しても一般性を損なわない
surv.time = c(rexp(n.each, rate = LB), rexp(n.each, rate = LA))
# 研究開始時刻(0)からエンドポイントまでの時間は entry + surv.time
# その時間が観察期間より短い場合は故障例(status = 1)
status = (entry + surv.time) < followuptime
# 故障例ならエンドポイントまでの時間は surv.time そのもの
# 打切り例 status = 0 なら観察期間終了までの時間
time = ifelse(status, surv.time, followuptime - entry)
# ログランク検定
# result = survdiff(Surv(time) ~ x1) # ここが間違い
result = survdiff(Surv(time, status) ~ x1) # こちらが正しい
# p値を直接取り出せないので、検定統計量を格納し、そこからp値を導出
pval[N] = pchisq(result$chisq, 1, lower.tail = FALSE)
}
mean(pval < 0.05)
}

n.each = 57
B = 0.05
A = 0.2
# 以下の2項目の単位は何でもよい
accrualtime = 3 * 365.25 # 登録期間
followuptime = 5 * 365.25 # 観察期間 ※登録期間を含む
REPETITION = 1000 # シミュレーション試行回数

> f(n.each, B, A, accrualtime, followuptime, REPETITION) # 実行結果は毎回異なる
[1] 0.823
> f(n.each, B, A, 3, 5, 10000)
[1] 0.8253

Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布
  投稿者:kamo 2018/05/07(Mon) 18:36 No. 22533
青木 繁伸 先生

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

>一番の間違いは,
>survdiff(Surv(time) ~ x1)
>としたところ。これは,
>survdiff(Surv(time, status) ~ x1)
>でしょう。
基本的なところを間違っておりました。大変失礼いたしました。
この点を修正したところ、理論式と概ね整合のとれる結果が得られました。

>もう一つは,
>x1[i] <- rbinom(size=1,n=1,prob=0.5) #ランダム化
>ですが,これでは群ごとのサンプルサイズが同じになりません。
乱数なので同じサンプルサイズになるとは限らない、といういことですね。ご指摘ありがとうございます。

シミュレーション回数10000回で、頂いたプログラムを何度か回してみましたが、検出力は全て82%と少し、という感じでございました。
80~80.5%程度の範囲に収まるのかなと勝手に想像しておりましたが、これはやはり、シミュレーションと理論式の違い(誤差?)と考えてよいのでしょうか?

Re: ログランク検定におけるサンプルサイズ算出時の打ち切り分布
  投稿者:青木繁伸 2018/05/07(Mon) 21:06 No. 22534
他のサンプルサイズ決定については,
http://www.st.nanzan-u.ac.jp/info/ma-thesis/2008/TANAKA/m07mm028.pdf
が,参考になるのではないかと思います。
私はとても追い切れませんが,必要とするあなたらならそれぞれを試してみるとよいのではないでしょうか?
その他にも,いくつかあるかも知れませんね(追求する情熱が私にはない)
何かわかったら,ここで教えてください。


R 記事No. 22299 について、csv出力の追加
  投稿者:明石 2018/05/04(Fri) 18:49 No. 22530
青木先生 様;

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

青木先生からご教示いただいたノウハウを、読み返して整理しています。

[22299] R ftable()関数      投稿者:明石 投稿日:2017/03/11(Sat) 19:03

タイタニックのデータを題材にして、多重クロス表をftable()関数で整形出力する方法について、ご教示を頂戴しました。
大変に重宝しています。
ありがとうございました。
ftable(dat, row.vars=c("Sex", "Age", "Class"), col.vars="Alive")

Alive 死亡 生還
Sex Age Class
女性 子ども 1等 0 1
2等 0 13
3等 17 14
乗務員 0 0
大人 1等 4 140
2等 13 80
3等 89 76
乗務員 3 20
男性 子ども 1等 0 5
2等 0 11
3等 35 13
乗務員 0 0
大人 1等 118 57
2等 154 14
3等 387 75
乗務員 670 192
上記コンソール出力と同じようなレイアウトで、csv出力したいと思って調べましたが、分かりませんでした。

出力したcsvファイルをEXCELで開くとコンソールへ出力と同じようなレイアウトになることを期待してます。

ご教示を頂戴できれば、大変に助かります。
毎々、お手数をおかけいたします。
何卒どうぞよろしくお願いいたします。

Re: R 記事No. 22299 について、csv出力の追加
  投稿者:青木繁伸 2018/05/04(Fri) 21:43 No. 22531
私は,Excel がないので,LibreOffice で確認したのですが,
a = ftable(dat, row.vars=c("Sex", "Age", "Class"), col.vars="Alive")
write.table(format(a, quote=FALSE), "temp.csv", sep=",", row.names=FALSE, col.names=FALSE)
で良いと思います。

Re: R 記事No. 22299 について、csv出力の追加
  投稿者:明石 2018/05/05(Sat) 08:51 No. 22532
青木先生 様;

お忙しいところを失礼いたします、明石と申します。
今回も助けていただきました。
ありがとうございました。

出力したcsvファイルをEXCELで開いて、所望したフォーマットであることを確認しました。

本体の機能にちゃんと備わっているのですね。
これを機に、関数の引数(ふだんはデフォルト使用で意識しません)をヘルプで調べるようにします。

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


R 記事No. 22000 について、出力項目の追加
  投稿者:明石 2018/05/03(Thu) 11:22 No. 22525
青木先生 様;

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

約2年前、青木先生から以下の有難いご教示を頂戴し、大変に重宝しています。

Re: R str()関数の結果の取り出し、整形出力
  投稿者:青木繁伸 2016/05/14(Sat) 05:12 No. 22001

a = data.frame(names(dat), sapply(dat, class))
colnames(a) = c("名前", "型")
print(a, row.names=FALSE)

などはどうでしょうか

ーーーーーーーーーーーーーーーーーー
名前、型の2列に、欠損値NAの個数 も追加した3列の出力 を検討しています。

まず、欠損値NAの個数 KAZU を計算します。
kazu <- numeric(ncol(dat))
for (i in 1:ncol(dat)) {
kazu[i] <-sum(is.na(dat[,i]))
}

次に、3列からなるデータフレームを作成します。
a = data.frame(names(dat), sapply(dat, class), kazu)
colnames(a) = c("名前", "型", "NA数")
print(a, row.names=FALSE)

data.frame()の中に、
欠損値NAの個数を算出する式を直接に書きたいと思っていますが、力不足で書くことができません。

ご教示いただけましたら助かります。
どうぞ、よろしくお願いいたします。

Re: R 記事No. 22000 について、出力項目の追加
  投稿者:青木繁伸 2018/05/03(Thu) 11:33 No. 22526
どうしても data.frame 関数中に書きたいならば,
a = data.frame(names(dat), sapply(dat, class), colSums(is.na(dat))
colnames(a) = c("名前", "型", "NA数")
print(a, row.names=FALSE)


Re: R 記事No. 22000 について、出力項目の追加
  投稿者:明石 2018/05/03(Thu) 14:26 No. 22527
青木先生 様;

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

お示しをしてくださいましたRコードは、よく理解できました。

どうしても data.frame 関数中に書きたいと思いますので、
勉強になりました。
ありがとうございました。


再現性や正確性の評価
  投稿者:Gummy 2018/05/01(Tue) 19:43 No. 22523
御世話になります。

3人が行った視覚的判定のバラツキについて、統計学的な検定方法を御教示いただけましたらと思い、投稿いたしました。

サンプルを30回自由に粉砕した後、その粉砕片を肉眼的に確認して評価しました(スコア:0-9)。3人がスコアを付けし、1000パターンのサンプルに対して行いました。なお粉砕片はその後溶媒に溶かし、物質の溶出濃度を測るため、実際の粉砕程度を後で確認しています。
スコアを付ける際は、スコア0-9の粉砕片の見本写真があり、その見本写真と実際の粉砕片を見比べてスコアを付けるため、バラツキは小さいということを証明したいのですが、それを証明する方法が分かりません。

級内相関とCV値を使用して、バラツキは小さかったと結論づけしましたが、スコアの正確性と再現性を証明するために十分なのか自信がありません。

恐れ入りますが御教示いただけましたら幸いです。よろしくお願い申し上げます.


サンプル数の異なる3群の検定
  投稿者:三好 2018/04/28(Sat) 13:24 No. 22518
クラスカル-ウォリス検定を考えております。しかし、1群はサンプル数が非常に小さい(20くらい)ことが予想され、もう1群も60前後と予想しています。このような場合、3群目は、200人など多いほうが良いのでしょうか。それとも、60前後など近い人数が良いのでしょうか。お教えいただければありがたいです。よろしくお願いいたします。

Re: サンプル数の異なる3群の検定
  投稿者:青木繁伸 2018/04/29(Sun) 07:36 No. 22519
実験デザインの段階では,各群のサンプルサイズは成るべく同じにするようにする
一般の統計調査の場合には,層別の母集団サイズが異なる場合などは特に,サンプルサイズを揃えることが難しい或いは意味がないこともある
標本が得られたあとで群を分けて違いを検討する場合には,サンプルサイズは揃えようがない

結論として,群ごとのサンプルサイズが違っても,そのまま検定する
間違っても取ったデータを捨てて,サンプルサイズを揃えたりしないように

Re: サンプル数の異なる3群の検定
  投稿者:三好 2018/04/29(Sun) 11:13 No. 22520
青木先生、
ご回答をありがとうございました。
確かに母集団のサイズが異なることを考慮しなければなりませんね。その点は、失念いたしておりました。
今回は第1群の母集団が小さく、サンプル数を増やすのが難しいと思われますので、第2群と第3群は揃えるように設計して見ます。
ご教示どうもありがとうございました。

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

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