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

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


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

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


コクラン・アーミテージ検定でのセルごとの数
  投稿者:木村 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 を教えてください。


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増分あたりの効果を、回帰係数として求めているので、負の値でも正の値でも使えるというわけです。


信頼性の検定におけるサンプルサイズ
  投稿者: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
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

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

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

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

ーーー

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

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

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


連続変数の自由度
  投稿者:統計初心者 2017/05/17(Wed) 17:56 No. 22352
お世話になります。

一般線形モデルで、従属変数に連続変数のパラメーターを入れて検定したとき(サンプル数は48)、連続変数のパラメーターの自由度が1になるのはなぜでしょうか?
サンプル数は48なので、48-1=47になると思うのですが。
統計ソフトはJMP13を使用しています。

お忙しいところ恐れ入りますが、ご回答よろしくお願いいたします。

Re: 連続変数の自由度
  投稿者:統計初心者 2017/05/18(Thu) 15:56 No. 22353
画像を追加します。
図の赤丸で囲った部分です。


Re: 連続変数の自由度
  投稿者:青木繁伸 2017/05/19(Fri) 08:07 No. 22354
この類の質問は,回答に困りますね。
変数の不偏分散の自由度(n-1)と勘違いされているのでは?

Re: 連続変数の自由度
  投稿者:統計初心者 2017/05/19(Fri) 11:46 No. 22355
青木先生

ご返答ありがとうございます。
この場合だと、線形モデルの式が y = xβ + e で、xが決まればyも決まるから、パラメーターの自由度は1という理解で大丈夫ですか?

Re: 連続変数の自由度
  投稿者:青木繁伸 2017/05/19(Fri) 13:19 No. 22356
大元から掘り起こして理解するということもたいせつとは思いますが,「そういう定義である」と割り切って先に進む方が建設的では?

http://aoki2.si.gunma-u.ac.jp/lecture/Regression/mreg/mreg4.html
にある分散分析表の構成法と解釈法を理解することのほうがたいせつでは?
なお,そこの表ではp個の独立変数を持つ重回帰分析なので,一変量の単回帰分析なら p=1 です。重回帰分析の場合も,個々の独立変数の自由度は1です。

Re: 連続変数の自由度
  投稿者:統計初心者 2017/05/19(Fri) 18:37 No. 22357
青木先生

大変参考になりました。ありがとうございます。


R ポアソン回帰分析 glm関数の引数offsetの使い方
  投稿者:明石 2017/05/08(Mon) 22:03 No. 22349
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

先日は、ポアソン分布のパラメータ計算方法について、分かりやすくご教示をいただき、誠にありがとうございました。

この投稿を契機に、Rのglm関数を用いて、ポアソン回帰分析の勉強を進めています。

glm関数のフォーミューラ(目的変数~説明変数)の記述方法は理解できました。

しかし、引数offsetの使い方が分かりません。

説明変数のうち、
どの変数をフォーミューラの右辺に書き、
どの変数を引数offsetに書くのか、
その使い分けが分かりません。

分かりやすい例示をしていただけると、助かります。
お手数をおかけいたします。
どうぞ、よろしくお願いいたします。

Re: R ポアソン回帰分析 glm関数の引数offsetの使い方
  投稿者:青木繁伸 2017/05/09(Tue) 06:20 No. 22350
example(glm) の中程以降に,family = gaussian ではありますが,offset を使う例が出てきますが...
glm> anorex.1 <- glm(Postwt ~ Prewt + Treat + offset(Prewt),
glm+ family = gaussian, data = anorexia)
"R glm poisson offset" の 4 語で検索すると,いくつか出てきますが,
講義のーと: データ解析のための統計モデリング
eprints.lib.hokudai.ac.jp/dspace/bitstream/2115/49477/5/kubostat2008d.pdf
久保拓弥 著 - ‎2008
などにも例があるようです。

Re: R ポアソン回帰分析 glm関数の引数offsetの使い方
  投稿者:明石 2017/05/09(Tue) 09:05 No. 22351
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

ご紹介してくださいました資料を参照します。
お手数をお掛けして、申し訳ございませんでした。

いつも有難いご教示をいただき、大変に感謝しています。


共分散構造分析の重決定係数につきまして
  投稿者:YUYU 2017/05/04(Thu) 20:04 No. 22341
青木先生

こんばんは。
はじめて質問させていただきます。YUYUと申します。

