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

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


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

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


ケース・コントロール群の設定について
  投稿者:くじら 2017/11/17(Fri) 16:09 No. 22476
お世話になります。

現在、既に調査が終わったデータがあり、論文を書くことになっております。

例えば、約5000人の回答があり、その中で、500人が介護をしていると答えました。
その人達と、介護をしていない人たちの健康状態の違いを解析したいのですが、コントロールとしては、5000人の中から、500人の介護をしていないと答えた人々をコントロール群としたいと考えております。

特に、性別、年齢をマッチさせたいと考えておりますが、クロス集計などで、このような方法で対照群を設定してよろしいでしょうか?

また、対照群は、介護していないと答えた人の中から、それぞれ性・年齢のマッチする人を、SPSSの無作為抽出の機能で選択しようと考えておりますが、何か他に方法がございますでしょうか?

お忙しい中大変申し訳ないのですが、教えていただければ幸いです。お手数をお掛け致しますが、どうぞよろしくお願い致します。

Re: ケース・コントロール群の設定について
  投稿者:青木繁伸 2017/11/17(Fri) 21:10 No. 22477
大変いいにくいことですが,あなたの調査方針・分析方針は不適切です。はっきりいって,涙が出ます。

仮説を立てて,その可否を確かめるために行うのが調査(データ収集)です。
調査をしてしまった後で,解析デザインを作り上げるのは本末転倒です。
対象が500人なら,対照も500人にすべきで,何のために無駄に4000人調査したのでしょうか。あなたは,500+500の調査データだけを利用しようとしているのですよ。
5000人の回答が得られる調査をするマンパワーがあるなら,介護をしている人,していない人それぞれを選んで(特に,なんとしても最初から性・年齢をマッチさせて)2500人ずつの調査をすればどれだけ素晴らしかったか

> 性別、年齢をマッチさせたいと考えておりますが

調査開始時に,調査対象を選ぶときにその条件を発動すべきでしたね。

そのような調査デザインならば,ケース・コントロール研究として分析できたのに(4000人の無駄もなかったのに)

まあ,ああだこうだ,死んだ子の年を数えても仕方ないので,横断研究ということで割り切って,全データを使って分析しましょ(別に,それが悪いことではなく,データ収集そのものが横断研究だったのだから仕方ないでしょう)
それしかいえません。
あえて,マッチングするためにデータを捨てるということだけはおすすめできません
だって,どのデータを捨てるか決められますか?
いい結果が出なかったら,データを選び直す?そんなこと,してはだめですからね。

Re: ケース・コントロール群の設定について
  投稿者:くじら 2017/11/18(Sat) 21:08 No. 22478
お返事、どうもありがとうございました。

たら、れば、で、2500名で性・年齢をマッチしたケース・コントロールができるとさえ考えられるお立場にいらっしゃる方には想像できない程のことが、現実には生じることがあり、その対応に苦慮する者の気持ちなどはご理解いただけないのでしょう。教科書どおりの理想的な調査ができる幸せな環境で研究をされているようで、とてもうらやましく思います。

ご指摘頂きましたように泣いてもどうしようもなく、また、同僚のやったこととはいえ5000名の方々に大変申し訳なく思い、こちらでもアドバイス頂ければと思い投稿させていただいた次第です。

Re: ケース・コントロール群の設定について
  投稿者:青木繁伸 2017/11/18(Sat) 22:17 No. 22479
運営上の問題で,No. 22477 の発言者名が別人になっておりましたことを,お詫びいたします。発言者は,青木繁伸 でした。申し訳ありません。
誤った発言者名が表示されてしまった方には,深くお詫び申し上げます。。
また,その記事に対して,発言者が違ったことにより本来の意図は違った発言を誘発してしまったかもしれないことについては,まことに申し訳ありません。
発言を消去いただいても差し支えありません。申し訳ありませんが,対処方よろしくお願いします。
不快な思いをさせてしまったたことについてはお詫び申し上げます。

Re: ケース・コントロール群の設定について
  投稿者:鈴木 康弘 2017/11/19(Sun) 09:02 No. 22480
 マッチングというのは原則としてデータ収集の段階でやるものですね。

 データを集めてしまってから性・年齢の交絡因子を調整したい、という話なら、
多変量解析じゃないでしょうか。

Re: ケース・コントロール群の設定について
  投稿者:くじら 2017/11/19(Sun) 22:52 No. 22482
青木先生

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

「いい結果が出なかったら,データを選び直す?そんなこと,してはだめですからね。」という助言であっても、5000名の人々に大変申し訳ないという思い、ネットですから詳細は書けないなどの制約の中で、それでも、何か、アドバイスが頂きたかったのです。

鈴木康弘様

コメント、どうもありがとうございました。
多変量解析は検討しております。詳細は説明できないのですが、年齢の幅が15歳から80歳とひろく、一般的な調査データではない為、また、介護を例として考えますと、介護している 500 (女性・年齢が50〜70歳台に偏る)と していない 4500(性別問わず15歳から80歳まで)という間抜けなデータになっており、横断研究として、ロジスティックなどで性・年齢をとにかく調整してもいいものやら、(実際には、その通りの横断研究ではありますが、)そうするしかないのかと思い始めておりますが、扱いに苦慮しております。


Rの'irr'パッケージ(新たなスレッドで)
  投稿者:0669SABON 2017/11/14(Tue) 13:10 No. 22474
再度、検定力分析について調べ直しましたがまだ十分に理解できておりません。
またirr パッケージのオンラインヘルプに引用されている文献ですが取り寄せることができておらず、まだ読んでおりません。
このような状況で質問するのは無礼だと承知していますが、ご教授いただけると幸甚です。

irrパッケージ のN2.cohen.kappaを使用したいのですが

オンラインヘルプの
mrg a vector of marginal probabilities given by raters
ということがわかりません。
そのため
Examples
require(lpSolve)
# Testing H0: kappa = 0.4 vs. HA: kappa > 0.4 (=0.6) given that
# Marginal Probabilities by two raters are (0.2, 0.25, 0.55).

0.2と0.25と0.55の数字が何を現しているかわかりません。
勉強不足で理解が及ばないためご教授頂けると大変助かります。
何卒よろしくお願いいたします。

Re: Rの'irr'パッケージ(新たなスレッドで)
  投稿者:青木繁伸 2017/11/14(Tue) 13:36 No. 22475
marginal probabilities というのは,「周辺確率」ということですが,統計学を知らないと難しいかも知れませんね。
集計表で marginal sum というのは「周辺和」といいますが,一般用語では「合計欄」とか「行和」,「列和」といったほうが分かりやすいのでしょう。
irr の kappa2 の example で扱われている diagnose データセットで例を挙げると
> data(diagnoses)
> # Unweighted Kappa for categorical data without a logical order
> kappa2(diagnoses[,2:3])
Cohen's Kappa for 2 Raters (Weights: unweighted)

Subjects = 30
Raters = 2
Kappa = 0.631

z = 7.56
p-value = 4.04e-14
この分析例では,rater4 と rater5 を対象としています
> head(diagnoses[,c(4, 5)])
rater4 rater5
1 4. Neurosis 4. Neurosis
2 5. Other 5. Other
3 3. Schizophrenia 3. Schizophrenia
4 5. Other 5. Other
5 4. Neurosis 4. Neurosis
6 3. Schizophrenia 3. Schizophrenia
集計表は table(diagnoses[,c(4, 5)]) で求められますが,簡単に表示すると
> t=table(diagnoses[,c(4, 5)])
> colnames(t) = rownames(t) = 1:5
> t
rater5
rater4 1 2 3 4 5
1 1 0 0 0 1
2 0 1 0 0 0
3 0 0 6 0 1
4 0 0 0 12 1
5 0 0 0 0 7
これに「周辺和」を付加すると,
> addmargins(t)
rater5
rater4 1 2 3 4 5 Sum
1 1 0 0 0 1 2
2 0 1 0 0 0 1
3 0 0 6 0 1 7
4 0 0 0 12 1 13
5 0 0 0 0 7 7
Sum 1 1 6 12 10 30
となります。Sum と書かれた行と列にあるのが「周辺和」です。
「周辺確率」はそれをサンプルサイズ(今の例では 30)で割ったものです
> rowSums(t)/30 # rater4 に対する周辺確率
1 2 3 4 5
0.06666667 0.03333333 0.23333333 0.43333333 0.23333333
> colSums(t)/30 # rater5 に対する周辺確率
1 2 3 4 5
0.03333333 0.03333333 0.20000000 0.40000000 0.33333333
それぞれの周辺確率は本来は同じ理論的な確率にしたがうものでしょう。N2.cohen.kappa で指定するのは,この「理論的な(周辺)確率」です。本来は得られたデータから決めるものではありませんが,取りあえずは 2 人の評価者の平均を参考にするということもあるでしょうし,このデータセットでは6人の評価者の平均を使うということもあるでしょう。先行研究の結果を使うこともあるでしょう。

Re: Rの'irr'パッケージ(新たなスレッドで)
  投稿者:0669SABON 2017/11/19(Sun) 10:37 No. 22481
青木先生

 ありがとうございます。
 返信が遅くなり申し訳ございません。
 また周辺確率という新しい言葉を聞いて再度調べ直しています。
 
 先のスレッドを見て、私自身もサンプルサイズなどを検討せず研究を進めていたため反省しています。

 研究方法、研究倫理を見直します。
 統計方法以前の問題ですね。

 
 ありがとうございました。
 またご相談させて頂きます。


Rの'irr'パッケージ
  投稿者:0669SABON 2017/11/06(Mon) 05:40 No. 22468
Cohen's Kappa のsample sizeを調べています。
Rの'irr'パッケージを取り込むことができたのですが、

> # Testing H0: kappa = 0.7 vs. HA: kappa > 0.7 given that
> # kappa = 0.85 and both raters classify 50% of subjects as positive.
> N.cohen.kappa(0.5, 0.5, 0.7, 0.85)

この0.5と0.7と0.85の数字の意味が分からなく、自験例ではどの数字を入れて良いか判断できません。

初歩的な質問で申し訳ないですがご教授頂けると幸甚です。

Re: Rの'irr'パッケージ
  投稿者:青木繁伸 2017/11/06(Mon) 08:08 No. 22469
オンラインヘルプは見ましたか?
Usage

N.cohen.kappa(rate1, rate2, k1, k0, alpha=0.05,
power=0.8, twosided=FALSE)
Arguments

rate1
the probability that the first rater will record a positive diagnosis

rate2
the probability that the second rater will record a positive diagnosis

k1
the true Cohen's Kappa statistic

k0
the value of kappa under the null hypothesis

Re: Rの'irr'パッケージ
  投稿者:0669SABON 2017/11/07(Tue) 05:29 No. 22470
青木先生
ありがとうございます。
オンラインヘルプを確認しました。

再度、質問をさせて頂きたいのですが
rate1
the probability that the first rater will record a positive diagnosis
というのは2択であれば50%→0.5という意味でしょうか。
7択であれば14%→0.14となるのでしょうか。

また
# Testing H0: kappa = 0.7 vs. HA: kappa > 0.7 given that
# kappa = 0.85 and both raters classify 50% of subjects as positive

2行目のkappa=0.85は080以上となるように求めたいので0.80としますが、
1行目のTesting H0:kappa =0.7とありますが、この0.7は何を指しているのでしょうか。

内容を理解しておらず、質問も分かり難いかと思いますがご教授頂けると幸いです。

Re: Rの'irr'パッケージ
  投稿者:青木繁伸 2017/11/07(Tue) 11:17 No. 22471
> rate1
> the probability that the first rater will record a positive diagnosis
> というのは2択であれば50%→0.5という意味でしょうか。
> 7択であれば14%→0.14となるのでしょうか。

N.cohen.kappa のオンラインヘルプの1行目に,

sample size estimator for the Cohen's Kappa statistic for a binary outcome
~~~~~~~~~~~~~~~~~~~~
とありますから前者ですね。(N2.cohen.kappa なら more than one category)

「ネガティブな反応」が多い状態とはどういうことか,つまりそのような状態は,評価者がやる気がない(反抗的)ということで,信頼できないデータである可能性があるということ(逆にポジティブな反応が多すぎても,いい加減な肯定と言うことで信頼できないであろう)。50-50 というのは,ポジティブ・ネガティブの基準が平均的(中央値!)ということで客観的で妥当な評価が行われていることを意味する。

> # Testing H0: kappa = 0.7 vs. HA: kappa > 0.7 given that
> # kappa = 0.85 and both raters classify 50% of subjects as positive
> 2行目のkappa=0.85は0.80以上となるように求めたいので0.80としますが、