子育て中の看護師のバーンアウト因果モデルを作成しています。
ある集団に、この因果モデルを当てはめたところ、標準化推定値の重決定係数が、0.97となりました。適合度は、GFI = .925, AGFI = .890, RMSEA = .049でした。

ほかの集団での重決定係数は、おおむね、0.60から0.70でした。
0.97は、高すぎると思うのですが、いかがでしょうか。
ぜひ、先生のご意見をお伺いいたしたく、メールをさせていただきました。

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

Re: 共分散構造分析の重決定係数につきまして
  投稿者:青木繁伸 2017/05/05(Fri) 08:39 No. 22342
0.97 という一つの数値だけ示されても,何ともいえないですね

Re: 共分散構造分析の重決定係数につきまして
  投稿者:YUYU 2017/05/06(Sat) 10:04 No. 22343
青木先生

おはようございます。
お返事をいただきましてどうもありがとうございます。

当該の因果モデルを添付させていただこうと思ったのですが、どうやってもアップロードできませんでした。

一点、質問させてください。相関係数が高すぎると同質のものであることが疑われるといわれるように、重決定係数にも、一般的におかしいと考えられる値があるのでしょうか?

Re: 共分散構造分析の重決定係数につきまして
  投稿者:青木繁伸 2017/05/06(Sat) 11:54 No. 22345
> どうやってもアップロードできませんでした。

本文中にテキストとして挿入するだけでよいのでは?

> 重決定係数にも、一般的におかしいと考えられる値があるのでしょうか?

何をもっておかしいとするか,に尽きるでしょう。

おかしな結果が出た理由は,分析結果ではなくデータそのものにある。
以下のような数値だけの結果を見せられても,おかしいかどうかわからないでしょう。
とても優れた予測式であるということかも知れない。
> x1 = c(3, 2, 1, 4, 2, 5, 6)
> x2 = c(4, 1, 2, 5, 4, 7, 6)
> y = c(17, 7, 9, 23, 15, 27, 35)
> a = lm(y ~ x1 + x2)
> summary(a)

Call:
lm(formula = y ~ x1 + x2)

Residuals:
1 2 3 4 5 6 7
-0.69528 -2.30901 1.96567 -0.02575 1.10730 -2.88412 2.84120

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1760 2.3277 0.076 0.9434
x1 3.8026 1.1733 3.241 0.0316
x2 1.5279 0.9979 1.531 0.2005

Residual standard error: 2.612 on 4 degrees of freedom
Multiple R-squared: 0.9545, Adjusted R-squared: 0.9318
F-statistic: 41.96 on 2 and 4 DF, p-value: 0.00207
以下に添付するのは,とてもよい(?)予測式が得られたというデータ例。標高の係数が絶妙。


Re: 共分散構造分析の重決定係数につきまして
  投稿者:YUYU 2017/05/06(Sat) 17:28 No. 22346
青木先生

ご回答いただきまして、どうもありがとうございます。
大変勉強になります。

因果モデルのテキストへの貼り付けの件、試してみましたが張り付きませんでした。
図として添付してみます。

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


Re: 共分散構造分析の重決定係数につきまして
  投稿者:青木繁伸 2017/05/06(Sat) 19:14 No. 22347
テキストの結果でなければテキストに挿入することはそもそも不可能ですが。

変数,潜在因子が多すぎるのではないかと思います。重回帰分析でも,独立変数をどんどん増やしていけば決定係数は大きくなりますので(自由度調整済み決定係数ならそのようなことはない)
解釈可能ならそれで良いかも知れませんが,より簡潔なモデルも検討してみたらいかがですか。

Re: 共分散構造分析の重決定係数につきまして
  投稿者:YUYU 2017/05/07(Sun) 09:03 No. 22348
青木先生

おはようございます。
先生のアドバイスを参考に、他モデルも検討してみようと思います。

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


ポアソン分布の「プロイセン陸軍で馬に蹴られて死亡した兵士数」の事例
  投稿者:明石 2017/05/03(Wed) 14:18 No. 22338
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

今回は、Rではなく、ポアソン分布について調べ物をしております。

Google先生に色々とお聞きしましたが分かりませんでしたので、投稿させていただきます。
どうぞよろしくお願いします。

−−−−−−−−−−
ポアソン分布の歴史的に有名な事例として、
「ロシアのプロイセン陸軍で馬に蹴られて死亡した兵士数」の事例があります。

Wikipedia他サイトでは、以下のように記載されています。
「プロイセン陸軍の14の騎兵連隊の中で、
 1875年から1894年にかけての20年間で馬に蹴られて死亡する兵士の数について
 調査しており、1年間当たりに換算した当該事案の発生件数の分布が
 パラメータ0.61のポアソン分布によく従うことを示している。」

これ以上に詳しい情報を得ることはできませんでした。

パラメータ0.61をどのように見つけたのか、計算したのか、適合度を確認したのか、
その方法に強い興味があります。

もし、これについて参考になる情報源をご存知でしたら、教えていただけましたらとても助かります。

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

Re: ポアソン分布の「プロイセン陸軍で馬に蹴られて死亡した兵士数」の事例
  投稿者:青木繁伸 2017/05/03(Wed) 17:40 No. 22339
ポアソン分布のパラメータは「期待値(平均値)」です。
> 死亡した兵士の数 = 0:4
> 軍団数 = c(109, 65, 22, 3, 1)
> mean(rep(死亡した兵士の数, 軍団数))
[1] 0.61

> rbind(死亡した兵士の数, 軍団数)
[,1] [,2] [,3] [,4] [,5]
死亡した兵士の数 0 1 2 3 4
軍団数 109 65 22 3 1
> sum(軍団数)
[1] 200
> (0*109+1*65+2*22+3*3+4*1) / 200
[1] 0.61

Re: ポアソン分布の「プロイセン陸軍で馬に蹴られて死亡した兵士数」の事例
  投稿者:明石 2017/05/03(Wed) 18:43 No. 22340
青木先生、
いつもお世話になり、ありがとうございます、明石と申します。

お休みのところ、ご丁寧にご教示をいただきまして、誠にありがとうございました。

一般化線形モデルのポアソン回帰モデルで、パラメータを推定するのかと思っておりました。

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


尺度開発 テストー再テストの2回の集団の代表性を示す意味
  投稿者:kashio 2017/04/28(Fri) 19:10 No. 22332
いつもご助言ありがとうございます

尺度開発の論文を投稿しました。査読者から以下4の指摘がありました。

4.研究結果
‖仂櫃瞭胆について、表1には1回目調査の結果のみですが、2回目調査に協力した集団(152名)は1回目調査の集団(774名)と特性が変わらないという前提でしょうか。1回目と2回目の調査対象に違いがないことを示す必要があります(サンプル代表性の担保)。つまり、表1には1回目と2回目の調査対象に関する記述集計を併記し、群間差があるか検定結果を示してください。検定方法については、p7のD.分析方法へ加筆してください。

これがどうもわかりません。

ほとんどの病院が2回(テスト―リテスト)はいやというので
1回目の調査の中の3病院が共同研究者の病院でここに頼みました。
300人ほどに頼みましたが、有効回答は半数以下でした。これもあり
テストーリテストは0・6で確かに低く課題が残りました。今後の課題だと書くつもりです。

2回目のテストは、テスト再テスト法にのみ使用します。
人数が少ないので、年齢、子供の有無、教育背景などにも差が出ますが、
示す意味が分からず、差があるから、なんだというのかもわかりません。
1回目の因子分析等の集団と再テストの集団と母集団が違うとだめ
違うのでしょうか?
ほかの尺度開発の論文をいくつか読みましたが、2回目の特性を示しているものは一つもありません。2回目は、1回目の一部です。そしてその同じ集団でテストーリテストの解析です。
代表性を示す意味も分かりません。

Re: 尺度開発 テストー再テストの2回の集団の代表性を示す意味
  投稿者:青木繁伸 2017/04/29(Sat) 22:39 No. 22334
極端な場合を考えてみましょう。

尺度構成は800人(男女400人ずつ)のデータで行った。
2回目の調査は(極端にもほどがあるが)女150人に行った。
再現性の検討は2回とも調査できた女150人分のデータで行った。...
これでは,再現性があるなしをいえませんね。

そのあたりの担保が欲しいのでは?つまり,あなたの2回の調査で,対象者の背景因子に大きな違いがないと言って欲しいのかな。

しかし,再現性の検討は,尺度構成に使った集団とは違う集団に対して(妥当性もチェックできるし),もう少しサンプルサイズが大きい調査を行う必要があると思います。サンプルサイズが小さいと,「再現性がある」のか「再現性がないのだけど,それを検出できないだけ」なのかわかりませんからね。再現性を言うためにはなるべく少ないサンプルサイズで調査するのが...などと。

Re: 尺度開発 テストー再テストの2回の集団の代表性を示す意味
  投稿者:kashio 2017/04/29(Sat) 23:36 No. 22335