「0.80以上となるように求めたい」のではなく,個々で指定する値は,実際にあなたのデータで得られたκの値

> 1行目のTesting H0:kappa =0.7とありますが、この0.7は何を指しているのでしょうか。
文字通り帰無仮説のκ。これもあなたが決めるもの。
つまり,あなたは「母数がこれ以上(0.7 以上)かどうかを検定したい」というわけでしょ。「いやいや,0.6 以上あれば十分」ということなら 0.6 を指定するわけです。

「κ統計量の検定・検定力分析」ということではなく,検定,検定力分析の「一般理論」をおさらいすべきですね。基礎がわかれば,応用はおのずとわかります。基礎がわかっていないと,応用はできません。
irr パッケージのオンラインヘルプに引用されている
Cantor, A. B. (1996) Sample-size calculation for Cohen's kappa. Psychological Methods, 1, 150-153.
Flack, V.F., Afifi, A.A., Lachenbruch, P.A., & Schouten, H.J.A. (1988). Sample size determinations for the two rater kappa statistic. Psychometrika, 53, 321-325.
を読むことも必要でしょう。


> N.cohen.kappa(0.5, 0.5, 0.7, 0.85)
[1] 96
これが例題
> N.cohen.kappa(0.5, 0.5, 0.6, 0.85)
[1] 38
帰無仮説の kappa を下げる(0.6)といことは,サンプルサイズは少なくてよいということになる
> N.cohen.kappa(0.5, 0.5, 0.8, 0.85)
[1] 753
kappa を上げる(0.8)ということは,サンプルサイズはより多く必要ということになる
> N.cohen.kappa(0.5, 0.5, 0.7, 0.75)
[1] 1142
実際のデータから得られたκが小さいと(帰無仮説のκとの差が小さくなる),サンプルサイズはより多く必要となる
> N.cohen.kappa(0.2, 0.2, 0.7, 0.85)
[1] 153
> N.cohen.kappa(0.8, 0.8, 0.7, 0.85)
[1] 153
ポジティブ(ネガティブ)率が 0.5 から離れるとサンプルサイズは多く必要

オンラインヘルプには思った以上に必要(有用)な情報が書いてあります。隅から隅まで,読むと良いでしょう。
例えば,デフォルトが twosided = FALSE であるのは,うっかり見逃すことも多いのではないかな?

Re: Rの'irr'パッケージ
  投稿者:0669SABON 2017/11/08(Wed) 06:00 No. 22472
青木先生

 基礎,内容も十分に理解していないまま使用方法をご質問し申し訳ございませんでした。
 先生の丁寧な解説にて使用方法は理解できましたが、やはり内容を十分に理解できていません。先生のご指摘の通り一般理論から勉強します。ご助言ありがとうございました。
 分からない場合はまたご質問させて頂きます。

Re: Rの'irr'パッケージ
  投稿者:0669SABON 2017/11/13(Mon) 05:45 No. 22473
再度、検定力分析について調べ直しましたがまだ十分に理解できておりません。
またirr パッケージのオンラインヘルプに引用されている文献ですが取り寄せることができておらず、まだ読んでおりません。
このような状況で質問するのは無礼だと承知していますが、ご教授いただけると幸甚です。

irrパッケージ のN2.cohen.kappaを使用したいのですが

オンラインヘルプの
mrg a vector of marginal probabilities given by raters
ということがわかりません。
そのため
Examples
require(lpSolve)
# Testing H0: kappa = 0.4 vs. HA: kappa > 0.4 (=0.6) given that
# Marginal Probabilities by two raters are (0.2, 0.25, 0.55).

0.2と0.25と0.55の数字が何を現しているかわかりません。
勉強不足で理解が及ばないためご教授頂けると大変助かります。
何卒よろしくお願いいたします。


R 計算結果オブジェクトからの要素の取り出し
  投稿者:明石 2017/10/29(Sun) 12:33 No. 22465
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

過日は、表の突合せ処理についてご教示いただき、誠にありがとうございました。
大変に良い勉強をさせていただきました。
改めて御礼を申し上げます。

Rについてお聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

-----------

質問は、計算結果オブジェクトからの要素の取り出しについてです。

Rの書籍に目を通しました、Google先生にもお聞きしましたが、
私の力では解決できないため、青木先生にご教示をいただきたいと思います。

NMFパッケージのマニュアルのexampleを例にご説明します。
#パッケージNMFの読み込み
library(NMF)

# random data
x <- rmatrix(20,10)

# run default algorithm with rank 2
res <- nmf(x, 2)

summary(res)

> summary(res)
rank sparseness.basis sparseness.coef silhouette.coef
2.0000000 0.1670822 0.4692352 1.0000000
silhouette.basis residuals niter cpu
1.0000000 13.6810699 440.0000000 NA
cpu.all nrun
NA 1.0000000
summary(res)から、以下の6つのスロットの値を取得したいと思います。
rank
sparseness.basis
sparseness.coef
silhouette.coef
silhouette.basis
residuals

resのクラスを確認しました。
> class(res)
[1] "NMFfit"
attr(,"package")
[1] "NMF"
NMFfitクラスということが分かりましたので、マニュアルを参照しました。
> ?NMFfit
http://127.0.0.1:10869/library/NMF/html/NMFfit-class.html
また、スロットの取得もしました。
> slotNames("NMFfit")
[1] "fit" "residuals" "method" "seed" "rng"
[6] "distance" "parameters" "runtime" "options" "extra"
[11] "call" "misc"
>
>
> getSlots("NMFfit")
fit residuals method seed
"NMF" "numeric" "character" "character"
rng distance parameters runtime
"ANY" ".functionSlotNULL" "list" "proc_time"
options extra call misc
"list" "list" "call" "list"
ここから先がまったく分かりません。

どのようにデータハンドリングすれば、
summary(res)の値を取り出すことができるのか、ご教示いただければ助かります。
rank
sparseness.basis
sparseness.coef
silhouette.coef
silhouette.basis
residuals

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

Re: R 計算結果オブジェクトからの要素の取り出し
  投稿者:青木繁伸 2017/10/29(Sun) 16:09 No. 22466
> a = summary(res)
> a
rank sparseness.basis sparseness.coef silhouette.coef silhouette.basis residuals
2.0000000 0.1547192 0.3082925 1.0000000 1.0000000 11.8603919
niter cpu cpu.all nrun
440.0000000 0.0940000 0.0940000 1.0000000
で,
> str(a)
Named num [1:10] 2 0.155 0.308 1 1 ...
- attr(*, "names")= chr [1:10] "rank" "sparseness.basis" "sparseness.coef" "silhouette.coef" ...
という名前付きのベクトルなので,
> a["rank"]
rank
2
> a["sparseness.basis"]
sparseness.basis
0.1547192
> a[c("residuals", "cpu")]
residuals cpu
11.86039 0.09400
のように,します。当然のことながら「何番目の要素」という指定もできます。
> a[2]
sparseness.basis
0.1547192
> a[10]
nrun
1
> a[7:9]
niter cpu cpu.all
440.000 0.094 0.094


【御礼】 Re: R 計算結果オブジェクトからの要素の取り出し
  投稿者:明石 2017/10/29(Sun) 20:37 No. 22467
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

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

ご丁寧に、わかりやすく、ご説明をしてくださいましたので、よくわかりました。

先に進めないで苦慮していましたが、仕組みが分かりました。

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


傾向スコアについて
  投稿者:課長 2017/10/17(Tue) 19:08 No. 22454
いつもお世話になっています。
独立変数が間隔尺度で、共変量がいくつかある2群について、傾向スコアで2群をペアにして対応のあるt検定をしようと思うのですが、この場合、マッチングは恣意的なものになるのではないかと思います。また、ペアの作り方で結果が変わると思うのですが、実際どうでしょう。

Re: 傾向スコアについて
  投稿者:青木繁伸 2017/10/19(Thu) 22:32 No. 22455
恣意的だし,一意には決まらないでしょう。
そもそも,群ごとのサンプルサイズが異なるとどうしようもないし。

Re: 傾向スコアについて
  投稿者:課長 2017/10/24(Tue) 20:58 No. 22464
ご返信ありがとうございます。返信が遅れて申し訳ありません。
やはり恣意的になるわけですよね。なんというか、回帰分析のP値については多重性を考慮しないとか、マッチングは恣意的とか、どうも医学系統計というか疫学系統計は私の頭では理解できない部分が多くて情けないです。
どうもありがとうございました。


R 表の突合せ集計
  投稿者:明石 2017/10/22(Sun) 18:08 No. 22459
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

昨日は、ネットワーク分析igraphパッケージの計算結果の取り出し、組み立てについて
ご教示をいただきまして、誠にありがとうございました。

改めて御礼を申し上げます。

また、Rについてお聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

-----------

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

表1と表2を突き合わせて、表3を出力します。

表1と表2はともに大規模データであることから、効率的な方法をご教示いただければ
大変に助かります。

表1は,wordとgroupの関係を示しています。 マスター表の位置づけです。
AとBは、g1に、
CとDは、g2に
EとFとGは、g3に属していることを示しています。

表2は、3つの文書があり、
各文書にwordが含まれるかどうかを、0/1の2値で示しています。

表2には、表1には含まれないwordも出現します。
(ここでは、H、P、X)

表1と表2を突き合わせて、
表2の各文書に出現するwordを、表1を参照して、wordが属するgroup単位に度数を集計して、表3を出力します。

表1と表2はともに大規模データであることから、効率的な方法をご教示いただければ
大変に助かります。

お手数をおかけいたします。
どうぞ、よろしくお願いいたします。


Re: R 表の突合せ集計
  投稿者:青木繁伸 2017/10/22(Sun) 21:36 No. 22460
仮想的な大規模データを生成するプログラムと,その大規模データを処理する「標準的」プログラムを提示いただければ,その「標準的」プログラムよりどの程度優秀なプログラムであるかを提案できると思います。こちらが思い浮かべたプログラムが「標準的」プログラムより劣っていれば,申し訳ないですから。

Re: R 表の突合せ集計
  投稿者:青木繁伸 2017/10/22(Sun) 23:05 No. 22461
表1のサイズが 20×2
表2のサイズが 1000000×26
のときの実行時間が 0.907 秒でした
example = FALSE で,あなたの提示したテストデータでの結果が表示されます
system.time({
set.seed(12345)
exampe = TRUE
if (exampe) {
l = 20
g = sort(sample(5, l, replace=TRUE))
tbl1 = data.frame(word=LETTERS[1:20], group=sprintf("g%02i", g))
o = order(tbl1[,2])
tbl1 = tbl1[o,]
} else {
l = 7
tbl1 = data.frame(word=LETTERS[1:7], group=paste("g", c(1,1,2,2,3,3,3), sep=""))
}
print(tbl1)
print(dim(tbl1))

if (exampe) {
m = 1000000
n = 26
tbl2 = cbind(1:m, matrix(sample(0:1, m*n, replace=TRUE), m, n))
colnames(tbl2) = c("ban", LETTERS)
head(tbl2, 20)
} else {
m = 3
n = 10
tbl2 = matrix(c(1,1,0,1,0,0,0,1,0,0,1,1,1,0,1,0,0,1,0,1,1,1,0,0,0,1,0,0,0,1), 3, 10)
dimnames(tbl2) = list(1:3, c(LETTERS[1:7], "H", "P", "X"))
}
head(tbl2)
print(dim(tbl2))

tbl = xtabs(~word+group, tbl1)
ans = matrix(0, m, ncol(tbl)+1)
for (i in 1:ncol(tbl)+1) {
name = rownames(tbl)[tbl[,i-1] == 1]
ans[, i] = rowSums(tbl2[,name])
}
ans[,1] = tbl2[,1]
colnames(ans) = c("ban", colnames(tbl))
print(head(ans))
})
おそらく,もっとシンプルな for ループを使ったプログラムでも,そこそこの実行時間ではないかと思います。

Re: R 表の突合せ集計
  投稿者:青木繁伸 2017/10/22(Sun) 23:28 No. 22462
おそらく,最も速いのは,データ構造を反映した
cibnd(rowSums(tbl2[,1:2]), rowSums(tbl2[,3:4]), rowSums(tbl2[,5:7])) などでしょうけど

【御礼】Re: R 表の突合せ集計
  投稿者:明石 2017/10/23(Mon) 06:10 No. 22463
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

色々な方法をご教示くださり、大変に良い勉強をさせていただきました。

考え方をしっかりと勉強をいたします。

ご丁寧に、ご親身に、ご教示をいただき、誠にありがとうございました。