ありがとうございました。男性が少なかったので女性だけを対象としてやり直しました。
青木先生の指摘通り、テスト再テストの人数が少ないうえ、共同研究者の大学の関連病院3つにしたため、特性も差が出ました。尺度の平均には差が出なかったのですが、年齢、実務職種、最終学歴に差が出ました。もうどうにも考察しようがないでしょうか?
探索的因子分析、確証的因子分析の結果は良好でした。クロンバックαは、0.69から0.84でした。テストリテストは、の相関係数は、0.6なので妥当性、信頼性はおおむね良好であった とは言えないでしょうか?

Re: 尺度開発 テストー再テストの2回の集団の代表性を示す意味
  投稿者:青木繁伸 2017/04/30(Sun) 17:46 No. 22336
まさかの「男性が少なかったので女性だけを対象としてやり直しました。」ですか。

調査デザインとして問題が多いでしょうね。

Re: 尺度開発 テストー再テストの2回の集団の代表性を示す意味
  投稿者:kashio 2017/05/01(Mon) 05:17 No. 22337
お返事ありがとうございます。対象が看護師なので最初から2000人以上対象にしないと男性看護師までは、無理かと危惧されたのですが、2回の調査の協力病院が得られず、尺度の平均にも性差が出たので、査読者に指摘されやり直しました。テストー再テストも0.6だったのでつらいですが、信頼性に課題が残ったが妥当性は確保されたとして査読者に返します。いつもありがとうございます。


対応のあるt検定?
  投稿者:いちご 2017/04/23(Sun) 15:07 No. 22329
お世話になっております。
初歩的な質問だと思うのですが、教えていただければ幸いです。

4件法の尺度(67項目)に対して、自分のことを回答した群(A群)とその人のことをよく知っている人がその人のことを回答した群(B群)の対応データが300ペア分あります。
このとき、A群の回答とB群の回答が等しいかそうでないかを検討するのに、どんな分析を行えばよいのでしょうか。
対応のあるt検定かなと思ったのですが、違うような気もします。

どなたか、ご教示いただけないでしょうか。
よろしくお願いいたします。

Re: 対応のあるt検定?
  投稿者:課長 2017/04/24(Mon) 18:52 No. 22331
4件法ということであれば、各項目についてウィルコクソンの符号付き順位和検定のほうがふさわしいかと思いますが、67項目あるということは67回も検定を行うということになりそうですね。それはそれで問題かと思いますねぇ。多変量分散分析か、因子分析して各因子の合計点でウィルコクソンの符号付き順位和検定をするか、ですね。t検定はあまりよろしくないかと思います。


二値データと連続変数の相関関係について
  投稿者:成瀬 2017/04/21(Fri) 16:11 No. 22325
初心者です。
基本的な質問で申し訳ございません。
二値データと連続変数の相関関係をみるには、二値データをダミー変数に置き換えて数値化して解析すると教わったのですが、これは正しい方法なのでしょうか?
二値データを数値化したところで、数値に意味がないわけなので、正しい方法とは思えないと考えています。

二値データと連続変数の相関関係をみる方法をご教示お願い申し上げます。

Re: 二値データと連続変数の相関関係について
  投稿者:課長 2017/04/22(Sat) 10:14 No. 22326
仰るとおり、1,0に意味は無いかと思います。計算上便利なだけでしょう。しかし、0,1を取るデータの数は意味があるのではないでしょうか?実際の計算では0や1を取る値の数、そして連続データの平均などの値を使うと思います。点双列相関係数やスピアマンの相関係数でいいんじゃないでしょうか?ちょっと自身がありませんので他の方の回答を期待します。

Re: 二値データと連続変数の相関関係について
  投稿者:青木繁伸 2017/04/22(Sat) 21:14 No. 22327
正しい方法です。
名前もあります。
点双列相関係数といいます。

しかし,それは,簡便化した計算方法を用いるものとして,「特別に?」定義づけられたものに過ぎず,実体は計算方法も含め「ピアソンの積率相関係数」にほかなりません。

「名義尺度を数値データとみて計算する」のは,うわべだけ見れば,「そんな馬鹿なことないよなあ」と思うかも知れませんが,名義尺度も二値データの場合は特別です。0/1 に限らず 1/2 であろうと50/1000 であろうと,標準化してしまえば同じ(標準化というのが,文字通り標準化ですから)。

相関係数を計算するには,間隔尺度,比尺度であること,というのがよく知られていて,名義尺度なんて間隔尺度でも比尺度でもないので,相関係数なんか計算できないよなあと思われるのでしょう。
しかし,間隔尺度の定義を考えてください。「データの間隔に意味がある」ですね。1と2の間隔と16と17の間隔は同じということです。
特別な場合として,データが2通りの値しか取らない場合,「1と2の間隔」しかないわけで,「じゃあ,だめだろ」といわれるかもしれないけど,「間隔尺度の取りうる値の個数」なんて,どこでも規定されているわけではなく,数学的には「間隔尺度の取りうる値の個数」なんて関係ないのです。二つしか取りうる値がないものも,間隔尺度なんですよ。

蛇足ですが,二値データは,簡単な一次変換でダミー変数になります。

なお,二値データを連続変数と同列に見て統計解析をするのは当たり前のように行われています。ダミー変数をつかっての重回帰分析などは広く行われています。
なお,k ≧ 3値データの場合も k-1 個のダミー変数に展開されて重回帰分析を初めとして多変量回帰に広く使われています。

二値データを数値化するというのには,明らかに意味があるのです。つまり,2つのカテゴリーがあれば,その二つは何らかの意味で大小関係があるのです。「何らかの意味」というのは,そのデータとそのデータを用いた分析の土俵上で決まるものですけど。

黒/白,男/女,高学歴/そうではない, 低収入/高収入,外国人/日本人,ある薬を投与された/されていない/,生産年齢人口/老年人口 ...
倫理的な問題も関与する場合もあるけど,なんだかんだいっても「少なくとも両者には何らかの意味で違いがある」ということは事実なんですから。

「統計学関連なんでもあり」の「全文検索」で「点双列相関係数」を検索すると,いろいろ興味深い記事を読めると思います。

Re: 二値データと連続変数の相関関係について
  投稿者:成瀬 2017/04/23(Sun) 12:46 No. 22328
課長様、青木先生

コメンどうもありがとうございました。
論文を読んでいてもニ値データと連続変数の相関解析の結果を良く目にしていたので、
解析方法に疑問を頂いていました。
2つしか取りうる値がないものも間隔尺度という考えは自分には全くなかったので、
学習になりました。

「点双列相関係数」で検索してさらに勉強させて頂きます。
ありがとうございました。

Re: 二値データと連続変数の相関関係について
  投稿者:課長 2017/04/23(Sun) 15:33 No. 22330
>成瀬様
中途半端な回答で申し訳ありません。

>青木先生
勉強になりました。今後とも宜しくお願い致します。


Cox回帰分析における独立変数の数とイベント数について
  投稿者:BOZ 2017/04/06(Thu) 11:23 No. 22314
連投ですみません。

Cox回帰分析で癌の再発に関連する臨床所見を検討したいと思っています。症例数は172例(うち再発が28例)、臨床所見(独立変数)は107個と非常に多いです。

RでAICによる変数選択(変数増減法)を行い、最終的に4個の独立変数がモデルに選ばれました。しかし以前この掲示板でご紹介いただいた下記の論文によると、イベント数は独立変数の10倍以上ないといけないとされています。ですので私のデータの場合、28/10=2.8なので、独立変数は多く見積もっても3個までとなってしまいます。

https://www.ncbi.nlm.nih.gov/pubmed/8970487?dopt=Abstract&holding=f1000

そこでお聞きしたいのですが、
4個の独立変数を含むモデルは、私のデータの場合ですと不適切でしょうか?
不適切としますと、例えば4個のうち最後に投入された変数1個を除き、独立変数を3個に絞ったモデルとするのは適切でしょうか?

Re: Cox回帰分析における独立変数の数とイベント数について
  投稿者:青木繁伸 2017/04/07(Fri) 22:46 No. 22315
3 はよくて 4 はだめだなんてことがないのは,明らかでしょう。
科学的に考えましょう。
選んだ 3 個も不適切かも知れないと疑うのが先でしょう。

Re: Cox回帰分析における独立変数の数とイベント数について
  投稿者:BOZ 2017/04/07(Fri) 23:13 No. 22316
青木先生、ありがとうございます。
モデルが適切か不適切かはどのようにして確認すればよろしいでしょうか?
モデルが他の集団にも適応可能かどうかバリデーションをすればよろしいですか?

Re: Cox回帰分析における独立変数の数とイベント数について
  投稿者:青木繁伸 2017/04/08(Sat) 19:23 No. 22317