R igraphパッケージのクラスタ計算結果の取り出し
  投稿者:明石 2017/10/21(Sat) 15:58 No. 22456
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rについてお聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

-----------

ネットワーク分析の igraph パッケージのクラスタ計算結果の取り出しに苦慮しています。

簡単なネットワークデータを、以下にお示しをします。

node1 node2 weight
2 5 0.03
1 3 0.04
1 4 0.08
3 5 0.20
5 6 0.15
4 7 0.25

このネットワークデータを作成するRスクリプト例です。

edge.list <- matrix( c(2, 5, 1, 3, 1, 4, 3, 5, 5, 6, 4, 7), nrow=6, ncol=2, byrow=T)

weight <- c(0.03, 0.04, 0.08, 0.20, 0.15, 0.25)

library(igraph)

g <- graph_from_edgelist(edge.list, directed=FALSE)

E(g)$weight <- weight

plot(g)

グラフクラスタリングをします。

cluster_optimal(g)

> cluster_optimal(g)
IGRAPH clustering optimal, groups: 2, mod: 0.46
+ groups:
$`1`
[1] 1 4 7

$`2`
[1] 2 3 5 6

7つのノードは、2つのグループにクラスタリングされることが分かりました。
この結果を、前記プロット図に追記した画像を、添付ファイルにてお示しをします。

以下のQ1、Q2の2つについて、Rプログラムをお示しいただけると、助かります

(Q1)
2つのグループがあり、
1番目のグループには、ノードの 1 4 7
2番目のグループには、ノードの 2 3 5 6
が属することを、どのように取り出せばよいのでしょうか

(Q2)
Q1ができないことから、Q2もできませんが、
クラスタリング結果を以下のデータフレームで整理したいのですが、どのように組み立てすればよいのでしょうか。

node group
1 1
2 2
3 2
4 1
5 2
6 2
7 1

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


Re: R igraphパッケージのクラスタ計算結果の取り出し
  投稿者:青木繁伸 2017/10/21(Sat) 19:27 No. 22457
この分析法については知りませんので,正しい解答か分かりませんが,以下のようにすればお望みの結果になるようですが?
> a = cluster_optimal(g)
> group = a$membership
> node = seq_along(group)
> d = data.frame(node, group)
> d
node group
1 1 1
2 2 2
3 3 2
4 4 1
5 5 2
6 6 2
7 7 1

Re: R igraphパッケージのクラスタ計算結果の取り出し
  投稿者:明石 2017/10/21(Sat) 20:50 No. 22458
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

ありがたいご教示をいただきまして、誠にありがとうございます。
助かりました。

> a = cluster_optimal(g)
> str(a)
List of 2
$ modularity: int [1:3] 1 4 7
$ membership: int [1:4] 2 3 5 6
- attr(*, "class")= chr "communities"
> a$membership
[1] 1 2 2 1 2 2 1

このmembershipを使うのですね。

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


母集団から抽出した複数の値の平均の下側確率
  投稿者:tosh 2017/10/05(Thu) 21:41 No. 22446
統計初心者です。

ある母集団から無作為抽出した、
サイズ 29、 平均 155.7、 SD 5.6 の標本があります。
ほぼ正規分布しています。

平均 - 2SD = 144.5 なので、
母集団から新たに1値を取り出した場合に、
それが 144.5より小さい値である確率は、
標準正規分布表の片側確率から約2.3% になると思います。

たとえば今後、母集団から新たに5値を取り出すことを考える場合、
その5値の平均が144.5より小さい値である確率も、
やはり約2.3%になると言えるでしょうか?

平均となると、取り出した値の大小によって変わるので、
単純に144.5より小さいか大きいか以外にも
考えるべきことがないのだろうか、と腑に落ちません。

Re: 母集団から抽出した複数の値の平均の下側確率
  投稿者:まさき君 2017/10/06(Fri) 08:18 No. 22447
正規分布の元である中心極限定理からそうなります。
「中心極限定理」について調べてみれば納得できるかと。

Re: 母集団から抽出した複数の値の平均の下側確率
  投稿者:tosh 2017/10/06(Fri) 09:38 No. 22448
アドバイスありがとうございます。
中心極限定理を調べ、
「標本平均の分布は、
 サイズを大きくするほど正規分布に近づく」
と理解しました。
おかげさまで、疑問点が明確になりました。

たとえば、新たにとる値が1000個で、
1000というのが十分「大きい」と言って良さそうなら、
取り出された1000個の値の平均(標本平均)はほぼ正規分布する、
と考えてよいことになり、
「標本平均が-2SD点(144.5)より小さい値になる確率は2.3%」…(a)
と言ってよいのかな、と思いました。

しかし今回、新たにとる値が5個のように「小さい」場合は、
正規分布とは異なる分布をしていることになるので、
(a)の主張はできないのではないか、という
もやもや感が残るのです。

Re: 母集団から抽出した複数の値の平均の下側確率
  投稿者:qq 2017/10/06(Fri) 18:37 No. 22450
1000個の標本を取り出して得られた平均は、何回やっても2SDの範囲に全て入ってしまうだろうな思います。
仮に母集団全体を測定できるとすると、何回数えても母集団の平均しか出てこないので、「標本平均が-2SD点(144.5)より小さい値になる確率」はゼロ%になるんじゃない?
5つの標本の平均、1000個の標本の平均、無限個の標本の平均を考える時、それぞれ異なる結果になると思うね。
最初の「ある母集団」と、そこからn個の標本を取り出し得られた平均からなる母集団は、別物になると思います。

Re: 母集団から抽出した複数の値の平均の下側確率
  投稿者:青木繁伸 2017/10/07(Sat) 20:18 No. 22451
母平均 = μ = 155.7,母標準偏差 = σ = 5.6 から,サンプルサイズ n の標本を取り出したとき,その標本平均は理論的に \bar{x} = μ, 平均の標準誤差は σ/ \sqrt{n}

シミュレーションの試行回数を trial = 1000000 として,

(1) サンプルサイズ = 1

> 母集団から新たに1値を取り出した場合に、
> それが 144.5より小さい値である確率は、
> 標準正規分布表の片側確率から約2.3% になると思います。

> x = rnorm(trial, mean=155.7, sd=5.6) # trial 回シミュレーション
> mean(x < 144.5) # mean の使い方に戸惑うかも知れないが, sum(s < 144.5)/trial のこと
[1] 0.022514

はい,正解です。

(2) サンプルサイズ = 5

> たとえば今後、母集団から新たに5値を取り出すことを考える場合、
> その5値の平均が144.5より小さい値である確率も、
> やはり約2.3%になると言えるでしょうか?

> x = matrix(rnorm(5*trial, mean=155.7, sd=5.6), 5, trial) # 5 行 trial 列のサンプルデータ
> y = apply(x, 2, mean) # trial 回ごとの平均値
> mean(y < 144.5) # その平均値が 144.5 未満の確率
[1] 3e-06 # おー!とっても小さい

なぜ?

> mean(y)
[1] 155.7009 # 平均値の平均値は母平均に近く,かつ,
> sd(y)
[1] 2.504464 # そのばらつきはとっても小さい!!平均が 144.5 以下になることはほとんどない。

だからこそ,標本平均で母平均を推定できるんです。

Re: 母集団から抽出した複数の値の平均の下側確率
  投稿者:tosh 2017/10/11(Wed) 13:41 No. 22453
青木先生、はじめまして。
ご指導ありがとうございます。

r をインストールして、ご教示いただいたコードを試してみました。
このような方法で、実際に大標本を作って、検証ができるのですね。
r scriptは初めてですが、記述がとてもシンプルで驚きました。

(確かに、mean(x < 144.5) という書き方は、シンプル過ぎて戸惑います)

最初の母集団Aから、(5値を取り出し、その平均値をとる)
という試行を、実際に100万回してみて、
100万個の要素からなる集団Bを作った。
集団Bのばらつきは、母集団Aのばらつきよりもずっと小さかった。

集団Bは、母集団Aよりも、とんがった形をしていて、
母平均を鋭敏に指し示す形になっているのが想像できました。

まさき君さん、qqさんもありがとうございました。
私のNo.22448での理解はかなり誤っていましたね。

もう少し、色々試してみます。


cox回帰分析の因子
  投稿者:jj 2017/09/07(Thu) 06:48 No. 22436
いつも勉強させていただいております。
非正規分布するある因子に関する研究で、正規分布とするために対数変換した際に1以下の因子は負の値となりますが、負の値を含んだデータでもcox回帰分析の説明変数として計算可能でしょうか?

Re: cox回帰分析の因子
  投稿者:青木繁伸 2017/09/07(Thu) 08:25 No. 22437
計算可能かどうかはやってみればすぐにわかること

それより,正規分布しない変数を分析に使えないとなると,どのように変数変換しても正規分布にならないダミー変数はどうすればよいのですか?

Re: cox回帰分析の因子
  投稿者:jj 2017/09/08(Fri) 03:16 No. 22438
お返事いただきましてありがとうございます。質問の仕方が不適格で失礼致しました。

計算はできました。負の値を含んだデータを説明変数に組み込むことは、解析方法として誤りではありませんでしょうか。

>それより,正規分布しない変数を分析に使えないとなると,どのように変数変換しても正規分布にならないダミー変数はどうすればよいのですか?
その場合は気にせずに分析するしかないのでしょうか?先生のお考えはございますでしょうか?

Re: cox回帰分析の因子
  投稿者:青木繁伸 2017/09/10(Sun) 21:06 No. 22439
> その場合は気にせずに分析するしかないのでしょうか?

あなたの信念で,「正規分布しない変数は分析に使えない」(だから,対数変換する)
ということなら,「ダミー変数は使わない」ということでしょ?

単純なことです。

Re: cox回帰分析の因子
  投稿者:jj 2017/09/12(Tue) 06:31 No. 22440
ご返答いただきましてありがとうございました。
大変参考になりました。また勉強させていただきます。

Re: cox回帰分析の因子
  投稿者:scdent 2017/09/13(Wed) 12:39 No. 22441
もう解決済みだとは思いますが、
説明変数の値に負の数値があっても重回帰分析に用いることができます。

このことは、
\睫席竸瑤涼佑防蕕凌値がない場合の解析結果と、
⊂紊寮睫席竸瑤ら定数(平均でも可)を引いたものを説明変数として置き換えた場合の解析結果を比較すれば、一目瞭然です。
´△修譴召譴寮睫席竸瑤硫鶺係数は同じになるはずです(切片は異なりますよ)。

つまり、重回帰分析は説明変数の値の1増分あたりの効果を、回帰係数として求めているので、負の値でも正の値でも使えるというわけです。

Re: cox回帰分析の因子
  投稿者:jj 2017/10/10(Tue) 15:25 No. 22452
scdent様

ご丁寧にお教えいただきまして大変ありがとうございました。


効果量と標準化応答平均の使い方
  投稿者:Mori 2017/10/06(Fri) 09:59 No. 22449
介入前後での治療効果を比較するためのResponsivenessの指標として、Cohenの「効果量(Effect Size: SE)」や「標準化応答平均(Standardized Response Mean: SRM)」がありますが、効果量にはパラメトリック(ES=Mean/SDBaseline)とノンパラメトリック(ES = Z/sqrt(N))があると記載されています。SRM = ME/SDDifferenceはこの式のみです。

しかし、論文を見渡すと、順序尺度であろうかなかろうがSRMは上記の式のみで計算されており、SEも順序尺度のデータであっても、パラメトリックの式だったり、ノンパラメトリックの式だったりしています。

教えていただきたいのですが、例えば、順序尺度のデータの場合、
 SEはパラメトリックまたはノンパラメトリックのどちらを使うのが正しいのでしょうか?
SRMの場合はパラメトリックやノンパラメトリックということを気にしなくてもいいのでしょうか。SRMのノンパラメトリックという式が見たことがありません。

どのようにして使い分ければ良いのか理解できず混乱しております。


コクラン・アーミテージ検定でのセルごとの数
  投稿者:木村 2017/09/19(Tue) 20:32 No. 22442
コクラン・アーミテージ検定で濃度を変えて症状(死亡など)が発生したデータを検定する時。
ある濃度(群)でたまたま極端な外れ値が出て、低濃度で発症した数が低いような場合もあります。

いろいろ頑張って調べたところ。
カイ二乗分布を使っているから発症数が5以上でなければ使えない(発症数が5になるまである濃度でデータ蓄積しなさい)または他の統計手法を使いなさいといったガイドを余所にて見つけました。
セルごとの値はあまり関係無いという海外フォーラムの回答も見つけましたけど、論拠がちゃんと説明されておりません。