> モデルが適切か不適切かはどのようにして確認

先行研究の知見を総合的に考慮する...

Re: Cox回帰分析における独立変数の数とイベント数について
  投稿者:BOZ 2017/04/08(Sat) 22:53 No. 22318
先行研究の結果と矛盾がなければ、仮にEPV(イベント数÷独立変数の数)≧10を満たさないモデルであっても、不適切とは限らないという事でよろしいでしょうか?

EPV≧10は多変量解析で守るべき指標なのか、あくまで参考程度の指標であるのか、その辺りが統計初心者の私にはよくわからず、悩んでおります。

Re: Cox回帰分析における独立変数の数とイベント数について
  投稿者:BOZ 2017/04/11(Tue) 11:55 No. 22324
EPVが5〜9程度でもいいのではないかという論文もあるようです。
https://www.ncbi.nlm.nih.gov/pubmed/17182981

あまりEPV≧10に固執しないで解析しようと思いました。


ロジスティック回帰分析(AICによる変数選択)の結果について
  投稿者:BOZ 2017/04/06(Thu) 10:36 No. 22313
ロジスティック回帰分析で腫瘍の良性/悪性に関連する画像所見を検討したいと思っています。症例数は70例(良性34例、悪性36例)、画像所見(独立変数)は21個です。

RでAICによる変数選択(変数増減法)を行ったのですが、下記のようにP値が0.05を上回る因子もモデルに入ってしまいます。この場合、有意でなかったFactor_CとFactor_D
を除いたFactor_A+Factor_Bを最終的な良悪性鑑別のモデルとするべきでしょうか?それとも4因子すべて含めるべきでしょうか?
Call:
glm(formula = x[, 1] ~ Factor_A + Factor_B + Factor_C + Factor_D, data = x)

Deviance Residuals:
Min 1Q Median 3Q Max
-0.93851 -0.22333 0.04443 0.19089 0.64343

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.427e+00 9.105e-01 -1.568 0.121822
Factor_A 8.127e-06 2.115e-06 3.843 0.000279 ***
Factor_B 2.625e-01 9.344e-02 2.810 0.006544 **
Factor_C -1.078e+00 6.143e-01 -1.755 0.083930 .
Factor_D -9.557e-01 5.801e-01 -1.648 0.104269
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.1024062)

Null deviance: 17.4857 on 69 degrees of freedom
Residual deviance: 6.6564 on 65 degrees of freedom
AIC: 45.947

Number of Fisher Scoring iterations: 2

Re: ロジスティック回帰分析(AICによる変数選択)の結果について
  投稿者:鈴木 康弘 2017/04/09(Sun) 09:13 No. 22319
 Factor_CとFactor_Dを含めないモデルだとAICが45.947より増えてしまう、ということですか?
 だったら4因子とも含めていいんじゃないでしょうか。(あまり自信なし)

Re: ロジスティック回帰分析(AICによる変数選択)の結果について
  投稿者:BOZ 2017/04/09(Sun) 12:36 No. 22320
>>Factor_CとFactor_Dを含めないモデルだとAICが45.947より増えてしまう、ということですか?
そうです。
4因子含めた場合、P値が0.05を上回るのはどう解釈したらよいのでしょうか・・・。

Re: ロジスティック回帰分析(AICによる変数選択)の結果について
  投稿者:青木繁伸 2017/04/09(Sun) 20:36 No. 22321
天はみずからたすくものをたすく
自己責任ってこと
それが研究者が研究者たること
その判断が広く受け入れられるかどうか,研究者としてやっていけるかどうかということ
だれかひとりのお墨付きをえられるかどうということではない
それじゃなきゃ論文一切を提示して私を共著者とする?私がレジェクトすれば,投稿をあきらめる?

Re: ロジスティック回帰分析(AICによる変数選択)の結果について
  投稿者:鈴木 康弘 2017/04/10(Mon) 07:10 No. 22322
AICはモデルの尤度から算出します。
個々の変数が有意かどうかはWald検定だと思います。
尤度比検定とWald検定の結果が一致しない場合、普通は尤度比検定の方を信用するみたいです。

Re: ロジスティック回帰分析(AICによる変数選択)の結果について
  投稿者:BOZ 2017/04/10(Mon) 12:37 No. 22323
青木先生、鈴木様

コメントありがとうございます。
自分の責任にて、4因子すべて含めたモデルで論文投稿してみます。
この度は大変ありがとうございました。

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

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