サンプルデータとして低濃度で発症数が5未満のもの、または高濃度で未発症数が5未満なものを使ってるサイトも多くみられます。
(たとえば青木さまのサイトでしたら、http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Armitage.html にて独立変数10 ケース数30 陽性数2のデータセットで解説されております)

セル内で5未満の値になるデータで、コクラン・アーミテージ検定を使っていいものでしょうか?
第一種過誤はデータの総数が少なかったとして扱うつもりです。

Re: コクラン・アーミテージ検定でのセルごとの数
  投稿者:青木繁伸 2017/09/19(Tue) 21:36 No. 22443
> カイ二乗分布を使っているから発症数が5以上でなければ使えない

これは,独立性のカイ二乗検定から来ているのかも知れませんが,それならば,「発症数」ではなく「発症数の期待値」のはずですが。

そのように書かれているページの URL を教えてください。

Re: コクラン・アーミテージ検定でのセルごとの数
  投稿者:木村 2017/09/20(Wed) 18:41 No. 22444
知人に相談したところ先日そのような趣旨でアドバイスされ、必死にぐぐって「セル内の数字が5未満では」と記載されたサイトを5つ見つけまして、いま必死に検索したものの全部は見つかりませんでした申し訳ありません。
ただただ、わが身の浅学と理解の浅さを恥じ入るばかりです。

医学統計勉強会 - 東北大学病院 - 循環器内科
>Cochran-Armitage 検定はχ2 分布を用いた近似検定であるので,セルの数が小さいときは他の性格検定を用いる必要があります.
http://www.cardio.med.tohoku.ac.jp/newmember/pdf/ms/16_5T.pdf

ちょっと違いますけど
Is the Cochran Armitage Trend Test valid, if any of the counts are less than 5?
https://www.researchgate.net/post/Is_the_Cochran_Armitage_Trend_Test_valid_if_any_of_the_counts_are_less_than_5

Re: コクラン・アーミテージ検定でのセルごとの数
  投稿者:木村 2017/09/27(Wed) 20:11 No. 22445
こちら以外にも友人に相談したりしましたけれど、駄目とも良いともいえず。

色々と考えて結局、これをもって後々文句つけられたりした際に、これでは私どもは批判者に対して適切な対応できませんし。
「セル内で5未満の値になるデータ」は追試して全部5以上になるように満たしてからまたコクラン・アーミテージ検定を用いる事といたしました。

解決に至りませんでしたけど、ご返信ありがとうございました。


信頼性の検定におけるサンプルサイズ
  投稿者:0669SAVON 2017/07/21(Fri) 05:39 No. 22430
いつも拝見させて頂いています。
ご質問が2点あります。
信頼性の検定を行っていますがサンプルサイズが分かりません。

2種類の評価方法の信頼性の検定を別々に行っています。

一つ目は7段階の指標で先行研究では中央値5(四分位範囲3-6)です。
12名の検査者にてFleissのκ統計量を用い検定します。

2つ目は先行研究では平均値15(標準偏差4)の評価方法です。
12名の検査者にてICC(2.1)を用い検定します。

これらのサンプルサイズの出し方をご教授頂けると幸甚です。

Re: 信頼性の検定におけるサンプルサイズ
  投稿者:青木繁伸 2017/07/21(Fri) 09:34 No. 22431
なんのために「サンプルサイズの出し方」が必要なのですか?
このような場合に「サンプルサイズ」という用語は適切でしょうか?
Number of subjects
Number of Judges(Raters)
のどちらが「サンプル」でしょうか?

Re: 信頼性の検定におけるサンプルサイズ
  投稿者:0669SAVON 2017/07/21(Fri) 11:01 No. 22432
用語を正しく説明できず申しわけありません。
Numver of subjectsです。
κ統計量が0.25と低値であったため信頼性が低いのではと考えましたが、査読者より検出力が低いためではとご指摘を頂ました。
今回の検定における12名の症例に対し、12名の検査者で信頼性の検定を行う際の検出力を求めることが本来の目的です。

Re: 信頼性の検定におけるサンプルサイズ
  投稿者:青木繁伸 2017/07/21(Fri) 13:53 No. 22433
「必要な検出力を得るためのサンプルサイズ」,「特定のサンプルサイズにおける検出力」の計算ということですね。

Cohen のκ統計量についてはよく知られているようですが,
Fleiss のκ統計量についてはなかなか見つからないですね。

icc については,ICC.Sample.Size パッケージにも関数が用意されているようです。
calculateIccPower
   Function to calculate post-hoc power for ICC studies
calculateIccSampleSize
   Function to calculate sample size required for
   studies where ICC is primary outcome.
その reference には
Zou, G. Y. (2012). Sample size formulas for estimating intraclass correlation coefficients with precision and assurance. Statistics in medicine, 31(29), 3972-3981.
が挙げられています。

Re: 信頼性の検定におけるサンプルサイズ
  投稿者:0669SAVON 2017/07/21(Fri) 17:37 No. 22434
青木先生
 
 情報ありがとうございます。
 早速、論文を確認します。

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


Rstudioのエラー関する質問です
  投稿者:y931 2017/07/17(Mon) 17:44 No. 22425
はじめまして
いつも拝見させていただいております.
統計学の質問ではありませんので、この掲示板に相応しくない場合には削除ください.

今回、Rstudioの利用でで非常に困った状況になり、投稿させていただきました.
もしご存知でしたら対処方法をお知らせください.
宜しくお願いいたします

RStudio Version 1.0.136
R Version 3.3.3

RStudio で次のようなデータフレーム作成した場合
A <- c ( 12, 10, 13 )
B <- c ( 20, 12, 16 )
C <- c ( 20, 15, 19 )
y <- data.frame (A,B,C)

以前はグラフのタブが表れて、エクセルのようなクロス表が表示されていました.
しかし最近、その表内の行(1 , 2 , 3)と列(A , B , C)のラベルが表示されず、
ラベルセルが空白になるエラーが発生しています.
(添付ファイル)
使用しているOSは、windows10 home, 64bit です.C)

なおRjpWikiにも同じ質問を投稿しております.


Re: Rstudioのエラー関する質問です
  投稿者:y931 2017/07/20(Thu) 11:13 No. 22429
Rstudioのサポートページ(https://support.rstudio.com/hc/en-us)に質問して解決しました.
バージョンを更新 (v1.0.153)することで不具合が解消されました.
お騒がせいたしました.


多重ロジスティック回帰分析
  投稿者:初学者 2017/07/18(Tue) 14:41 No. 22426
いつも大変お世話になっております。
多重ロジスティック回帰分析で、単変量のロジスティック回帰分析と同じように、カットオフ値を求めることはできるのでしょうか。ご教示の程、よろしくお願い申し上げます。

Re: 多重ロジスティック回帰分析
  投稿者:青木繁伸 2017/07/18(Tue) 15:26 No. 22427
その点については,多重ロジスティック回帰分析も単変量のロジスティック回帰分析も同じでしょう

Re: 多重ロジスティック回帰分析
  投稿者:初学者 2017/07/19(Wed) 14:25 No. 22428
大変お世話になっております。
お返事をいただき、ありがとうございました。

ロジスティック回帰分析について、あまりにも不勉強でした。
また、過去ログのほうでロジスティック回帰について検索し、
http://aoki2.si.gunma-u.ac.jp/lecture/mb-arc/arc041/05383.html
の記事で、解決することができました。
確認が不足しており、申し訳ありません。
今後とも何卒よろしくお願いいたします。


ポアソン回帰におけるパラメータ推定
  投稿者:kamo 2017/07/09(Sun) 16:15 No. 22415
お世話になります。
Rのoptim関数を用いて、ポアソン回帰(説明変数2つ)の係数を推定したいのですが、うまくいきません。

> x1 <- as.factor(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
> x2 <- as.factor(c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50))
> y <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1)
>
>
> logL <- function(par){
+ beta0 <- par[1]
+ beta1 <- par[2]
+ beta2 <- par[3]
+ sum(dpois(y,exp(beta0+beta1*x1+beta2*x2),log = TRUE))
+ }
>
> opt <- optim(c(1,1,1),logL,control = list(fnscale=-1),method = "BFGS",hessian = TRUE)
Error in optim(c(1, 1, 1), logL, control = list(fnscale = -1), method = "BFGS", :
initial value in 'vmmin' is not finite
In addition: Warning messages:
1: In Ops.factor(beta1, x1) : * not meaningful for factors
2: In Ops.factor(beta2, x2) : * not meaningful for factors
>

ご教示頂ければ幸いです。よろしくお願い致します。

Re: ポアソン回帰におけるパラメータ推定
  投稿者:青木繁伸 2017/07/09(Sun) 18:15 No. 22416
> x1 <- as.factor(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
> x2 <- as.factor(c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50))
x1, x2 が factor なので,その数値演算は無理です。

任意の beta0, beta1, beta2 とすると,
> beta0 = 0
> beta1 = -0.1
> beta2 = -0.11
> beta0+beta1*x1+beta2*x2
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ は因子に対しては無意味です
2: Ops.factor(beta2, x2) で: ‘*’ は因子に対しては無意味です
3: Ops.factor(beta1, x1) で: ‘*’ は因子に対しては無意味です
4: Ops.factor(beta2, x2) で: ‘*’ は因子に対しては無意味です
> exp(beta0+beta1*x1+beta2*x2)
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ は因子に対しては無意味です
2: Ops.factor(beta2, x2) で: ‘*’ は因子に対しては無意味です
> dpois(y,exp(beta0+beta1*x1+beta2*x2),log = TRUE)
[1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ not meaningful for factors
2: Ops.factor(beta2, x2) で: ‘*’ not meaningful for factors
> sum(dpois(y,exp(beta0+beta1*x1+beta2*x2),log = TRUE))
[1] NA
警告メッセージ:
1: Ops.factor(beta1, x1) で: ‘*’ not meaningful for factors
2: Ops.factor(beta2, x2) で: ‘*’ not meaningful for factors
となるので,optim では求まりません。

独立変数が factor のとき,lm や glm などでは自動的にダミー変数が作成されるのですが,今回のような場合は自分でダミー変数を作って対処しないとだめでしょう。
> # x1 は既にダミー変数
> x1 = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)
> x2 = c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50)
> x2.50 = (x2 == 50)+0 # かっこは必要, +0 で logical を numeric に
> x2.60 = (x2 == 60)+0
> # x2 == 40 の場合は,x2.50 = 0, x2.60 = 0
>
> logL <- function(par){
+ sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE)) # par を直接使って問題ない
+ }
> opt <- optim(rep(1, 4),logL,control = list(fnscale=-1),method = "BFGS",hessian = TRUE)
> opt
$par
[1] -0.6495380 -1.4754486 0.9360592 0.3696557

$value
[1] -38.44917

$counts
function gradient
32 19

$convergence
[1] 0

$message
NULL

$hessian
[,1] [,2] [,3] [,4]
[1,] -28.998906 -4.0006148 -1.799968e+01 -9.999199e+00
[2,] -4.000615 -4.0006148 -3.350026e+00 -1.728501e-01
[3,] -17.999681 -3.3500262 -1.799968e+01 1.776357e-09
[4,] -9.999199 -0.1728501 1.776357e-09 -9.999199e+00
> # 検証
> par = opt$par
> exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60)
[1] 1.3317864 1.3317864 1.3317864 0.7558728 0.7558728 0.7558728 0.7558728 0.5222870 1.3317864 1.3317864 1.3317864 0.7558728 0.7558728
[14] 0.7558728 0.7558728 0.7558728 0.7558728 0.7558728 0.7558728 1.3317864 1.3317864 1.3317864 1.3317864 0.7558728 1.3317864 0.1194346
[27] 0.1194346 0.1194346 0.1194346 0.3045477 0.3045477 0.3045477 0.3045477 0.3045477 0.3045477 0.3045477 0.1728500 0.3045477 0.3045477
[40] 0.3045477 0.3045477
> dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE)
[1] -1.3317864 -1.3317864 -1.3317864 -0.7558728 -0.7558728 -0.7558728 -0.7558728 -1.1718250 -1.0452652 -1.0452652 -1.0452652 -1.0357550
[13] -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.0357550 -1.4518912 -1.4518912 -1.4518912 -1.4518912 -2.0087844
[25] -2.2639823 -0.1194346 -0.1194346 -0.1194346 -0.1194346 -0.3045477 -0.3045477 -0.3045477 -0.3045477 -0.3045477 -0.3045477 -0.3045477
[37] -0.1728500 -1.4934752 -1.4934752 -1.4934752 -1.4934752
> sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE))
[1] -38.44917

Re: ポアソン回帰におけるパラメータ推定
  投稿者:kamo 2017/07/09(Sun) 22:03 No. 22417
詳細なご回答、誠にありがとうございます。
ダミー変数・デザイン行列まで作成してから、関数に放り込むべきであったのですね。
とてもよく理解できました。
ありがとうございました。

Re: ポアソン回帰におけるパラメータ推定
  投稿者:kamo 2017/07/10(Mon) 15:06 No. 22418
度々失礼します。
上記質問の続きとしまして、
「残差逸脱度Residual Deviance」の算出方法を知りたいです。

ご回答の最終行
> sum(dpois(y,exp(par[1]+par[2]*x1+par[3]*x2.50+par[4]*x2.60),log = TRUE))
[1] -38.44917
がこのモデルにおける対数尤度と理解しておりますが、
ここから残差逸脱度というものが計算できますでしょうか?

ご教示頂ければ幸いです。

Re: ポアソン回帰におけるパラメータ推定
  投稿者:青木繁伸 2017/07/11(Tue) 15:04 No. 22419
「対数尤度の -2 倍がデビアンス」では?

Re: ポアソン回帰におけるパラメータ推定
  投稿者:kamo 2017/07/15(Sat) 15:17 No. 22420
書籍等で調べると、先生のおっしゃるように残差逸脱度=対数尤度の-2倍との記述がございました。
しかし、

> x1 <- as.factor(c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))
> x2 <- as.factor(c(50,50,50,60,60,60,60,40,50,50,50,60,60,60,60,60,60,60,60,50,50,50,50,60,50,40,40,40,40,50,50,50,50,50,50,50,60,50,50,50,50))
> y <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,3,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1)
> fm <- glm(y~x1+x2,family="poisson"(link="log"))
> summary(fm)

Call:
glm(formula = y ~ x1 + x2, family = poisson(link = "log"))

Deviance Residuals:
Min 1Q Median 3Q Max
-1.6321 -0.7804 -0.3009 0.5385 1.2393

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.6495 1.0356 -0.627 0.53058
x11 -1.4757 0.5638 -2.617 0.00886 **
x250 0.9360 1.0405 0.900 0.36832
x260 0.3697 1.0805 0.342 0.73225
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

Null deviance: 40.539 on 40 degrees of freedom
Residual deviance: 28.838 on 37 degrees of freedom
AIC: 84.898

Number of Fisher Scoring iterations: 5

Residual deviance: 28.838と出力され、(これが残差逸脱度ですよね?)
-2*(-38.44917) = 76.89833
と一致しません。
どこかで誤解しているでしょうか?
ご教示いただければ幸いです。

Re: ポアソン回帰におけるパラメータ推定
  投稿者:青木繁伸 2017/07/15(Sat) 21:55 No. 22421
glm でできるなら,わざわざ変なことして optim でやる必要ないでしょう

full model の 対数尤度は
> sum(dpois(y, lambda=y, log=TRUE)
[1] -24.03019

したがって,その -2 倍は full model の deviance
> sum(dpois(y, lambda=y, log=TRUE))*-2
[1] 48.06037

当該モデルの residual deviance が 28.838 なので,deviance は
> sum(dpois(y, lambda=y, log=TRUE))*-2+28.838
[1] 76.89837

逆にたどれば,当該モデルの deviance が 76.89837 で full model の deviance が48.06037 なので,residual deviance は
76.89837 - 48.06037 = 28.838

Re: ポアソン回帰におけるパラメータ推定
  投稿者:kamo 2017/07/16(Sun) 19:07 No. 22422
青木先生
ご回答ありがとうございます。
glm関数におけるポアソン回帰の計算手順を理解することが、このプログラミングの目的でした。お手数おかけして、すみません。
とてもよく理解することができました。ありがとうございました。

Re: ポアソン回帰におけるパラメータ推定
  投稿者:青木繁伸 2017/07/16(Sun) 19:14 No. 22423
> glm関数におけるポアソン回帰の計算手順を理解すること

そのような場合には,ソースコードを読むことをお勧めします
一点の曖昧さもなく,完璧に記述されていますから

Re: ポアソン回帰におけるパラメータ推定
  投稿者:kamo 2017/07/17(Mon) 15:02 No. 22424
ご返信ありがとうございます。
ソースコードを読解するのは中々ハードルが高そうで、敬遠していましたが...
やはり、計算手順の理解には、その方法が確実ということですよね。
チャレンジしてみます。
またつまづくようなことがあったときは、この掲示板にお世話になるかと思いますが、今後ともご指導よろしくお願い致します。


相関分析と単回帰分析と決定係数
  投稿者:いし 2017/07/08(Sat) 17:35 No. 22412
はじめて質問させていただきます。
同一人物が実施した2種類のテストの関係性を観たく、
相関、および単回帰をかけました。
その結果、回帰のほうでは有意差が出ませんでした。
この場合は、相関係数 R も R^2の値にかかわらず、
有意差なし、という点を優先して、結果の解釈をしたほうが良いでしょうか。

おそれいりますが、おしえていただけませんか?

Re: 相関分析と単回帰分析と決定係数
  投稿者:青木繁伸 2017/07/09(Sun) 08:47 No. 22414
相関係数の検定も,単回帰分析の係数(傾き)の検定も,決定係数の検定も,p 値は同じです


計測不能の検査の統計
  投稿者:てんぽ 2017/07/05(Wed) 15:36 No. 22408
青木先生

お世話になります。

ある疾患を診断するための検査について調べています。
神経を伝わるスピードを計測する検査で、
刺激してから、刺激により筋肉が収縮するまでの時間を計測します。
通常であれば4.5ミリ秒以下で、異常となると検査値が上昇し、
重症例では計測できなくなります。
例えば、
  3.8ms…正常
  5.2ms…軽症
  7.8ms…重症
  計測できず…最重症
などです。

この検査値と、他のパラメーター(年齢や他の検査値など)との
関連を調べようとしています。
その際、計測できればその計測値を用いれば良いのですが、
最重症の、計測できなかった症例の値をどのようにして統計するべきでしょうか。

なにか大きな数字を当てはめるのか、
そもそも統計はできないのか、困っています。
ご教示のほどよろしくお願いします。

Re: 計測不能の検査の統計
  投稿者:青木繁伸 2017/07/05(Wed) 23:21 No. 22409
3.8,5.2,7.8 も直線関係ではないので,順序尺度として扱うのが妥当かとおもいます。つまり,正常,軽症,重症,最重症の4段階ですね。
関連を見る他の変数が比尺度,間隔尺度,順序尺度の場合は順序相関係数,名義尺度変数の場合は,関連性係数(クラメール係数など)
関連を見る他の変数が比尺度,間隔尺度の場合は,正常,軽症,重症,最重症ごとの平均値,中央値などの比較ということになるかと。
また,多変量解析の場合は,正常,軽症,重症,最重症をダミー変数として用いれば普通の重回帰分析,ロジスティック回帰,判別分析などを適用できます。

Re: 計測不能の検査の統計
  投稿者:てんぽ 2017/07/08(Sat) 08:30 No. 22411
どうもありがとうございます。
重症度を定義して分類するべきなんですね。
その手法で検討させていただきます。

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


直線性の適合度検定
  投稿者:やまが 2017/06/17(Sat) 09:43 No. 22396
青木先生

突然のメール失礼いたします。
私は現在、発がん性物質のMedian Toxic Dose (TD50)を、
Carcinogenic Potency Database (CPDB, https://toxnet.nlm.nih.gov/cpdb/)に
記載の方法に従って計算しようとしています。
⇒計算方法(https://toxnet.nlm.nih.gov/cpdb/td50.html)

投与量[mg/kg-bw]: 0, 232, 464, 927
腫瘍発生確率[-] :0.08, 0.62, 0.84, 0.92

この方法では、まず、物質の投与量に対する腫瘍発生率の変化の直線性について
適合度検定を行っています。
a χ2 goodness-of-fit statistic tests the validity of the linear relationship between dose and tumor incidence

この直線性の適合度検定をRで行いたいと考えています。
以下のページなども参考にさせていただき、
http://www.tamagaki.com/math/Statistics606.html

帰無仮説:”理論値(線形式から得られる値)と比較してずれがない。”
対立仮説:”理論値(線形式から得られる値)と比較してずれがある。”

とすればよいのかな?と思ったのですが、具体的な計算方法がわかりません。
御教示いただけますと大変ありがたいです。
よろしくお願い申し上げます。

Re: 直線性の適合度検定
  投稿者:青木繁伸 2017/06/17(Sat) 15:10 No. 22397
もう40年も前にやっていたことで,細かいことは忘れてしまいましたね。

直線性の適合度検定も何通りかあったはずで,「この方法では、まず、物質の投与量に対する腫瘍発生率の変化の直線性について適合度検定」というのがどれか定かではないですが,たとえば,コクラン・アーミテージの検定は
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Armitage.html
で説明してあります。そのページにしたがってプログラムしてもよし,
http://aoki2.si.gunma-u.ac.jp/R/Cochran-Armitage.html
そこからリンクしている(R で計算してみる)ように R にすでにある prop.trend.test を使ってもよいでしょう。
http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Armitage-r.html

http://www.tamagaki.com/math/Statistics606.html の記事とは違います(一様性の分解という意味では関係あり)。

Re: 直線性の適合度検定
  投稿者:やまが 2017/06/17(Sat) 21:38 No. 22398
青木先生

ご回答いただきありがとうございます。
教えていただいたリンク先を確認し試してみます。
結果が得られ次第再度記入させていただきます。

Re: 直線性の適合度検定+最尤法
  投稿者:やまが 2017/06/30(Fri) 23:00 No. 22402
青木先生

直線性の適合度検定では大変お世話になりました。
さて、現在、上述のTD50を求めるため、
腫瘍発生確率: p=1-exp(-(b0+b1x))のパラメータ(b0, b1)を推定するため最尤法を行おうとしていますが、うまくいっておりません。元々、
http://www.yukms.com/biostat/haga/download/archive/likelihood/Likelihood.pdf
を参照し、Excelでは計算できていたのですが、
Rでも計算できるように取り組んでおります。

http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf
のロジスティック回帰を参考に、
p=1-exp(-(a + b・x))
Ln(1/(1-p))=a + b・x
とすればglmで同様に計算できるのではないかと考え、以下を試してみました。
> x <-c(0,1,3,9)
> r <-c(0,1,3,38)
> n <-c(49,50,50,50)
> p <- r/n
> plot(x,p)
>
> tumor<-cbind(n,n-r)
> glm.tumor<-glm(tumor~1+x,family=binomial)
>
> glm.tumor$coefficients["x"]
x
0.1539377
> glm.tumor$coefficients["(Intercept)"]
(Intercept)
-0.1495911
>
> x <- seq(min(x),max(x),by=0.01)
> eta.pred<- glm.tumor$coefficients["x"]*x+glm.tumor$coefficients["(Intercept)"]
> p.pred <- 1-exp(-eta.pred)
>
> lines(x, p.pred)
しかしながら、正しくない値が得られてしまいました。
そこで、対数尤度関数を最大化するパラメータを求めるために、
http://abrahamcow.hatenablog.com/entry/2016/04/09/085412
http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf
を参考にoptimを用いて以下を試みました。

> x <-c(0,1,3,9)
> r <-c(0,1,3,38)
> n <-c(49,50,50,50)
>
> LL<-function(par){
+ beta0<-par[1]
+ beta1<-par[2]
+ eta <- beta0+beta1*x
+ p<- 1-exp(-eta)
+ NLL<- sum(dbinom(r,n,p,log=TRUE))}
>
> opt1<-optim(c(1,1),LL,control=list(fnscale=-1))
There were 50 or more warnings (use warnings() to see the first 50)
> #opt1<-optim(c(1,1),LL,control=list(fnscale=-1),hessian=TRUE)
> opt1$par
[1] 5.405072e-09 9.328301e-02
> opt2<-optim(opt1$par,LL,control=list(fnscale=-1),method="BFGS",hessian=TRUE)
Error in optim(opt1$par, LL, control = list(fnscale = -1), method = "BFGS", :
non-finite finite-difference value [1]
In addition: Warning message:
In dbinom(r, n, p, log = TRUE) : NaNs produced
> my.std<-sqrt(-diag(solve(opt1$hessian)))
Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), :
'data' must be of a vector type, was 'NULL'
> up<-qnorm(0.995,opt1$par,my.std)
Error in qnorm(0.995, opt1$par, my.std) : object 'my.std' not found
> low<-qnorm(0.005,opt1$par,my.std)
Error in qnorm(0.005, opt1$par, my.std) : object 'my.std' not found
> round(matrix(c(low,up),2,2),3)
Error in matrix(c(low, up), 2, 2) : object 'low' not found

今度は正しい値が得られたのですが、
完全に収束していないのかhessianが得られませんでした。
(ちなみに以下のようにしてもだめでした)
 opt1<-optim(c(1,1),LL,control=list(fnscale=-1),hessian=TRUE)

大変申し訳ございません。
何らかの解決法がございましたらご教示いただけないでしょうか?
よろしくお願い申し上げます。

Re: 直線性の適合度検定
  投稿者:青木繁伸 2017/07/01(Sat) 10:10 No. 22403
いくつか間違いがありますね

tumor <- cbind(n, n-r)

ではなく,

tumor <- cbind(r,n)

です。しかし,

glm(tumor~1+x, family=binomial)

で得られるのは,「ロジスティック回帰」の結果なので,

eta.pred<- glm.tumor$coefficients["x"]*x+glm.tumor$coefficients["(Intercept)"]
p.pred <- 1-exp(-eta.pred) # ロジスティック曲線の式ではない

で予測値は計算できません。

ロジスティック回帰の結果を示すには,

y = 1/(1+exp(-predict(glm.tumor)))
points(x, y, col=2, pch=19)

x2 <- seq(min(x),max(x),by=0.01)
eta.pred<- glm.tumor$coefficients["x"]*x2+glm.tumor$coefficients["(Intercept)"]
p.pred <- 1/(1+exp(-eta.pred))
lines(x2, p.pred)

とすれば,「ロジスティック回帰」はうまくいっていることがわかるでしょう。

Re: 直線性の適合度検定
  投稿者:青木繁伸 2017/07/01(Sat) 15:50 No. 22404
optim を使う場合,対数尤度は
http://www.ma.utexas.edu/users/mks/RA/1hit.html
によれば,
LL2 = function(par) {
p = 1-exp(-par[1]-par[2]* x)
valid = (p != 0) & (p != 1)
p = p[valid]
n = n[valid]
r = r[valid]
sum(r*log(p) + (n-r)*log(1-p))
}
のようで,sum(dbinom(r, n, p, log = TRUE)) とはちがうのだけど,得られる答えは同じになる。

hessian = TRUE を指定すると
optim(c(1, 1), LL, control = list(fnscale = -1), hessian = TRUE)
解が得られない。

LL2 を使った
ans = optim(c(1, 1), LL2 , control= list(fnscale=-1), hessian=TRUE)
ans

では,エラーメッセージ付きで hessian は求まるが,1×1 である。

なお,1-exp(-par[1]-par[2]* x) は原点を通らないので,1-exp(-par[1]* x) とすると,
LL3 = function(q) {
p = 1-exp(-q*x)
valid = (p != 0) & (p != 1)
p = p[valid]
n = n[valid]
r = r[valid]
sum(r*log(p) + (n-r)*log(1-p))
}
として,optimize を使うのがよさそうではある。
ans2 = optimize(LL3, lower=0, upper=1, maximum=TRUE)
ans2

Re: 直線性の適合度検定
  投稿者:やまが 2017/07/02(Sun) 21:44 No. 22405
青木先生

丁寧なご回答をいただきありがとうございます。
大変助かりました。
確かに、Log(0)はよくないと気づかされました。
次は、このパラメータの信頼区間を求める必要があるのですが、
もう少し自分で考えてみようと思います。

結局再度質問させていただくことになるかもしれませんが、
今後ともよろしくお願いいたします。

Re: 直線性の適合度検定
  投稿者:やまが 2017/07/07(Fri) 22:23 No. 22410
いつもお世話になっております。
先日教えていただいた方法をもとに、
以下の資料のp.3の方法を参考に信頼区間を求めました。
http://www.agr.nagoya-u.ac.jp/~seitai/document/R2009/t090614glmIntro.pdf

x <-c(0,1,3,9)
r <-c(0,1,3,38)
n <-c(49,50,50,50)

LL3 <- function(q) {
p <- 1-exp(-q*x)
valid = (p != 0) & (p != 1)
p = p[valid]
n = n[valid]
r = r[valid]
sum(r*log(p) + (n-r)*log(1-p))
}
ans2 <- optimize(LL3, lower=0, upper=1, maximum=TRUE)
ans2

q<-ans2$maximum
n.sim <- 10000
beta <- data.frame(beta1 = rep(NA, n.sim))
eta <- q*x
p <- 1 - exp(-eta)

for(i in 1:n.sim){
r <- rbinom(length(p), n, p)

ans22 <- optimize(LL3, lower=0, upper=1, maximum=TRUE)
ans22$maximum
beta[i, ] <- ans22$maximum
}
beta.ci <- apply(beta, 2, quantile, prob = c(0.005, 0.995))
beta.ci

log(2)/ans2$maximum
log(2)/beta.ci[1,1]
log(2)/beta.ci[2,1]

一応うまくいっているようなのですが、
この方法だとどうしても計算のたびに異なる結果になります。
そこで、以下の資料のp.36,37を参考に求められないか検討しております。
http://www.yukms.com/biostat/haga/download/archive/likelihood/Likelihood.pdf
しかしながら、p.36に、
”#最尤推定値π^の対数尤度と帰無仮説πの対数尤度との差の2 倍は,近似的に,自由度1のカイ2 乗分布に従う”
という記述があるのですが、なぜそうなるのか基本的な考え方がわかっておりません。
もう少し詳しく解説されている書式やwebページ、または考え方のヒントがございましたらご教示いただけないでしょうか?
また、この方法以外に信頼区間を求める良い方法がございましたら併せてご教示いただけないでしょうか?
お手数をおかけして申し訳ございません。
よろしくお願い申し上げます。


入試問題の答え合わせ
  投稿者:大学院受験生 2017/07/03(Mon) 03:59 No. 22406
過去 20 年間の記録から、 A 地域の洪水の 1 年間における平均発生回数は、 4 (回 / 年)で あった。洪水の発生はポアソン過程と仮定する。以下の問に答えよ(ただし、 exp(-1)=0.368 , exp(-2)=0.135 , exp(-4)=0.018 として計算せよ)。
(問 1 )翌年に A 地域で 1 回も洪水が発生しない確率を求めよ。
(問 2 )翌年に A 地域で洪水の発生が 2 回以上となる確率を求めよ。
0.018
0.91
であってますか?

Re: 入試問題の答え合わせ
  投稿者:青木繁伸 2017/07/03(Mon) 10:59 No. 22407
あっているようです


多重ロジスティック回帰分析のモデル適合度について
  投稿者:明石 2017/06/28(Wed) 18:54 No. 22401
青木先生

ご回答いただきありがとうございます。
集計表を添付させていただきました。
ご指摘いただいた通り、正判別率については前例の予測が方群に偏っており、発現率がそのまま判別率となっております。

教えていただいたリンク先のURLを確認してみます。



多重ロジスティック回帰分析のモデル適合度について
  投稿者:明石 2017/06/27(Tue) 18:44 No. 22399
ロジスティック回帰分析のモデル適合度をホスマー・レメショウ検定で求めたのですが、P<0.01、正判別率90%という結果になりました。
モデルの適合度が良いとは言えなくても、正判別率が良いという状態をどう解釈すればいいのかわからずにいます。
分割表の25%以上のセルで期待値が5を下回っているので、カイ二乗分布に従わず、ホスマー・レメショウ検定が使えないということはあるのでしょうか。
サンプルサイズは70vs70。事象の発現率は10%です。
ご教授いただければ幸いです。

Re: 多重ロジスティック回帰分析のモデル適合度について
  投稿者:青木繁伸 2017/06/28(Wed) 07:25 No. 22400
せめて,ホスマー・レメショウ検定が適用された集計表(区分ごとの期待値と観察値),判別結果の集計表(2群×発現・非発現数)を示してください。

発現率が10%なら場合によっては90%という正判別率もあり得るでしょうし,期待値と発現数の相違もある程度大きいサンプルサイズのせいということもあるでしょう。

以下も参考になるでしょう。
内田治:ロジスティック回帰分析におけるモデルの適合度指標に関する考察と提案
http://www.iic.tuis.ac.jp/edoc/journal/ron/r8-1-2/index.html
要約すると,ホスマー・レメショウ検定は不適切ということ


R リストのファイル出力
  投稿者:明石 2017/06/11(Sun) 12:41 No. 22393
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

昨日の質問「R ファイル出力」では、大変にお世話になり、ありがとうございました。
おかげさまで、問題解決できました。

ファイル出力について、お聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

-----------

テキストマイニングの出力結果の一例です。
リスト構造の変数BOWには、文字列が代入されています。
(先頭の3要素のみ、お示しをします)

> BOW

[[1]]
[1] "わが国 景気 緩やか 拡大"

[[2]]
[1] "公共 投資 減少 傾向 輸出 増加 企業 収益 高水準 業況 感 良好 水準 推移"

[[3]]
[1] "輸出 海外 経済 拡大 背景 増加 国内 民間 需要 高水準 企業 収益 緩やか 増加"

このリストを読み取り、テキストファイルとして出力をしたいと考えています。

以下、期待する、テキストファイルの出力結果(先頭3行)です。

わが国 景気 緩やか 拡大
公共 投資 減少 傾向 輸出 増加 企業 収益 高水準 業況 感 良好 水準 推移
輸出 海外 経済 拡大 背景 増加 国内 民間 需要 高水準 企業 収益 緩やか 増加

しかしながら、私が検討しているプログラムでは以下の問題があります。

リストの要素ごとに、テキストファイルでは改行したいのですが(ここでは3行)、
1行になってしまい、苦慮しています。
改行コードの挿入方法、ファイルのappend等の知識が不足していることが原因だと思われます。

また、BOWの要素数でループを回していますが、apply関数を使うスマートな方法を知りたいと思います。

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

Re: R リストのファイル出力
  投稿者:青木繁伸 2017/06/11(Sun) 20:59 No. 22394
invisible(lapply(BOW, cat, "\n", file="test", append=TRUE))

【御礼】 Re: R リストのファイル出力
  投稿者:明石 2017/06/12(Mon) 09:02 No. 22395
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

ご教示をいただき、ありがとうございます。

listでapply関数を使う方法を、勉強させていただきました。

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


R ファイル出力
  投稿者:明石 2017/06/10(Sat) 20:46 No. 22390
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

先日の質問「R ベクトルの計算」では、大変にお世話になり、ありがとうございました。

また、Rについて、お聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

-----------

テキストマイニングのRパッケージを複数利用しています。
どのパッケージでも同じ問題にぶつかっており苦慮しています。

以下、パッケージ"textmineR"のマニュアルに掲載されている例題を例にとり、やりたいことをご説明します。

library(textmineR)

data("nih_sample")

# Convert a character vector to a document term matrix

dtm <- CreateDtm(doc_vec = nih_sample$ABSTRACT_TEXT,
doc_names = nih_sample$APPLICATION_ID,
ngram_window = c(1,1))

> class(dtm)
[1] "dgCMatrix"
attr(,"package")
[1] "Matrix"

> dim(dtm)
[1] 100 5210

document term matrix は、行方向は文書名、列方向は単語名からなる、行列です。

この計算結果では、100行*5210列のサイズをもつという情報が表示されています。

この行列をcsvファイル出力したいと考えています。

私が使う(知る)ファイル出力は write.csv ですので、それを使うと。

> write.csv(dtm, "dtm.csv")
as.data.frame.default(x[[i]], optional = TRUE) でエラー:
cannot coerce class "structure("dgCMatrix", package = "Matrix")" to a data.frame
>

行列なのでcsv出力できると思っていましたが、できません。

データ構造は、以下のようになっています。
> str(dtm)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
..@ i : int [1:14073] 99 99 99 99 99 99 99 99 99 99 ...
..@ p : int [1:5211] 0 1 2 3 4 5 6 7 8 9 ...
..@ Dim : int [1:2] 100 5210
..@ Dimnames:List of 2
.. ..$ : chr [1:100] "8693991" "8693362" "8607498" "8697008" ...
.. ..$ : chr [1:5210] "consumer" "attempt" "voucher" "allocated" ...
..@ x : num [1:14073] 1 1 1 1 1 1 1 1 1 1 ...
..@ factors : list()
>

この dtm の、100行*5210列の行列をcsvファイル出力できる方法について、ご教示いただけましたら助かります。

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

Re: R ファイル出力
  投稿者:青木繁伸 2017/06/11(Sun) 08:06 No. 22391
"R dgCMatrix write" で検索すると
https://stackoverflow.com/questions/4558277/write-a-sparse-matrix-to-a-csv-in-r
なんかが出てきて,
as.matrix で行列にできる(as.matrix(dtm) とする)というのが書かれているけど,
それは資源の無駄遣いなので writeMM を使うべき(?)とかも書かれています
https://stackoverflow.com/questions/17574099/in-r-how-do-you-write-a-sparse-matrix-to-a-file
には writeMM と readMM について書かれている。
その他のページも見てください。

Re: R ファイル出力
  投稿者:明石 2017/06/11(Sun) 10:46 No. 22392
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

いつもありがたいご教示をいただきまして、誠にありがとうございます。

助かりました。

さっそく勉強させていただきます。
ありがとうございました。


R ベクトルの計算
  投稿者:明石 2017/06/03(Sat) 15:18 No. 22383
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

先日の質問「R 文字列の抽出、突合せ」では、大変にお世話になり、ありがとうございました。

特に、先生からのご教示(No. 22373)で、
table関数の引数に文字列ベクトル与えると、要素の文字列の出現回数を得られることについては、大変に驚きました。
table関数の懐の大きさに感銘しました。Rの凄さに驚きました。

また、Rについて、お聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
何卒どうぞよろしくお願いいたします。

ーーー

表1は、単語のラベル、特徴ベクトルを示しています。
単語のラベルは w1,w2,w3,w4,w5 の5種類、特徴ベクトルは7次元とします。
(特徴ベクトルの数字自体には特に意味はありません)

表1 単語(ラベル、特徴ベクトル)
----------------------
ラベル 特徴ベクトル
----------------------
w1 1 2 3 4 3 2 1
w2 2 2 4 3 2 1 2
w3 1 3 3 2 1 2 3
w4 3 2 1 2 2 2 1
w5 1 4 2 2 3 4 1
----------------------

表2は、文に出現する単語の出現回数を示しています。
文の数は7(7行)で、各文に出現する単語のラベルは表1とリンクします。

表2 文(単語の出現回数)
-------------------------
w1 w2 w3 w4 w5
-------------------------
2 1 0 3 2
1 0 0 3 2
1 3 3 2 1
3 2 0 2 0
0 0 3 0 1
3 0 3 3 2
1 2 3 4 5
----------------------

表1と表2から、以下のルールで、文の特徴ベクトルを作成します。

【文の特徴ベクトルの作成ルール】
表2の各文について、
出現する各単語(表2)について、その特徴ベクトル(表1)に、出現回数(表2)を乗じたものを総計して、単語の出現総数で割る

【例示】
表2の1行目の文について、文の特徴ベクトルを計算すると
(2+w1ベクトル + 1* w2ベクトル + 3*w4ベクトル + 2*w5ベクトル)/(2+1+3+2)

表1、2は、行、列ともに大規模であることから、効率的な計算方法を検討しています。

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

Re: R ベクトルの計算
  投稿者: 2017/06/03(Sat) 21:10 No. 22384
tab1 <- matrix(c(1, 2, 3, 4, 3, 2, 1,
2, 2, 4, 3, 2, 1, 2,
1, 3, 3, 2, 1, 2, 3,
3, 2, 1, 2, 2, 2, 1,
1, 4, 2, 2, 3, 4, 1), 5, byrow = TRUE)

tab2 <- matrix(c(2, 1, 0, 3, 2,
1, 0, 0, 3, 2,
1, 3, 3, 2, 1,
3, 2, 0, 2, 0,
0, 0, 3, 0, 1,
3, 0, 3, 3, 2,
1, 2, 3, 4, 5),,5, byrow = TRUE)

t(sapply(seq(nrow(tab2)),
function(i) colSums(tab1 * tab2[i,])/sum(tab2[i,])))

[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1.875000 2.500000 2.125000 2.625000 2.500000 2.375000 1.125000
[2,] 2.000000 2.666667 1.666667 2.333333 2.500000 2.666667 1.000000
[3,] 1.700000 2.500000 2.800000 2.500000 1.900000 1.900000 1.900000
[4,] 1.857143 2.000000 2.714286 3.142857 2.428571 1.714286 1.285714
[5,] 1.000000 3.250000 2.750000 2.000000 1.500000 2.500000 2.500000
[6,] 1.545455 2.636364 2.272727 2.545455 2.181818 2.363636 1.545455
[7,] 1.666667 2.866667 2.266667 2.266667 2.200000 2.533333 1.533333

でお望みの結果が得られているでしょうか。

表2の1行目を例にすると
> tab1 * tab2[1,]
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 2 4 6 8 6 4 2
[2,] 2 2 4 3 2 1 2
[3,] 0 0 0 0 0 0 0
[4,] 9 6 3 6 6 6 3
[5,] 2 8 4 4 6 8 2

> colSums(tab1 * tab2[1,])
[1] 15 20 17 21 20 19 9

> sum(tab2[1,])
[1] 8

> colSums(tab1 * tab2[1,])/sum(tab2[1,])
[1] 1.875 2.500 2.125 2.625 2.500 2.375 1.125

をsapply関数でtab2の行ごとに計算し
t関数で行と列を入れ替えています。

確認をお願いします。

Re: R ベクトルの計算
  投稿者: 2017/06/04(Sun) 05:11 No. 22386
> (tab2 %*% tab1)/rowSums(tab2)
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 1.875000 2.500000 2.125000 2.625000 2.500000 2.375000 1.125000
[2,] 2.000000 2.666667 1.666667 2.333333 2.500000 2.666667 1.000000
[3,] 1.700000 2.500000 2.800000 2.500000 1.900000 1.900000 1.900000
[4,] 1.857143 2.000000 2.714286 3.142857 2.428571 1.714286 1.285714
[5,] 1.000000 3.250000 2.750000 2.000000 1.500000 2.500000 2.500000
[6,] 1.545455 2.636364 2.272727 2.545455 2.181818 2.363636 1.545455
[7,] 1.666667 2.866667 2.266667 2.266667 2.200000 2.533333 1.533333

の方が簡単ですね。

【御礼】 Re: R ベクトルの計算
  投稿者:明石 2017/06/04(Sun) 07:58 No. 22387
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

荒 様、
貴重なご教示をいただきまして、ありがとうございます。

確認させていただきました。
ここで使われている考え方は、大変に勉強になりました。

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


順位で並べて分けた2群の差の検定
  投稿者:朱鷺 2017/06/02(Fri) 12:28 No. 22381
お世話になります。

10名の研究協力者のスコアをもとに、スコアの低い5名と、スコアの高い5名の2群に区分し、ある実験をおこないその結果を比較しようとしていたところ、
まず「低スコア群」と「高スコア群」の2群間に統計的に有意な差があることを示すよう指摘を受けました。

この場合、標本数が少ないのでt検定ではなく、マン・ホイットニーのU検定になるのではないかと思ったのですが、
そもそもスコアの低い5名と高い5名に分けた2群ですので、順位尺度を用いるU検定では2群の差が有意になるのは当たり前のような気がします。

そこで質問ですが、(1) まず上記の解釈は間違っていますでしょうか。 (2) この場合、U検定を使うべきでしょうか。それともt検定などの他の方法を用いるべきなのでしょうか。

初歩的な質問ですが、ご教示のほど、どうぞよろしくお願いいたします。

Re: 順位で並べて分けた2群の差の検定
  投稿者:青木繁伸 2017/06/03(Sat) 14:27 No. 22382
> 「低スコア群」と「高スコア群」の2群間に統計的に有意な差があることを示すよう指摘

がおかしいのは言うまでもないですが,

> スコアの低い5名と、スコアの高い5名の2群に区分し、ある実験をおこないその結果を比較

もよくないでしょう。

スコアと,実験結果の相関をみるほうが自然で妥当ではないですか?

Re: 順位で並べて分けた2群の差の検定
  投稿者:朱鷺 2017/06/03(Sat) 23:09 No. 22385
青木先生、
ご教示くださいましてありがとうございました。

>> 「低スコア群」と「高スコア群」の2群間に統計的に有意な差があることを示すよう指摘
> がおかしいのは言うまでもないですが,

そのような指摘を受け、ずっともやもやしていたのですが、すっきりしました。
ご教示いただいた方向で分析することにします。感謝申し上げます。


グラフ
  投稿者:課長 2017/05/31(Wed) 22:58 No. 22378
ポアソン分布や二項分布が正規分布で近似できるところをグラフで視覚的に見てみたくて、Rでグラフを書こうとしたのですが、うまくゆきません。
正規分布のグラフは
curve(dnorm(x,mean=1,sd=1))
などでうまくゆくのですが、ポアソン分布などはうまくゆきません
curve(dpois(x,lambda=2))などとすると謎の図形が出てきてエラーも50くらい出てきます。
これまで統計をうわべだけ学んでいて基礎から学んだらわからないことだらけであることに気づきました。程度の低い質問ですが、どうかご教示いただければと思います。

Re: グラフ
  投稿者:青木繁伸 2017/05/31(Wed) 23:57 No. 22379
参考までに
# 二項分布
n = 20 # サンプルサイズ
p = 0.3 # 母比率
y = dbinom(0:n, n, prob=p)
x = barplot(y, space=0, ylim=c(0, max(y)*1.1))
text(x, -0.01, 0:n, xpd=TRUE)
Mean = n*p
SD = sqrt(n*p*(1-p))
y = dnorm(0:n, mean=Mean, sd=SD)
lines(0:n+0.5, y, col="red")
points(0:n+0.5, y, pch=19, col="red")

# ポアソン分布
n = 20 # サンプルサイズ
lambda = 10 # ポアソン定数
y = dpois(0:n, lambda)
x = barplot(y, space=0, ylim=c(0, max(y)*1.1))
text(x, -0.01, 0:n, xpd=TRUE)
Mean = lambda
SD = sqrt(lambda)
y = dnorm(0:n, mean=Mean, sd=SD)
lines(0:n+0.5, y, col="red")
points(0:n+0.5, y, pch=19, col="red")

Re: グラフ
  投稿者:課長 2017/06/01(Thu) 03:29 No. 22380
青木先生、いつもこんな未熟者のためにご丁寧な回答ありがとうございます。いただいたプログラムをまさに解剖するかのように分析して学習しようと思います。統計の基礎をやりなおすうえで、Rは強力な武器になると思ってますので、これからももしかしたら初歩的な質問をさせてもらうこともあるかもしれませんが、なるべく十分調べてからの質問にしようと思います。ありがとうございます。


共分散
  投稿者:課長 2017/05/30(Tue) 04:13 No. 22375
互いに独立な確率変数X1,X2,X3があり
COV(X1,(X1+X2+X3)/3)=1/3*COV(X1、X1)

とテキストにありますが、これを証明できません

なんでこうなるのかわからず、こんな未熟な質問をしてすみません

Re: 共分散
  投稿者:青木繁伸 2017/05/31(Wed) 08:57 No. 22376
LaTeX で書くといいとは思うのですが,プログラミング言語の数式も厳密なものにちがいないので,R で書いてみます。sum をΣ,mean を $\bar{x}$ に置き換えて普通の数式にしてみてください。
Cov の定義は,a, b の二変数について,
Cov = function(a, b) {
n = length(a)
(1/n)*sum((a-mean(a))*(b-mean(b)))
}
1/n を 1/(n-1) にすれば R の Cov と同じ
Cov = function(a, b){
n = length(a)
cov(a, b)*(n-1)/n
}

w = (x+y+z)/3 とすれば
Cov(x, w) = (1/n)*sum((x-mean(x))*(w-mean(w)))
= (1/n)*sum(x*w) - (1/n)*sum(x*mean(w)) + (1/n)*sum(mean(x)*w) + (1/n)*sum(mean(x)*mean(w))
= (1/n)*sum(x*w) - mean(x)*mean(w) + mean(x)*mean(w) + mean(x)*mean(w)
= (1/n)*sum(x*w) + mean(x)*mean(w)
= (1/n)*sum(x*((x+y+z)/3)) + mean(x)*mean((x+y+z)/3)
= (1/n)*sum(x*(x/3)) + (1/n)*sum(x*(y/3)) + (1/n)*sum(x*(z/3)) + mean(x)*mean(x/3) + mean(x)*mean(y/3) + mean(x)*mean(z/3)
= (1/n)*sum((x-mean(x))*((x/3)-mean(x/3))) + (1/n)*sum((x-mean(x))*((y/3)-mean(y/3))) + (1/n)*sum((x-mean(x))*((z/3)-mean(z/3)))
= Cov(x, x/3) + Cov(x, y/3) + Cov(x, z/3)
= Cov(x, x/3) + Cov(x, y)/3 + Cov(x, z)/3
= Cov(x, x/3) + 0/3 + 0/3   ∵ x, y, z は独立

Re: 共分散
  投稿者:課長 2017/05/31(Wed) 18:43 No. 22377
本当にご丁寧な解説本当にありがとうございます。プログラムから指揮を書き出して理解できました。すっきりしました。本当にありがとうございます。もう少し統計の基礎をやり直してゆきます。


R 文字列の抽出、突合せ
  投稿者:明石 2017/05/28(Sun) 07:22 No. 22370
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

昨日の質問では、大変にお世話になり、ありがとうございました。

また、Rについて、お聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
どうぞよろしくお願いいたします。

ーーー

以下に例示します。
使うデータは、(1)(2)の2つです。

(1) 文字列

  " 政治 経済 政治 社会 経済 国民 "

複数の単語を含む、文字列です。
単語の区切りは、1個以上の半角空白文字です。
文字列の先頭、最後にも、1個以上の半角空白文字があります。

(2) 単語の表(表1とよびます)

ban word
1 政治
2 経済
3 社会
4 外交
5 教育
6 福祉
7 防衛

求めたい出力は、以下にお示しをします表2です。

(1)の文字列から単語を抽出して、
その単語を、(2)表1と突合せをして、表2を作成します。

表2
word freq ban
政治 2 1
経済 2 2
社会 1 3
国民 1 NA

項目の意味は、以下の通りです。
word  抽出した単語
freq  文字列の中の出現度数
ban   表1における位置(行番号)
     もし、表1に単語がない場合には、NA

(1)の文字列から単語を抽出する際に、文字区切りが半角空白1文字ならば、
text <- " 政治 経済 政治 社会 経済 国民 "
s <- strsplit(text, " ")
のように、簡単に取り出しができるのですが、
1個以上の半角空白文字があることから、苦慮しています。

お示しをしました例示は、小さいデータですが、
やりたいデータでは、表1のサイズ、文字列ともに大規模なので、
効率のよい処理も期待しています。

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

Re: R 文字列の抽出、突合せ
  投稿者:青木繁伸 2017/05/28(Sun) 08:40 No. 22371
先頭,末尾の空白は sub 関数で除去できます
strsplit の区切り文字にも,正規表現が使えます
" +" は,「連続する一個以上の半角空白」を意味します
> s = "      政治      経済    政治 社会   経済 国民          "
> s = sub("^ +", "", s) # 先頭の空白を除去
> # s = sub(" +$", "", s) # 末尾の空白は除去しておかなくてもよい
> unlist(strsplit(s, " +"))
[1] "政治" "経済" "政治" "社会" "経済" "国民"

Re: R 文字列の抽出、突合せ
  投稿者:明石 2017/05/28(Sun) 11:06 No. 22372
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

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

strsplit関数で、区切り文字の指定に正規表現が利用できることを教えていただきましたので、
今後の活用範囲が広がります。
ありがとうございました。

ご教示いただきましたコードを活用させていただき、
目的の表2の作成に向けて、自分でコードを作成しました。

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

s <- " 政治 経済 政治 社会 経済 国民 "

s <- sub("^ +", "", s) # 先頭の空白を除去
# s = sub(" +$", "", s) # 末尾の空白は除去しておかなくてもよい

s <- unlist(strsplit(s, " +"))
# [1] "政治" "経済" "政治" "社会" "経済" "国民"

d <- unique(s)
# [1] "政治" "経済" "社会" "国民"

tb1 <- data.frame(
ban = 1:7,
word = c("政治", "経済", "社会", "外交", "教育", "福祉", "防衛"), stringsAsFactors=F )

n <- length(d)

word <- character(n)
freq <- numeric(n)
ban <- numeric(n)

for (i in 1:n) {
word[i] <- d[i]
freq[i] <- length(which(s == d[i]))
ban[i] <- ifelse( any(tb1$word == d[i]), which(tb1$word ==d[i]) ,"NA")
}

tb2 <- data.frame(word, freq, ban)
tb2

word freq ban
1 政治 2 1
2 経済 2 2
3 社会 1 3
4 国民 1 NA

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

Re: R 文字列の抽出、突合せ
  投稿者:青木繁伸 2017/05/28(Sun) 14:58 No. 22373
> s = " 政治 経済 政治 社会 経済 国民 "
> word = c("政治", "経済", "社会", "外交", "教育", "福祉", "防衛")

> s = sub("^ +", "", s)
> s = unlist(strsplit(s, " +"))

> freq = table(s) # 集計は table で
> word2 = names(freq)

> ban = unname(sapply(word2, function(w) ifelse(w %in% word, grep(w, word), NA)))
> tb2 = data.frame(word = word2, freq = as.integer(freq), ban=ban)
> tb2
word freq ban
1 経済 2 2
2 国民 1 NA
3 社会 1 3
4 政治 2 1
> # 並べ替えが必要なら以下の2行
> tb2 = tb2[order(ban),]
> rownames(tb2) = seq_along(ban)
> # 結果
> tb2
word freq ban
1 政治 2 1
2 経済 2 2
3 社会 1 3
4 国民 1 NA

Re: R 文字列の抽出、突合せ
  投稿者:明石 2017/05/28(Sun) 15:17 No. 22374
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

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

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


ベクトルが作れない
  投稿者:課長 2017/05/26(Fri) 05:02 No. 22358
Rでベクトルを作ろうと思ったら、次のようなメッセージが出てしまいます。

> x<-c(1,1,1,1,1,1)
c(1, 1, 1, 1, 1, 1) でエラー:
6 個の引数が 'exp' に渡されましたが、1 が必要とされています

どうして急にベクトルが作れなくなったのか、ご教示いただければ幸いです。

Re: ベクトルが作れない
  投稿者:青木繁伸 2017/05/26(Fri) 11:35 No. 22359
どこかで,c という関数を再定義してしまったのではないでしょうか。

コンソールに c と入力して
> c
function (...) .Primitive("c")
以外のものが表示されるようならば,確認の上で,
> rm(c)
と入力して,その再定義されたものを消去すればよいでしょう。

Re: ベクトルが作れない
  投稿者:課長 2017/05/26(Fri) 20:53 No. 22363
青木先生のおっしゃる通りでした。なぜかCに関数が定義されていました。そしてご指導通りにしたらなおりました。ありがとうございます

Re: ベクトルが作れない
  投稿者:課長 2017/05/26(Fri) 22:49 No. 22364
青木先生の言われるとおりにしたらなおったのですが、Rを再起動するとまたcに謎の関数が定義されてしまいます。これはどうしたら直るかご教示いただけませんか。

> c
function (x) .Primitive("exp")

が毎回出ます。

Re: ベクトルが作れない
  投稿者:青木繁伸 2017/05/27(Sat) 07:00 No. 22365
R を起動してすぐに rm(c) をしたあと,quit(save="yes") でワークスペースを保存して R を終了してください。
そのあと R を再起動すると c の定義は永久になくなると思います。
もしそれでもなくならないなら,実行されるファイルのどこかに c <- exp というような c を再定義する箇所がないか捜してください。
なお,蛇足ですが,c を exp のような関数ではなく定数を代入しても,c 関数は置き換えられません。
> c <- 1
> c(9, 3, 1)
[1] 9 3 1

Re: ベクトルが作れない
  投稿者:課長 2017/05/27(Sat) 21:14 No. 22369
ありがとうございます、青木先生!
おっしゃるとおりにいたしましたところ、直りました。
また、勉強にもなりました。
心より感謝いたします。
今後ともよろしくお願いいたします・


R レイアウトが異なるテキストファイルの読み込み
  投稿者:明石 2017/05/26(Fri) 14:33 No. 22360
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

Rについて、お聞きしたいことがでてきました。
ご教示をいただければ、大変に助かります。
どうぞよろしくお願いいたします。

ーーー

添付の画像ファイルをご覧ください。

読み込むテキストファイルは、
1行目
2行目〜
と列の長さが異なることから、ファイルの読み込みで苦慮しています。

データの区切りは、半角空白文字です。

1行目の2つの数字(300000 200)は、2行目以降の表のサイズを表しています。
題材によりサイズは異なります。

2行目以降に、ラベル(文字型)、数値(200列)の組が、300000行 並んでいます。

出力したいものは、
1行目に記載されている表のサイズ(300000 200)をもつ、2行目以降のデータのmatrixです。
・ラベル(文字型)を、rownamesに格納
・colnames は、paste0("V",1:200)
・数値を数値型で、200列

2行目以降の表のサイズが大きいので、効率的な処理もできればと思います。

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


Re: R レイアウトが異なるテキストファイルの読み込み
  投稿者:青木繁伸 2017/05/26(Fri) 18:27 No. 22361
1 行目を読み飛ばせば,あとは,1行目に書いてあるサイズ情報を利用せずとも,read.table がちゃんと読んでくれます。
read.table(ファイル名, skip=1)
colnames もデフォルトで付けてくれます。

Re: R レイアウトが異なるテキストファイルの読み込み
  投稿者:明石 2017/05/26(Fri) 18:36 No. 22362
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

ご教示をいただきまして、大変に助かりました。
ありがとうございました。

Re: R レイアウトが異なるテキストファイルの読み込み
  投稿者:明石 2017/05/27(Sat) 15:23 No. 22366
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

昨日はご教示をいただき、大変に助かりました。
ありがとうございました。

ーーー

昨日の画像ファイルのデータ例について、追加のご質問がございます。 
どうぞよろしくお願いいたします。

1行目の 300000 200 は、2行目以降の表データのサイズ(行 列)を示しています。

「文字」+「数値(200列)」の組が、300000行の繰り返しがあることを示しています。

(表のサイズは、題材により異なります。 固定ではありません)

2行目以降の表データのサイズ(行 列)を、1行目で明示している理由は、
ケースによっては、数値(200列)の後半に欠損がある場合を想定してからです。

件名の「レイアウトが異なるテキストファイルの読み込み」の真意は、ここにあります。

先生にご教示いただきたいことは、以下です。
お手数をおかけいたします。

・2行目以降に読み込む表データを、matrix型で格納する

・読み込んだ「文字」は、rownamesに格納する

・読み込んだ「数値(200列)」は、1列〜200列に数値データとして格納する。
 
・ケースによっては、数値(200列)の後半に欠損がある場合への対応として、
 データに欠損がある行については、列数が合うように(この場合は 200)
 その箇所をNA で補う。

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

ご教示いただければ、大変に助かります。

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

Re: R レイアウトが異なるテキストファイルの読み込み
  投稿者:青木繁伸 2017/05/27(Sat) 16:05 No. 22367
何度かお願いしておりましたが,添付図ではなく,仕様の全てを含む簡単なテストデータと,期待する結果をそえて質問してください。

説明文は曖昧さを含みますし,読み取りが不十分なこともありますから。
その点,データと結果があれば,仕様を満たすプログラムは一意に定まりますから。

「ケースによっては、数値(200列)の後半に欠損がある」というのは,一行ごとに数値の個数が違うことがあるということですか?
以下のようなものでよいですか?
file: "test.dat"
=====次行から=====
4 6
a 1 2 3 4 5
b 11 12 13
c 21 22 23 24 25 26
d 31 32 33 34
=====前行まで=====

func = function(fileName) {
s = readLines(fileName, 1)
s = as.numeric(unlist(strsplit(s, " ")))
ncol = s[2]
d = read.table(fileName, skip=1, row.names=1, header=FALSE, col.names=c("name", paste0("V", 1:ncol)), fill= TRUE)
as.matrix(d)
}
func("test.dat")

結果
> a = func("test.dat")
> a
V1 V2 V3 V4 V5 V6
a 1 2 3 4 5 NA
b 11 12 13 NA NA NA
c 21 22 23 24 25 26
d 31 32 33 34 NA NA
> class(a)
[1] "matrix"

【御礼】 Re: R レイアウトが異なるテキストファイルの読み込み
  投稿者:明石 2017/05/27(Sat) 17:00 No. 22368
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

今回もご丁寧にご教示いただき、誠にありがとうございます。
御礼を申し上げます。

先生からご指摘を受けましたように、
今後、質問をさせていただく際には、お書きになられたことに沿うようにいたします。

不愉快な思いをさせてしまい、大変に申し訳ございませんでした。
深くお詫びをいたします。

ーーー

今回の質問の内容については、先生がお書きになられた通りです。
私が期待していました通りです。

青木先生からのご教示は、大変に良いお手本となります。
さっそく、勉強をさせていただきます。

今回もありがたいご教示に感謝いたします。
ありがとうございました。

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

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