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

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

[トップに戻る] [利用上の注意] [ワード検索] [過去ログ] [統計学自習ノート] [管理用]
latest article: No. 22847, 2019/09/20(Fri) 23:21
おなまえ
タイトル
コメント
ＵＲＬ
添付ファイル
暗証キー (英数字で8文字以内)
自分の投稿記事を修正・削除する場合 --> 処理 記事No 暗証キー

複数のヒストグラムの取り扱いについて
　　投稿者：平沢 2019/09/20(Fri) 16:35 No. 22841
青木先生

初めての投稿になります。
大学の研究で行き詰ってしまったためご教授お願い致します。

試験体に含まれた薬剤量と、その薬剤の主成分が試験体のどこにどれほど入っているかを分析しています。データとしては、1つの試験体をおおよそ25000分割し(試験体毎に分割数は異なる)、その1点1点で主成分濃度を得ています。
例えば、試験体Aでは薬剤量149kg/m3、主成分濃度のサンプル数26510個といったものです。

以上のようなデータが30試験体分あり、それぞれで主成分濃度のヒストグラムを作成致しました。ヒストグラムの階級数はスタージェスの公式を用いて算出しています。
ヒストグラムの形について、主観では、薬剤量が 1)一定値以下 2)一定値付近 3)一定値以上の3タイプに分類されるように見えておりますが、これを客観的に言える方法が無く困っております。
2)一定値付近では比較的正規分布しているように見えたためコルモゴロフ＝スミルノフ検定を行いましたが、帰無仮説は棄却されました。

稚拙な説明で申し訳ございませんが、ご助言いただければ幸いです。
Re: 複数のヒストグラムの取り扱いについて
　　投稿者：青木繁伸 2019/09/20(Fri) 23:06 No. 22843
データのイメージがはっきりつかめないので，シミュレーションデータを作って分析してみました。

n 試験体（貴方の例では n = 30) のそれぞれについてかなりのデータがあり，それを class （貴方の例では class = 25000） 区間に分割してヒストグラム（度数分布）を作った。このヒストグラムについて，似ているもののグループ分けをしたい？
set.seed(123)
# class = 2500
class = 250 # ヒストグラムの階級数
N = 10000 # データ数
# n = 30
n = 15 # 試験体
data0 = matrix(0, n, N) # 大元のデータ（試験体数 x データ数）
data = matrix(0, n, class) # 試験体ごとの度数分布
Mean = runif(n, -1, 1) # 試験体ごとの平均値（取りあえず）
for (i in 1:n) { # 試験体ごとの元データ（取りあえず正規分布に従うとして）
  data0[i, ] = rnorm(N, mean=Mean[i])
}
Min = min(data0) # データ全体の最小値（度数分布を作るために）
Max = max(data0) # データ全体の最大値（度数分布を作るために）
for (i in 1:n) { # 試験体ごとの度数分布（ヒストグラムは描かない）
  data[i, ] = hist(data0[i,], breaks=seq(Min, Max, length=class+1), plot=FALSE)$counts
}
さて，このような，試験体 x 度数分布の集計結果をデータとみなして主成分分析を行うとどうなるか？
a = prcomp(data) # scale.=TRUE が必要かどうか？
主成分分析の結果を図示する。試験体の配置をみるために，主成分得点を図示する。
図をクリックすると拡大表示

次の発言に続く

Re: 複数のヒストグラムの取り扱いについて
　　投稿者：青木繁伸 2019/09/20(Fri) 23:09 No. 22844
続き

第1主成分の順番で度数分布を図示してみる。
plot(c(Min, Max), c(1, n), type="n", xlab="class", ylab="sample no.", yaxt="n")
t = order(a$x[,1])
for (i in 1:30) {
	points(seq(Min, Max, length=class+1), rep(i, class+1), pch = 19, cex=sqrt(data[t[i],] )/10)
	text(-4, i, t[i], col="red")
}
先の図からも似ているもの同士の位置関係がわかるが，実際どのように似ているかを試験体ごとの度数分布を示してみる。
第一主成分の順番で表示してみると，平均値の大きさにしたがって配置されていることがわかる。

Re: 複数のヒストグラムの取り扱いについて
　　投稿者：青木繁伸 2019/09/20(Fri) 23:15 No. 22846
第2主成分の順でプロットしてみる。
plot(c(Min, Max), c(1, n), type="n", xlab="class", ylab="sample no.", yaxt="n")
t = order(a$x[,2]) # 第2主成分を指定する
for (i in 1:30) {
	points(seq(Min, Max, length=class+1), rep(i, class+1), pch = 19, cex=sqrt(data[t[i],] )/10)
	text(-4, i, t[i], col="red")
}
右上，左上，下部の3パターンに分類できる。

Re: 複数のヒストグラムの取り扱いについて
　　投稿者：青木繁伸 2019/09/20(Fri) 23:21 No. 22847
結局の所，試験体ごとの度数分布をデータと見なして似たもの同士を集める統計手法を使えばどうだろうかということ。

因子分析，対応分析，クラスター分析でもよいだろうと思う。（私が理解したデータ構造があたっていれば）


VBAでの購買間隔の平均を計算
　　投稿者：波音 2019/09/20(Fri) 12:27 No. 22839
業務上の事情があってExcel VBAで「購買月間隔の平均」を計算するユーザー定義関数を作っています（作りました）。
自分でチェックしたところ、正しく動作はしているようですが。。。
記憶をたよりに書いていることもあり、何かもう少しスマートに書ける気がするようなしないような。。。

という状況ですので、もう少しスマートに書く方法があればご教示いただきたいです。

やりたいことは以下のようなものです。

例えば、私が毎月購買していれば4月〜3月までの12セルすべてに「1」が入ります。
もし購買していない月があれば「空白」が入ります。
仮に4,5,6；9,10；3に購買しなかった場合、空白は3, 2, 1ということになるので、その平均は(3 + 2 + 1) / 3という計算です。

この計算結果は2.00となりますが、意味としては「この人は平均的に2ヶ月間隔で買っていない(2ヶ月間隔でブランクがある)ということ。

VBAコード

Option Explicit

' 指定したセル範囲において、入力値と入力値の間の空白をカウントし、
' 　空白数の平均を計算して返す関数

' ID×日付の購買頻度クロス集計表で、IDごとの購買間隔の日数の平均を求めるために作った関数

Function betweenAve(datRng As Range)

Dim dat As Variant                  ' 関数の引数で取得したアドレスにある数値を格納する配列
Dim n As Integer, i As Integer      ' 配列の長さn、forの繰り返し用i

Dim blankCount As Integer           ' 空白セルをカウントして格納していくハコ
Dim sumN As Integer                 ' 購買日間隔がいくつかるかカウントして格納していくハコ

Dim mySum As Double, ans As Double  ' 空白セルを累計して格納するハコ

' まずは引数で指定されたセル範囲のアドレスから、セルに入力されている値をdatへ格納
' 二次元配列として格納されるので、行ベクトルの長さを取得
dat = Range(datRng.Address).Value
n = UBound(dat, 2)

' 各種変数を初期化しておく
blankCount = 0
sumN = 0
mySum = 0

' セルの値を1つずつチェックしていき、空白があればblankCountに「1」を格納していく
' 　それをmySumへ累計していく
For i = 1 To n
    If dat(1, i) = "" Then              ' i番目のセルが空白がどうかチェック
        blankCount = blankCount + 1     ' 空白なら「1」を代入
    Else
        If i <> 1 Then                  ' 空白でない場合、そのセルが最初のセルかどうかチェック
            If dat(1, i - 1) = 0 Then   ' 最初のセルでなければ、1つ前のセルをチェックして空白なら、
                sumN = sumN + 1         ' 　sumNに+1カウント
            End If
        Else
            '何も処理しない
        End If
        
        mySum = mySum + blankCount      ' mySumへ空白の数を累計していく
        blankCount = 0                  ' 1ループの処理が終わるたびにblankCountを初期化する
    End If
Next i

' 最終セルが空白の場合、ループ内で数がカウントされないのでループ後にカウントする
mySum = mySum + blankCount

' 最終セルが空白の場合、購買日の間隔数もカウントされないのでループ後にカウントする
If dat(1, n) = "" Then
    sumN = sumN + 1
End If

' 購買日数間隔の平均を計算
 ans = mySum / sumN

' 戻り値
 betweenAve = ans

End Function

Re: VBAでの購買間隔の平均を計算
　　投稿者：波音 2019/09/20(Fri) 12:29 No. 22840
Excelのワークシートのスクリーンショットも添付します。
このようなデータイメージです。



分散分析で大丈夫でしょうか
　　投稿者：コロン 2019/09/18(Wed) 08:06 No. 22833
お世話になっております。ご教授ください。

３クラスに対して、授業前後の意識の変化を5件法で問うています。項目は5項目あり、それらの平均値を特性値にします。

この分析手法ですが、２元配置分散分析でよろしいでしょうか？それとも、五件法ということて、ノンパラの手法、フリードマンで行うべきでしょうか？

よろしくお願いいたします。
Re: 分散分析で大丈夫でしょうか
　　投稿者：コロン 2019/09/18(Wed) 09:06 No. 22834
追加で書かせてください。

5項目の平均値を特性値にすると書きましたが、足し算した値を特性値にすることも考えましたが、大きな違いはないのでしょうか。

基本的な質問で申し訳ございません。
Re: 分散分析で大丈夫でしょうか
　　投稿者：青木繁伸 2019/09/18(Wed) 19:05 No. 22835
> 平均値を特性値にすると書きましたが、足し算した値を特性値にすることも考えました

5 で割るか，そのままかで，scale の違いだけなのでどちらも同じ結果になるのでは？
Re: 分散分析で大丈夫でしょうか
　　投稿者：コロン 2019/09/18(Wed) 19:18 No. 22836
青木先生

ありがとうございます。このデータは、どの手法を用いて差があるかどうかを見ればよろしいのでしょうか。
Re: 分散分析で大丈夫でしょうか
　　投稿者：青木繁伸 2019/09/19(Thu) 08:14 No. 22837
5 件法でしたね。だとすると合計を取るのも平均値をとるのも不適切。
ということで，フリードマン検定か。
Re: 分散分析で大丈夫でしょうか
　　投稿者：コロン 2019/09/19(Thu) 10:13 No. 22838
青木先生

お返事ありがとうございます。ということは、一つの概念を表していると考えている5項目を、それぞれにフリードマンという理解でよろしいでしょうか。


数量化�粁爐砲弔�まして
　　投稿者：こだま 2019/09/11(Wed) 17:30 No. 22830
青木先生

お世話になります。
こだまと申します。

本掲示板の　2019/08/14(Wed) 13:13 No. 22795　にて
統計分析の手法（判別分析・数量化�粁燹砲砲弔い銅遡笋気擦討い燭世い燭發里任后�

先生からのご助言やその他の参考書等を読み、最終的には数量化�粁爐砲栃�析を進めた
のですが、ご教示いただきたい事があり、再度書き込みさせていただきます。

数量化�粁爐砲董�公園のあり/なしを目的変数として、その要因を分析したのですが、分析に使用したサンプル数が、あり：144、なし：5　でした。
相関比や判別的中率については、それぞれ0.49、96.0％となり、精度としても良さそうなので、この分析に基づいて文章を書こうかと思っていたのですが、今回のようにサンプル数に偏りがある場合に問題はないのだろうか（分析するデータとしてはよろしくない？）
と思いました。

ある書籍では、「数学的にはサンプルが3以上あればよい」というような記載していたのですが、本当にそれで良いのか、分析上問題がある場合はどのような分析ができそうかについてご教示いただけますと幸いです。
Re: 数量化�粁爐砲弔�まして
　　投稿者：青木繁伸 2019/09/11(Wed) 21:08 No. 22831
「数学的にはサンプルが3以上あればよい」 という，「サンプル」が何を意味するか？そもそも，統計学で「サンプル」という用語を使うだけで，怪しさ満載ですが。

しっかりした分析プログラムで分析できれば，「数学的に問題なし」と理解してよいでしょう。「サンプル」がなにかを置いておいても，「サンプルの3個が全く同じ値を取ったら」数学的には問題ありというか，分析プログラムはエラーで停止するでしょう。

それはさておき，「数学的に問題なし」と「実質的に妥当性あり」とは違います。例えばあり144，なし0なら，数学的に問題あり。あり144，なし1なら「数学的に問題なし」ですが，「妥当性がある」とは言えないでしょう。あり144，なし2なら？あり144，なし3なら？...
どのあたりから妥当性があるといえるでしょうね？これは統計学で答えることはできない問題でしょう。大方の人が妥当だと思ってくれないと，論文としてアクセプトするのは難しいと言うことになるでしょう。
Re: 数量化�粁爐砲弔�まして
　　投稿者：こだま 2019/09/11(Wed) 21:33 No. 22832
青木先生

お世話になっております。

御返事いただきありがとうございます。

ご指摘いただいたように、色々と当方の専門とする分野の論文をレビューしている中で、統計の面で本当に妥当なのか？というものがいくつかあり、妥当性というところでもやもやしておりました。

ここで質問してお答えいただく事ではなかったにも拘らず、ご教授いただきありがとうございます。大変勉強になりました。


教育心理　重回帰　ロジスティック
　　投稿者：佐藤 2019/07/01(Mon) 16:11 No. 22770
教育心理分野で、ある実習後の学習意欲を目的変数とし、実習中の環境要因２つ、実習生の個人要因２つを説明変数として重回帰分析をしようと考えていました。

学習意欲に対して正規性の検定を行ったところ正規性が認められなかったため、重回帰分析を適応できませんでした。

1. この場合、ノンパラメトリックの相関分析までしか適応できないのでしょうか？

2. 学習意欲のデータを見たところ二峰性があったため、学習意欲の中央値で2群に分割し、高学習意欲群1と低学習意欲群0にダミー変数化してロジスティック回帰分析をしても良いのでしょうか？（それともこれはデータの恣意的な変更や情報量の削減という観点で許されないことなのでしょうか？）

統計に明るくなくおかしなことを聞いているかもしれませんが、ご教授頂けますと助かります。
Re: 教育心理　重回帰　ロジスティック
　　投稿者：青木繁伸 2019/07/01(Mon) 21:10 No. 22771
> 学習意欲に対して正規性の検定を行ったところ正規性が認められなかったため、重回帰分析を適応できませんでした。

重回帰分析において，目的変数が正規分布するかどうかは関係ありません。というか，目的変数がどのような分布をするかによって，適切な分析方法があります。

例えば，一つ前の質問のような，目的変数が 0/1 の二値変数の場合にはロジスティック回帰分析とか，その他にも glm 関数は目的変数が様々な場合に対して分析を行うことができるように対処しています。
binomial(link = "logit")
gaussian(link = "identity")
Gamma(link = "inverse")
inverse.gaussian(link = "1/mu^2")
poisson(link = "log")
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")
また，単に重回帰分析の場合であっても，従属変数が正規分布に従うことを要件にしてはいません。
Re: 教育心理　重回帰　ロジスティック
　　投稿者：佐藤 2019/07/02(Tue) 14:18 No. 22772
青木先生

ご返信有難うございます。
不勉強でした重回帰分析は従属変数の正規分布を前提としていないんですね。
ネット上でK-S検定を行って正規性を確認してから、
重回帰分析を行っている例があったので、重回帰分析は正規分布を仮定している統計法と誤解をしていました。

従属変数は間隔尺度のデータなので、無理に2値化してロジスティック回帰分析をするよりも、
重回帰分析を行なうほうがデータに沿った分析方法のように考えられます。

重回帰分析もう一度勉強し直してみます。

青木先生、今回は誠にありがとうございました。
Re: 教育心理　重回帰　ロジスティック
　　投稿者：佐藤 2019/08/23(Fri) 11:48 No. 22812
青木先生

以前、丁寧にご回答頂きまことにありがとうございます。
以前の質問の続きなのですが、目的変数、説明変数ともに正規性の無いデータで重回帰分析を行なったところ、学会にて正規性の無いデータで重回帰分析は使用できないとの指摘を受けました。

「すくできる！リハビリテーション統計」という書籍や他医療統計の書籍に
「重回帰分析では独立変数も従属変数も正規分布した数量データ（間隔尺度か比率尺度）である必要があります。」

という記述があります。
以前のやり取りの中で、正規分布しないデータでも重回帰分析が適用できるという、
当方の理解が間違っていたのでしょうか？
Re: 教育心理　重回帰　ロジスティック
　　投稿者：青木繁伸 2019/08/23(Fri) 14:21 No. 22813
その本の著者も，学会（査読者？）も間違えています。間違ったことを書いているクソページもたくさんありますね。

前に説明したことに基づいて，反論しなかったのですか？

説明変数は正規分布しなくても差し支えありません。ダミー変数を説明変数に使う場合を考えて見ればわかるでしょう。ダミー変数は 0/1 の二値データで，正規分布などしません。カテゴリーデータは重回帰分析に使えない（使うなら数量化 I 類を使いなさい）などという，とんでもない時代錯誤の考えを持っている人もまだいるのだろうか？また，実験計画によるデータの場合，説明変数は特定の値を取ります（薬剤投与量と効果の場合，薬剤投与量は例えば 10mg, 20mg, ..., 100mg などとされるでしょう。だれも薬剤投与量を正規乱数から選ぼうなどとは思わないでしょう）。

目的変数は，正規分布しなくても構いません。正規分布しないといけないのは残差の分布です（さらに，それさえも必須条件ではない）。
誤差が正規分布しないと，予測値の信頼区間を求めることができないからです。誤差（測定誤差も含む）は普通は正規分布に従うものです。なので，統計分析できるのです。

日本語のページは信頼できないということなら，英語のページを検索してはいかがでしょうか。

https://www.statisticssolutions.com/assumptions-of-multiple-linear-regression/
Multiple regression assumes that the residuals are normally distributed.

http://www.restore.ac.uk/srme/www/fac/soc/wie/research-new/srme/modules/mod3/3/index.html
Normally distributed residuals: The residuals should be normally distributed.

https://www.m3.com/open/clinical/news/article/604122/
回帰式の形は、被説明変数のタイプによって変えることができます。被説明変数が連続変量の場合は重回帰分析、2値変数（時間依存性がないイベントの発生または非発生、死亡または生存、治癒または非治癒など）の場合はロジスティック回帰分析、打ち切りのある2値変数の場合はCox回帰分析を用います。

重回帰分析では、残差（回帰分析による予測値と実測値の差分）の分布が正規分布である必要があります。

https://bellcurve.jp/statistics/course/9700.html
回帰モデルを考えるにあたって、誤差 μi にはいくつかの仮定している条件があります。
1. μi の期待値は0である
2. μiの分散は常にσ2 である
3. 異なる誤差 μi、μj は互いに独立である
これらの3条件から、μi は互いに独立に正規分布 N(0, σ2) に従うと仮定されます。

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

具体例でも示してみましょう。正規分布すべきは「残差」です。

理論モデル y = 2.5 + 4.8 * x

重回帰モデルは y = a + b * x + ε
εが正規分布

データを用意する

n = 1000
x = runif(n, min=0, max=100) # 独立変数
EPSILON = rnorm(n, mean=0, sd=10) # ε：正規分布
y = 2.5 + 4.8 *x + EPSILON # 誤差を含む従属変数

y は正規分布しない！！

hist(y, main="図1")

予測は十分

plot(y ~ x, pch=19, col="#aa000055", alpha=0.3, main="図2")
ans = lm(y ~ x)
summary(ans)
abline(ans)
text(10, 400, sprintf("y = a + b x\na = %.3f\nb = %.3f", ans$coef[1], ans$coef[2]), pos=4)

誤差（予測誤差）は正規分布する

plot(ans$residuals ~ x, pch=19, col="#aa000055", main="図3")

hist(ans$residuals, main="図4")

図は，クリックすることで拡大表示

Re: 教育心理　重回帰　ロジスティック
　　投稿者：佐藤 2019/08/28(Wed) 09:41 No. 22827
青木先生

出先で返信書き込み遅くなり大変申し訳ございません。
詳細な回答いただき誠にありがとうございます。

当方で医療関係の統計処理の書籍を確認したところ、

1. 目的変数が正規分布している必要がある。
2. 目的変数、説明変数とも正規分布している必要がある。

という記載の本が非常に多く見られました。
医療系の人間は統計をあくまでツールとしてしか使用していないので、
学会でも上記の説が主流となっているようです。
正しい記述の本も、正規分布している必要がないものをわざわざ「正規分布している必要はない」と記載するわけではないので、上記の説が蔓延していると理解しました。

２値化してロジスティック回帰すると明らかに情報量の損失が発生しますので、
青木先生にご教授頂いたことを用いて、反論をしてみます。

ご丁寧にご教授頂き誠にありがとうございました。


計算関係にあるデータ間の相関検定の可否
　　投稿者：コダマ 2019/08/25(Sun) 17:48 No. 22816
身長(A)と足の長さ(B)のデータが300人分あります。
胴長度(C)=1-B/Aで計算されるとき
BとCの相関を一次回帰分析で検定することは妥当でしょうか

背が高いほど胴長率が大きいのかを知りたいと思っているのですが、考えてみたら、ひとつひとつのデータ自体は一次計算式の代入変数と計算結果という関係にあるわけです。検定の可否をご教授いただけると幸いです。
Re: 計算関係にあるデータ間の相関検定の可否
　　投稿者：青木繁伸 2019/08/25(Sun) 20:57 No. 22820
> ひとつひとつのデータ自体は一次計算式の代入変数と計算結果という関係

データ変換が一次変換であれば同じですが，この場合は一次変換ではないので，結果は異なるでしょう。
layout(1:2)
A = c(150, 160, 180, 175)
B = c(80, 75, 95, 82)
plot(B ~ A)
C = 1-B/A
plot(B, C)
layout(1)
predict(lm(B ~ A))
predict(lm(B ~ C))
の結果
> layout(1:2)
> A = c(150, 160, 180, 175)
> B = c(80, 75, 95, 82)
> plot(B ~ A)
> C = 1-B/A
> plot(B, C)
> layout(1)
> predict(lm(B ~ A))
       1        2        3        4 
75.71429 80.19780 89.16484 86.92308 
> predict(lm(B ~ C))
       1        2        3        4 
87.51340 78.87029 86.76991 78.84639 







Re: 計算関係にあるデータ間の相関検定の可否


　　投稿者：コダマ  2019/08/26(Mon) 18:06  No. 22825   


ありがとうございました。大いに参考にさせていただきます。











どう扱えば良いか困っています








　　投稿者：りゅうへい  2019/08/25(Sun) 12:53  No. 22814   


※答えを求めているわけでは断じてありません。ご理解下さい。
以下のデータが得られたとして　行う分析とそれにより得られる結論を述べよというものです。

会社Ａ 資本金1000万 従業員100名
業績平均値：総務部７０ 営業部５０ 開発部５０

会社Ｂ 資本金3億 従業員400名
業績平均値：総務６０ 営業３０ 開発６０
いずれの平均値に対してもＳＤ＝5とする

分散分析を行うのだとは思うのですが，会社間では資本金と人数が異なるので
それをどう分析の中で考慮すればいいのかが分かりません。

ご教授のほどよろしくお願いします。







Re: どう扱えば良いか困っています


　　投稿者：青木繁伸  2019/08/25(Sun) 17:25  No. 22815   


業績平均値は各部の部員の業績の平均値ということでしょうかね？
資本金の違いは関係ないのでは？
全従業員数の違いも関係ないと思いますが，総務部，営業部，開発部 の人数がわからないので，分析不可能では？
つまり，会社全体の業績平均値を知ろうにも，各部の人数がわからなければ計算しようがない。

「これだけのデータではどのような分析もできない。問題自体が不適切という結論が得られる」が答えかな？







Re: どう扱えば良いか困っています


　　投稿者：りゅうへい  2019/08/25(Sun) 18:36  No. 22818   


ご回答ありがとうございます。
問題には　仮想の業績得点の平均値と書かれていました。
平均値の認識は　ご指摘の通りで合っていると思います。
3部署合計人数＝従業員数なのか　3部署以外の職員を含め　従業員数なのかが記載されていませんでした。　
3部署合計が全従業員数　という前提であれば　各部の人数が分からずとも
会社全体の平均値は　Ａは　(70+50+50)/3=56.67 B (60+30+60)/3=50　
と　求めることは可能でしょうか？

資本金の違いは　結論が出た際の背景要因として触れる項目なのでしょうか？

ご教授のほどよろしくお願いします。







Re: どう扱えば良いか困っています


　　投稿者：青木繁伸  2019/08/25(Sun) 20:40  No. 22819   


> 3部署合計が全従業員数　という前提であれば　各部の人数が分からずとも会社全体の平均値は　Ａは　(70+50+50)/3=56.67 B (60+30+60)/3=50 と　求めることは可能でしょうか？

不可能です。平均値の意味をわかっていれば，不可能であることは自明です。

たとえば，会社 A の合計人数が 100 人であっても，
総務部，営業部，開発部がそれぞれ，20, 30, 50 の場合は (70*20 + 50*30 + 50*50) / (20+30+50) = 54 だし，
22, 33, 45 の場合なら (70*22 + 50*33 + 50*45) / (22+33+45) = 54.4 だし，
1, 2, 47 の場合なら (70*1 + 50*2 + 50*97) / (1+2+97) = 50.2 だ。
こんな例を挙げる必要があるとは思わなかった（^_^;）

(70+50+50)/3=56.67　が成り立つのは， 総務部，営業部，開発部がそれぞれ，100/3, 100/3, 100/3 と等しいとき（あり得ないが）

(70*100/3 + 50*100/3 + 50*100/3)/(100/3 + 100/3 + 100/3) = 56.66666666







Re: どう扱えば良いか困っています


　　投稿者：青木繁伸  2019/08/25(Sun) 21:00  No. 22821   


> 資本金の違いは　結論が出た際の背景要因として触れる項目なのでしょうか？

たとえ他の条件が適切に示されていたとしても，資本金の違いがこのような違いをもたらしたかどうかについては，情報が「全く不十分」なので，何ともいえません。

この問題，インターネット上で提示されてるものですか？それなら，URL を示してください。







Re: どう扱えば良いか困っています


　　投稿者：りゅうへい  2019/08/25(Sun) 21:15  No. 22823   


すみません　データ数が等しいケースだと勝手に解釈しておりました。
宿題で出された問題です。











どのような統計をつかうのでしょうか








　　投稿者：こういち  2019/08/18(Sun) 03:28  No. 22802   


青木先生
以下の疑問に関して教えてください

「疾患A型の患者」および「疾患B型の患者」における左室形態（4つタイプがありNormal, concentric remodeling(CR), concentric hypertrophy(CH), eccentric hypertrophy(EH))の分布を明らかにしました。結果の全体像は表1のようになりました。

ここで私が示したいことはA型ではCRが有意に多く、B型ではCHが有意に多い、
ということです。

そのために私がどのようにしたかというと
表2、表3のような2x2表をつくって（この作り方が正しいのか疑問ですが）
Fisher法で計算をしてそれぞれp=0.0427, p=0.0378という結果を得ました。
まず質問はこの統計処理が正しいのかを教えてください。

次にこれがもし正しいとしたら、下位比較をするという意味で、たとえば有意水準をBonferroniによって0.05/4C2=0.0083以下、とするのでしょうか？

別の考え方として4x2のクロス集計をおこなってP=0.1169を得たとして、
有意差なしでA型とB型の左室の型の分布に差がなし、で終了でしょうか？

以上、教えてください
よろしくお願いいたします









Re: どのような統計をつかうのでしょうか


　　投稿者：青木繁伸  2019/08/19(Mon) 10:11  No. 22803   


最初から CR，CH 型の比較だけを行うと決めていたのなら，(CR, CH)×(A,B) の 2×2 分割表の検定をすればよいのですが，今回はそうではなさそうなので，全ての組み合わせでの 2×2 分割表の検定をしてボンフェローニ法を適用する。(CR, CR 以外)×(A,B)，(CH, CH以外)×(A,B)などとするのではないので注意。つまり，NormalとCR, NormalとCH，NormalとEH, CRとCH，CRとEH，CHとEH  の 6 通りの比較を行う。
なお，「全体で有意である場合に下位検定として行える」という条件が必要かどうかは諸説ありますが，「この条件は不要」との見解が優勢かもしれない。

Fisher の正確検定でなければ，http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Tukey.html に書いてある方法を使えるでしょう（http://aoki2.si.gunma-u.ac.jp/lecture/Hiritu/Pmul-Ryan.html でもよいが，前者の法が通りがよいかな？）。







Re: どのような統計をつかうのでしょうか


　　投稿者：こういち  2019/08/20(Tue) 00:56  No. 22804   


青木先生

お返事をいただきありがとうございました。
responseが遅くなりお詫び申し上げます。
先生からいただいた説明を何度も何度も読み直してみましたが
私の説明が悪くて質問の意図が伝えられていないと思いました

私の疑問は
「たとえばＡ型という疾患においてCR, CH型の頻度を比較をする」ということではなく
「4つある形態の一つであるCRという形態を呈する患者が、Ａ型では全体の56%, B型では全体の38%で、これらを比較するときに、どの様な統計を用いれば、CRという形態はA型がB型よりも有意に多いことを示せるか？」というA型とB型の比較のための統計処理がわからない、ということです。

そして同様に上記と独立して、
「4つある形態の一つであるCHという形態を呈する患者が、Ａ型では全体の34%, B型では全体の52%で、これらを比較するときに、どの様な統計を用いれば、CHという形態はB型がA型よりも有意に多いことを示せるか？」というA型とB型の比較ののための統計処理がわからない、ということです。

私のいっている比較そのものがそもそも統計処理にそぐわず、有意に多い、ことは示せないのでしょうか？。
基本的なことが分かっていないとんちんかんな質問かもしれないと、
だんだん心配になってきましたが、どうか教えていただければ、本当にありがたいです。
よろしくお願いいたします。







Re: どのような統計をつかうのでしょうか


　　投稿者：青木繁伸  2019/08/20(Tue) 08:28  No. 22805   


No. 22802 
> 私が示したいことはA型ではCRが有意に多く、B型ではCHが有意に多い、ということです。

と

No. 22804 
> 「たとえばＡ型という疾患においてCR, CH型の頻度を比較をする」ということではなく...CRという形態はA型がB型よりも有意に多いことを示せるか

相反する説明がありますが。まあ，それは行方向のパーセントと列方向のパーセントのどちらを評価対象にするかということで，どちらが原因で，どちらが結果なのかで自ずと決まるものですが。でも，chisq.test でも，Fisher.test でも行・列の違いは無視されます（結果は同じ）。

> 「4つある形態の一つであるCRという形態を呈する患者が、Ａ型では全体の56%, B型では全体の38%で、これらを比較するときに、どの様な統計を用いれば、CRという形態はA型がB型よりも有意に多いことを示せるか？」というA型とB型の比較のための統計処理

ならば，二項検定でよいでしょう。これだと，行ごとの検定になります。CR が A:66, B:19 で帰無仮説では母比率=0.5
> binom.test(66, 85)

	Exact binomial test

data:  66 and 85
number of successes = 66, number of trials = 85, p-value = 3.041e-07
同様に CH では

> binom.test(40,66)

	Exact binomial test

data:  40 and 66
number of successes = 40, number of trials = 66, p-value = 0.1089
ただし，この二つの検定は独立ではないです（独立性の検定という意味ではない）。
CR の 66 も 29 も Normal, CH, EH の度数による制約があります。

(Normal,CR,CH,EH)x(A,B) はそれらの制約のもと，Normal, CR, CH, H それぞれで，A, B の割合に違いがあるかを検定するものですよね。
検定ファミリーの一つに (CR,CH)x(A,B) があり，それは，CR, CH で A, B の割合に違いがあるかを検定する。
> fisher.test(matrix(c(66, 40, 19, 26), 2))

	Fisher's Exact Test for Count Data

data:  matrix(c(66, 40, 19, 26), 2)
p-value = 0.03111
CR では A，CH では B が有意に多い（A では CR，B では CH  が有意に多い）といえるわけでしょう。

なお，また最初に戻りますが，「A型ではCRが有意に多く、B型ではCHが有意に多い」ということなら，A, B ごとに一様性の検定
> chisq.test(c(7, 66, 40, 5))

	Chi-squared test for given probabilities

data:  c(7, 66, 40, 5)
X-squared = 86.407, df = 3, p-value < 2.2e-16

> chisq.test(c(2, 19, 26, 3))

	Chi-squared test for given probabilities

data:  c(2, 19, 26, 3)
X-squared = 34, df = 3, p-value = 1.981e-07
でも，これはこの場合は Normal と EH がCR, CH より格段に少ないことが原因。
> binom.test(66, 106)

	Exact binomial test

data:  66 and 106
number of successes = 66, number of trials = 106, p-value = 0.01478

> binom.test(19, 45)

	Exact binomial test

data:  19 and 45
number of successes = 19, number of trials = 45, p-value = 0.3713
こんないろいろなことを考えなきゃいけないけど (CR,CH)x(A,B) でやればよいのでは？







Re: どのような統計をつかうのでしょうか


　　投稿者：こういち  2019/08/21(Wed) 21:43  No. 22806   


青木先生

丁寧に教えていただきまして、本当にありがとうございます。
教えていただいたことが恐らくは十分に理解されていないと思いますので
時間をかけて考えてみようと思います。

最後にひとつだけ教えてください、
先生が薦めてくださる「(CR, CH)x(A, B)をやる」、というのは
下のような２ｘ２の表を作成してp=0.0311を求め、
（この表の作り方が正しいのかわかりません）
この場合「A型ではCRが多く、B型ではCHが多い」
（←この言い方が正しいかどうかにも自信がありません）ことが
0.0311>0.008(=0.05/6)であり有意とは言えない

という理解でよろしいでしょうか？
最後に教えていただければ本当にありがたいです









Re: どのような統計をつかうのでしょうか


　　投稿者：青木繁伸  2019/08/22(Thu) 08:13  No. 22807   


> A型ではCRが多く、B型ではCHが多い
> 0.0311>0.008(=0.05/6)であり有意とは言えない

でよいと思います。







Re: どのような統計をつかうのでしょうか


　　投稿者：こういち  2019/08/22(Thu) 18:18  No. 22811   


青木先生

お忙しい中いろいろ教えていただき
本当にありがとうございました
また、よろしくお願いいたします











3Dグラフの作成方法について








　　投稿者：さくらとお城  2019/08/09(Fri) 14:19  No. 22791   


青木先生

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

　複数のチャンネルで計測した経時データを，plot3Dなどで図示しようとしているが，なかなかうまく行きません。

　データとエクセルで作成したグラフは添付ファイルのようになります。実際は，もっとチャンネルが多いので，ここでは2つだけ表示しています。

　きれいに図示できる方法を教えてください。
Download:22791.zip 22791.zip








Re: 3Dグラフの作成方法について


　　投稿者：青木繁伸  2019/08/09(Fri) 20:20  No. 22792   


3次元グラフはあまりよい表示方法ではないでしょう。
添付する図を見ても，二次元グラフ（折れ線グラフ）のほうが二つのデータの相違・同等がよく分かるでしょう。
確かにたくさん重なるとわかりにくいでしょうが，3次元折れ線グラフが重なったら，なにがなんだかわからなくなります。









Re: 3Dグラフの作成方法について


　　投稿者：青木繁伸  2019/08/09(Fri) 20:30  No. 22793   


ggplot の例によくあるような，縦にずらっと折れ線グラフを並べるのはどうですか？
一つのグラフに，傾向の似通った複数の折れ線グラフを描けばなおよいかも。









Re: 3Dグラフの作成方法について


　　投稿者：さくらとお城  2019/08/10(Sat) 10:56  No. 22794   


青木先生

　ご返事ありがとうございます。添付のようなグラフを考えていました。ちなみに，これが別ソフトで作ったもので，細かい設定などが大変で，Rでもできるかなと思って質問しました。

　先生のggplotの例も良さそうで，試してみます。









Re: 3Dグラフの作成方法について


　　投稿者：さくらとお城  2019/08/22(Thu) 10:22  No. 22808   


青木先生

　その後，色々と調べたけど，恥ずかしながら，なかなか先生のような図を描けません。ggplotの例を書いたコードを見せて頂いて宜しいでしょうか？

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







Re: 3Dグラフの作成方法について


　　投稿者：青木繁伸  2019/08/22(Thu) 12:50  No. 22809   


私自身，ggplot は使わないので，細かな設定は知りません。
以下のようなプログラム例を示しておきますので，追加，修正の上ご利用ください。
ポイントは通常のデータフレームを時間とチャンネルとデータ値の3列から成るデータフレームに変換して，それを用いて折れ線グラフを描くということです。
set.seed(123)

m = 10
n = 50

# n x m 行列：列ベクトルがチャンネルごとの時系列データ

simdata0 = matrix(0, n, m)
for (i in 1:m) {
	simdata0[, i] = cumsum(rnorm(n))
}

# チャンネルごとのデータを，
# time:時間（1:n），value:データ値，チャンネル名の，n x 3 のデータフレームとして rbind で縦に繋いだデータフレームに変換

simdata = NULL
for (i in 1:m) {
	simdata = rbind(simdata, data.frame(time = 1:n, value = simdata0[, i], 
		channel = paste0("x", i)))
}

g = ggplot(simdata, aes(x = time, y = value))
g = g + geom_line(aes(colour = channel), size = 0.5)
g = g + geom_point(aes(colour = channel), size = 1, alpha = 0.5)
g = g + facet_wrap(~channel, ncol = 2) # ncol で何列のグラフにするか
g = g + theme_classic() # または g = g + theme_bw() でも
g
以下のようなグラフになります（部分）









Re: 3Dグラフの作成方法について


　　投稿者：さくらとお城  2019/08/22(Thu) 12:53  No. 22810   


青木先生

　早速のご返事ありがとうございます。先生のコードを勉強させていただきます。











統計分析手法の選択につきまして








　　投稿者：こだま  2019/08/14(Wed) 13:13  No. 22795   


青木先生

はじめまして。

自身の周りに統計手法を専門的に理解している人間があまりおらず書き込みさせて
いただきます。

私は現在、大学院生（博士後期課程）に在学中で現在GISを用いた緑地の分布に
ついての分析結果学を会論文投稿に向け作業をしているところです。

その中で、現在の緑地の分布要因は何に規定されているのか？を統計手法による
分析を用いて明らかにしたいと考えております。

その際にどの分析手法が良いのかに行き詰っているところです（判別分析なのか、
数量化2類なのかと悩んでいます。）。

データはメッシュで集計しており，メッシュ毎に以下のデータを格納しています。

目的変数：緑地の分布タイプ（別途分析により7タイプに分類済）
説明変数：用途地域（13タイプ），区画整理事業の有無，傾斜角度，標高，駅の有無等

目的変数が質的データ、説明変数が質的データと量的データが混在しているため、単純
に判別分析や数量化2類を使えない？のかと思っています（理解不足なのかもしれませんんが）。

お忙しいところお手数をおかけしますが、どうぞご回答いただけますと幸いです。







Re: 統計分析手法の選択につきまして


　　投稿者：青木繁伸  2019/08/14(Wed) 13:20  No. 22796   


> 目的変数が質的データ、説明変数が質的データと量的データが混在しているため、単純
に判別分析や数量化2類を使えない？のかと思っています

の点についてのみ。
質的データはダミー変数として取り扱うことで，量的データと一緒にして判別分析に使用できます。（ちなみに，数量化II類はダミー変数を使う判別分析と等価です）







Re: 統計分析手法の選択につきまして


　　投稿者：こだま  2019/08/14(Wed) 13:58  No. 22797   


青木先生

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

判別分析を採用してみたいと思います。

標高等の量的なデータを区分して数量化2類かな？などと考えて
いて、その区分の方法が難しいなと思っていたところでした。

今後、また相談させていただくことがあるかもしれませんが、
どうぞよろしくお願いいたします。







Re: 統計分析手法の選択につきまして


　　投稿者：こだま  2019/08/14(Wed) 17:12  No. 22798   


青木先生

度々申し訳ございません。

色々な文献（多変量解析の本やインターネット、判別分析を行っている研究論文）を
調べているのですが、多くが目的変数が2択（例えば、公園がある/ない のような者）
の事例でした。

今回、分析しようとしているように目的変数が7つのタイプがあっても分析手法として
採用しても問題はないのでしょうか。

可能であればで結構ですので、ご教授いただけますと幸いです。







Re: 統計分析手法の選択につきまして


　　投稿者：青木繁伸  2019/08/14(Wed) 21:32  No. 22799   


ロジスティック判別も線形判別も，3群以上であっても判別はできます。
どこでもよいですが，判別分析のページを読んでみてください。
一番簡単な2群判別の話しか書いていないのは，くそページです。

7群判別ともなると理論的には可能ですが，実際に分析して得られる解の解釈は難しいかも。（貴方の用意した説明変数で十分に判別できるかどうかということですが）







Re: 統計分析手法の選択につきまして


　　投稿者：こだま  2019/08/14(Wed) 23:06  No. 22800   


青木先生

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

論文提出まで時間があるので判別分析含め、その他の多変量解析について
色々な参考書やページを読んで勉強いたします。

説明変数の選択が難しいところですが、こちらは既往文献等を参考にしつつ
選択していきたいと思います。

実際に分析を進めてみて何かありましたら相談させていただきたいと思います。

夜分遅くまでお手数をおかけしました。











【R】リストに格納された複数データフレームのEXCEL保存








　　投稿者：明石  2019/08/05(Mon) 16:12  No. 22787   


青木先生　様；

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

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

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

リストに格納された複数データフレームを、
拡張パッケージ openxlsx を用いて、
EXCELの複数のシートに格納したxlsxファイルで保存することを考えます。

以下、irisデータでお示しをします。

data(iris)

dfs <- split(iris, iris$Species)

library(openxlsx)　　#　拡張パッケージ

write.xlsx(dfs, file="iris_Species.xlsx")

setosa、versicolor、virginica　の３つのシートからなる
"iris_Species.xlsx"ファイルは保存されますが、以下のメッセージがでます。

「Note: zip::zip() is deprecated, please use zip::zipr() instead」

私が書いたコマンドには、zip()関数はありませんので、
このようなメッセージが表示された場合、どのように対処すればよいのでしょうか。

ご教示をいただけましたら、大変に助かります。
お手数をおかけいたします。







Re: 【R】リストに格納された複数データフレームのEXCEL保存


　　投稿者：青木繁伸  2019/08/05(Mon) 19:40  No. 22788   


write.xlsx から呼ばれる saveWorkbook の中で zip を呼ぶようですが（明示的には見えない），私にはそれ以上追求できませんでした。
deprecated というのは，将来的に廃止されるということで，無効になっているわけではないです。ということで，しばらくはそのまま使えば良いと思います。そのうち，使えなくなったらまずは openxlsx パッケージのバージョンを確認して，新しくなっていたらインストールし直すということになるでしょう。それもできないという事態になったら，パッケージ開発者に連絡するしかなくなるかな。







Re: 【R】リストに格納された複数データフレームのEXCEL保存


　　投稿者：青木繁伸  2019/08/05(Mon) 19:43  No. 22789   


Excel シートにする必要性がどの程度あるのか，他の人と共有するためには Excel シートでなければならないとか，いきなり Excel を書き込むのではなく，CSV ファイルとして書き出すのではだめなのか，とか，私なら Excel を使わない方向で行きたいですね。







Re: 【R】リストに格納された複数データフレームのEXCEL保存


　　投稿者：明石  2019/08/06(Tue) 08:05  No. 22790   


青木先生　様；

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

今回も有難いご教示をいただき、大変に助かりました。

青木先生がお書きになられたことを、戒めたいと思います

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











【R】リストに格納された複数データフレームのファイル保存（ループを使わない）








　　投稿者：明石  2019/07/29(Mon) 17:25  No. 22784   


青木先生　様；

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

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

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

リストに格納された複数データフレームのファイル保存を、
ループを使わない書き方について、ご教示いただければ大変に助かります。

以下、例示いたします。

データフレーム dat を、カラム名"支店名"でファイル分割したとします。

dfs <- split(dat, dat$支店名)

支店名ごとに分割されたデータフレーム群を、csvファイル保存します。

n <- length(dfs)
names <- names(table(dat$支店名))

for (i in 1:n) {
　fn <- paste(names[i],".csv", sep="")
　write.csv(dfs[[i]],fn, row.names=F)
}

ループを回さないで、apply関数群などを用いて、シンプルに書けるようならば、
ご教示をいただきたいと思います。

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







Re: 【R】リストに格納された複数データフレームのファイル保存（ループを使わない）


　　投稿者：青木繁伸  2019/07/29(Mon) 22:16  No. 22785   


iris データセットでやってみましょう

Species で 3 分割

iris2 <- split(iris, iris$Species)

書き出す。junk への代入は，ゴミを表示させないためだけ。

junk <- sapply(names(iris2), function(n) write.csv(iris2[[n]], paste(n, ".csv", sep="")))

あまりわかりやすいものではないですね。
貴方が示した for ループを使うほうがずっと分かりやすい。







Re: 【R】リストに格納された複数データフレームのファイル保存（ループを使わない）


　　投稿者：明石  2019/07/30(Tue) 06:13  No. 22786   


青木先生　様；

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

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

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











因子分析の因子の差について








　　投稿者：コロン  2019/07/16(Tue) 13:38  No. 22778   


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

早速ではございますが，「4つの因子」が出てきたという例で説明させて頂きます。

数学力が高いもの，中程度のもの，低いものという3つのグループの間で，因子において差があるかどうかを見る場合なのですが...

（1）因子毎に，一要因分散分析（被験者間）を行い，差があれば多重比較で見ていく
（2）4つの因子をすべて用いて，二要因分散分析（被験者「間」×被験者「内」）を行い，必要に応じて下位検定を行う

どちらが正しいのでしょうか？古い論文を見ていると（1）が圧倒的に多いのですが...。

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







Re: 因子分析の因子の差について


　　投稿者：青木繁伸  2019/07/17(Wed) 09:06  No. 22779   


理論的には (2) が正しい（妥当）と思います

なぜ (1) が圧倒的に多いかは，いくつか理由があるでしょう

サンプルサイズが小さく，(2) に耐えない
モデルが複雑になり，解釈が難しくなる
コンピュータプログラムがなく（利用しにくく），計算できない
「だって，みんな (1) をやっているんだもの」と長いものに巻かれる
そもそも (2) という選択肢を知らない







Re: 因子分析の因子の差について


　　投稿者：コロン  2019/07/17(Wed) 09:36  No. 22780   


青木先生

お忙しい中，いつもありがとうございます。

やはり，そうですよね....。かなり検索をしてみたのですが，すべて（1）で処理されていたのです。（2）が適切ではないかなとずっと考えていましたので，すっきり致しました。

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







Re: 因子分析の因子の差について


　　投稿者：コロン  2019/07/20(Sat) 17:19  No. 22781   


青木先生

もう一つ重要なことを聞き忘れていました。

各因子の項目の数が異なるのですが、尺度得点あるいは因子得点を利用して２要因分散分析を行ってもよろしいでしょうか。

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







Re: 因子分析の因子の差について


　　投稿者：青木繁伸  2019/07/24(Wed) 11:04  No. 22782   


項目の数が異なるというのが二元配置分散分析の障害になるのですか？







Re: 因子分析の因子の差について


　　投稿者：コロン  2019/07/25(Thu) 16:34  No. 22783   


青木先生

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

【要因１】
因子１　７項目
因子２　５項目
因子３　４項目

【要因２】
学力上位群
学力下位群

項目数が異なるため，各因子の平均値をそのまま解釈して良いのかなと悩んでおりました。ひょっとしたら，過去の研究，項目数が異なるために因子毎に比較していたのかなと考えたり...。

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











固辞，数字の置換








　　投稿者：コロン  2019/07/08(Mon) 13:44  No. 22773   


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

基本的な質問で申し訳ございませんが，以下のデータがあるとします。

class	score
1	10
1	11
2	13
2	14
3	12
3	15
4	17
4	14
5	13
5	10

classの 3 を 5 に修正して新しいデータフレーム形式で表示させるにはどのようにすればよろしいのでしょうか？







Re: 固辞，数字の置換


　　投稿者：青木繁伸  2019/07/08(Mon) 15:37  No. 22774   


データフレームを a とすれば，
a$class[a$class==3] <- 5
で置換できます。







Re: 固辞，数字の置換


　　投稿者：コロン  2019/07/08(Mon) 16:40  No. 22775   


青木先生

お忙しい中申し訳ございません。

ネットで検索をしていたら，置換は gsub というサイトが多かったため，gsubでやっていたのですが，何度やっても，データフレーム形式にならず，かつ，すべての数字が５に変換されたので，わからなくなり，投稿した次第です。

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







Re: 固辞，数字の置換


　　投稿者：青木繁伸  2019/07/08(Mon) 18:07  No. 22776   


gsub でもできないことはないですが，gsub は本来は文字列の置換ですから
a$class <- as.numeric(gsub(3, 5, a$class))
as.numeric してやるなんて余計なことをしないといけません。







Re: 固辞，数字の置換


　　投稿者：コロン  2019/07/08(Mon) 18:24  No. 22777   


青木先生

as.numericなんですね！gsubでもできました。ただ先生から最初に提示いただきましたものがわかりやすいです。

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











hclust()の距離行列とメソッドの指定








　　投稿者：波音  2019/06/28(Fri) 15:03  No. 22764   


Rのhclust()について教えてください。

引数にward.Dとward.D2が指定できるようになっていますが、例えば「d」という距離行列(ユークリッド距離)が存在したと時にいくつかの方法がネット上や書籍で示されています。

・距離行列をそのまま指定する
hclust(d, method="ward.D")

・距離行列の自乗したものを指定する
hclust(d^2, method="ward.D")

・川端ほか「Rにおる多変量解析入門」オーム社（2019）P374記載
hclust((1/2)*d^2, method="ward.D")

結局のところ「解釈ができれば何でもいい」というくらいの問題なのかもしれませんが、アルゴリズム的（？）にはどういう指定が適切なのでしょうか。

川端ほか「Rにおる多変量解析入門」オーム社（2019）P313にward.Dとward.D2の違いが説明されていますが、どういう距離行列とメソッドの組み合わせ指定をすればよいのでしょうか、、、







Re: hclust()の距離行列とメソッドの指定


　　投稿者：青木繁伸  2019/06/29(Sat) 12:05  No. 22765   


川端さんが具体的にどのように説明されているのか分かりませんが，私の理解は
http://aoki2.si.gunma-u.ac.jp/lecture/misc/clustan.html
に示してあります。これは奥野忠一らの多変量解析，多変量解析入門にあったものだと思います。まだパソコンなんてものがなく，コンピュータはFORTRAN,COBOLでプログラムするという時代のものです。もともとは，CLUSTAN というプログラムのアルゴリズムの説明にあるものです。実際，私もこれに基づいてクラスター分析のプログラムを書いて使っていました。
当初の R には今の ward.D が ward として実装されていました。推測するに，「バグ」だったのだと思います。
私の書いた本でも，「ウォード法と重心法を採用する場合にはユークリッド平方距離を使うべきである。ユークリッド平方距離を選択できるオプションがないので，後述の hclust の第 1 引数としてユークリッド距離を二乗したものを指定することを忘れないようにしよう。」と書きました。
その後指摘があったかなんかで，本来あるべき ward.D2 が実装されましたが，過去のしがらみを捨てきれなくて（バグだったとも言えなくて）両方を残したのだと邪推しています。







Re: hclust()の距離行列とメソッドの指定


　　投稿者：波音  2019/07/01(Mon) 10:41  No. 22769   


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

たしかに青木先生の本にも「hclust の第 1 引数としてユークリッド距離を二乗したものを指定」と書かれており、それ以前にもこのように教わっていたので、私の中ではある意味ルール的にd^2としていましたが、、、

やはり（使う人の立場によっては程度問題ではあると思いますが）アルゴリズムを理解しておくことは重要だと感じさせられました。

少し話はそれてしまいますが、誰でも使えるものというのは「考えなくても（理解していなくても）使えるもの」ともいえるので、ただでさえ色々と誤解の多い統計話題は今後さらに注意しないといけないと思いました。

※パッケージ文化に恵まれた我々は特に注意せねばならないのかもしれません。











ダミー変数を従属変数にした場合の切片








　　投稿者：山本  2019/06/29(Sat) 23:47  No. 22766   


はじめて書き込みます。山本です。
講義の関係で重回帰分析を行なっているのですが、ダミー変数を従属変数とした際に切片が負という結果がでました。
ダミー変数は0か1かですので正になるはずなのですが…なぜでしょうか
この結果はどう解釈すればよいのでしょうか
それとも説明変数を変更する必要があるのでしょうか
お願いいたします







Re: ダミー変数を従属変数にした場合の切片


　　投稿者：青木繁伸  2019/06/30(Sun) 07:19  No. 22767   


従属変数が 0/1 の二値データであることと，切片が負の値を取ることは何の矛盾もありません。
というか，独立変数のとる値により切片は負の値になることも，[0, 1] の値になることも，1より大きな値になることもあります。
たとえば，
d <- data.frame(y = c(1, 0, 1, 1, 0, 0, 0, 1, 0, 0),
	x1 = c(58, 37, 38, 50, 36, 59, 49, 66, 47, 58),
	x2 = c(54, 53, 50, 46, 46, 65, 37, 69, 37, 44))
の場合は，切片は負の値になります。
> lm(y~x1+x2, data=d)
(Intercept)           x1           x2  
  -0.646540     0.005094     0.015825  
ここで，x1 から 250 を引いた変数 x1.250 を作り，それを用いて予測式を計算すると
> d$x1.250 <- d$x1 - 250
> lm(y~x1.250+x2, data=d)
(Intercept)       x1.250           x2  
   0.627073     0.005094     0.015825  
ここで，x1 から 435 を引いた変数 x1.435 を作り，それを用いて予測式を計算すると
> d$x1.435 <- d$x1 - 435
> lm(y~x1.435+x2, data=d)
(Intercept)       x1.435           x2  
   1.569546     0.005094     0.015825  
それぞれの偏回帰係数は全く同じになっていることに注意。
つまり，独立変数の平行移動をした場合，その影響を打ち消すために定数項が調整されるということ（ちなみに独立変数を定数倍したときは偏回帰係数が定数分の1になる）。

で，今回の問題の本質は，なぜ従属変数が 0/1 の二値データなのに予測値が0より小さくなったり，1より大きくなったりしてしまうのかということ。

それは直線回帰をするからやむを得ないということです。
図を見れば明らかです。

二値データはロジスティック回帰分析を行うのが定石です。
図の赤で示したのがロジスティック曲線です。
> glm(y~x1+x2, data=d, family=binomial(link="logit"))

Call:  glm(formula = y ~ x1 + x2, family = binomial(link = "logit"), 
    data = d)

Coefficients:
(Intercept)           x1           x2  
   -5.08999      0.02083      0.07214  

Degrees of Freedom: 9 Total (i.e. Null);  7 Residual
Null Deviance:	    13.46 
Residual Deviance: 11.93 	AIC: 17.93









Re: ダミー変数を従属変数にした場合の切片


　　投稿者：山本  2019/06/30(Sun) 12:02  No. 22768   


素早く回答していただきありがとうございます。
そういった分析もあったのですね…
よく分かりました











介入研究の因子分析








　　投稿者：春  2019/06/13(Thu) 00:34  No. 22746   


青木先生

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

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

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







Re: 介入研究の因子分析


　　投稿者：春  2019/06/17(Mon) 15:24  No. 22758   


青木先生

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

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

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







Re: 介入研究の因子分析


　　投稿者：青木繁伸  2019/06/17(Mon) 19:45  No. 22759   


上の方にあるものばかり見ていたようで，気づきませんでした。

さて，先行研究で使われていた尺度は因子分析で作成されたものではないと言うことでしょうか。

とすると，まずいですね。あなたの調査データで一因子性がないという結果になっても，先行研究を否定するということになってしまいますし，対象が違うから違う因子が得られたのだといわれるかも知れませんし。

介入研究だからといって，因子分析のやり方に特徴や注意点があるとは思えませんが。

介入前後で因子が異なることは想像するに難くないと思いますよ。

実際の因子分析の結果だけでなく，基礎的な統計解析を見ずにああだこうだいうのは難しいです。もう少しはっきりいえば，あなたの持っているデータを自由にいろいろ分析してみてその結果を総合して何かをいうということでなければ（貴方の分析結果を見ただけではなく），なんもいえない。ということです。そういうことをするのは，共同研究者以上の立場ですよね。査読者はある意味そのような立場にあるのでしょうが，実際に「データを見せろ，私が解析してみて貴方の言うことが正しいかどうか示してやる」とは言えない訳です。

データ解析というものは，そういうものです。







Re: 介入研究の因子分析


　　投稿者：青木繁伸  2019/06/17(Mon) 21:33  No. 22760   


「十分な分析を行い，以下のような結果を得ました」というとき，どのように十分な分析だったかを「詳細に」示し，それを査読者が十分だったと認めるかどうかということです。

蛇足ではありますが，査読者が十分だったと認め論文掲載となっても，読者（第三者）の十分だったかどうかの判断とは必ずしも一致しない（投稿論文への異議申し立てはあり得る）。

その昔，ある学会誌に「奇病」の投稿があったが，老臨床医が「それは，脚気だろ」と指摘してチャンチャン！となったことがある（そんなに大昔のことではない）

厳しく言えばそういうこと。







Re: 介入研究の因子分析


　　投稿者：春  2019/06/18(Tue) 00:41  No. 22761   


青木先生

ご助言ありがとうございます。

そもそも、因子分析されていない尺度を使う事自体、無理があるということですね。
もう1つ、気がかりなことがあります。Nの数が30に満たないので、それも、因子分析するには、少なすぎるかなあと思っておりました。

このまま、因子分析を行うことは、諦めた方が良いのかとも思いますが、以下のように、進めるのは、間違いでしょうか？

介入前の因子分析の結果では、4因子に分かれ、α係数も1つだけ0.7に足りませんでしたが、3つは0.8以上あり、内的整合性は高いと判断できます。
介入後も、介入前に抽出された因子に合わせて、α係数を出して充分な値であれば、その後に介入前と介入後で対応のあるｔ検定をする。

・・・というのは、分析としてはどうでしょうか？
ご意見を伺えると助かります。







Re: 介入研究の因子分析


　　投稿者：青木繁伸  2019/06/18(Tue) 08:38  No. 22762   


> Nの数が30に満たないので、それも、因子分析するには、少なすぎるかなあと思っておりました。

そこは最初に書いておいてくれなくては。
圧倒的に少なすぎます。後で述べる対応のある t 検定であっても，十分とは言いがたい。

> その後に介入前と介入後で対応のあるｔ検定をする。

何の対応のある t 検定ですか？先行研究の尺度の合計値ですか？それとも，4因子それぞれに含まれる項目の合計点ですか？

そもそも，因子分析の結果の4因子はそれぞれ何項目が含まれるのか，それぞれの因子の意味づけ（解釈）はできるのか。

もともと20項目あるなら，項目別に検定を行ったらどうなるのか？







Re: 介入研究の因子分析


　　投稿者：春  2019/06/18(Tue) 12:13  No. 22763   


青木先生

おはようございます。
ご回答いただいて、ありがとうございます。

項目別の検定は、既に行っていて、効果は確認済みです。
しかし、ご指摘いただいたように、Nの数が少ないので、今回は因子分析は諦めようと思います。
しかし、今後、Nを増やして再度、分析を行いたいと思っているので、このまま、続けさせていただきたいと思います。

対応のあるt検定については、抽出した各因子の下位尺度得点として、介入前、介入後の各平均値を算出し、比較をしようと思っています。

しかしながら、今後、因子分析を行うときには、各因子内の項目の数や、各因子の意味づけ（解釈）に気を付けていきたいと思います。
そして、因子分析を用いた方が良いのか、項目別に検定だけの方が良いのか判断しようと思います。

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











【R】do.call()関数を用いたリスト（複数のデータフレーム）のmerge








　　投稿者：明石  2019/06/16(Sun) 09:50  No. 22754   


青木先生　様；

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

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

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

リスト（複数のデータフレーム）のmergeを、forループで順次mergeするのではなく、
do.call()関数を用いて、もし書けるようならば、ご教示をお願いいたします。

do.call(rbind, list)　にならい、
do.call(merge, list)　のように、もし書けたとしたら、どのように書けるのか興味があります。

# tempフォルダ下にある、csvファイルのファイル名を取得

path <- setwd("C:/temp") 
names <- list.files(path, pattern="csv")
nn <- length(names)

d <- list(nn)

loop <- 1: nn

for (i in loop) {
  d[[i]] <- read.table(names[i], header=TRUE, sep=",", stringsAsFactors=F) 
}

#　リスト d に格納されている nn個のデータフレームは、いずれも同じ、mergeの属性 "支店名"をもっています。
#　要素の添字に沿って、完全外部結合を順次、繰り返したいと思います。

for (i in loop) {
  if( i == 1) {
    dat <- d[[1]]
  } else {
    dat <- merge(dat, d[[i]], by.x="支店名", by.y="支店名" , all.x = TRUE, all.y = TRUE ,no.dups = F)
  }   
}   

もし、do.call(merge, d) のように書けるようならば、ご教示をいただきたいと思います。

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







Re: 【R】do.call()関数を用いたリスト（複数のデータフレーム）のmerge


　　投稿者：青木繁伸  2019/06/16(Sun) 12:34  No. 22756   


別に do.call でなくても，簡単に書きたいということなんでしょうか？

Reduce を使うのはいかが？
set.seed(11111)
d <- vector("list", 3)
d[[1]] <- data.frame(a=1:5, b=letters[1:5])
d[[2]] <- data.frame(a=c(1,3,5,6,8), c=LETTERS[1:5])
d[[3]] <- data.frame(a=c(2,4,5,7,9), d=rnorm(5))

merge2 <- function(x, y) merge(x, y, all.x = TRUE, all.y = TRUE, no.dups = FALSE)  
Reduce(merge2, d)
 一行で書くと
Reduce(function(x, y) merge(x, y, all.x = TRUE, all.y = TRUE, no.dups = FALSE), d)
なお，for を使う場合も，以下のようにすれば十分簡潔に書けるのでは？
x <- d[[1]] # 結果の初期化は，最初のデータフレームを代入
for (i in 2:length(d)) { # 2 番目以降のデータフレームをマージ
  x <- merge(x, d[[i]], all.x = TRUE, all.y = TRUE, no.dups = FALSE)
}
x
どのように書いても，以下のような結果になると思いますが。
  a    b    c            d
1 1    a    A           NA
2 2    b <NA>  0.003630422
3 3    c    B           NA
4 4    d <NA>  0.798122717
5 5    e    C  0.893397342
6 6 <NA>    D           NA
7 7 <NA> <NA> -1.195893897
8 8 <NA>    E           NA
9 9 <NA> <NA> -1.269878611








Re: 【R】do.call()関数を用いたリスト（複数のデータフレーム）のmerge


　　投稿者：明石  2019/06/16(Sun) 16:41  No. 22757   


青木先生　様；

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

先生がお書きになられたように、簡潔に書きたいというのが意図です。
そこで、do.call()関数に着目しました。

幾つかの方法をお示しくださいましたので、よく理解できました。
reduse()関数を初めて知りましたので、これを機会に勉強させていただきます。

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











低値・高値の基準について








　　投稿者：H  2019/06/14(Fri) 11:49  No. 22750   


青木先生
はじめまして。

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

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

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







Re: 低値・高値の基準について


　　投稿者：青木繁伸  2019/06/14(Fri) 21:04  No. 22751   


> 他の研究でもホルモン剤内服後のホルモン値が妊娠率に影響するかについての報告がありますが

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

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

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







Re: 低値・高値の基準について


　　投稿者：H  2019/06/15(Sat) 01:33  No. 22752   


返信ありがとうございます。
他の報告ではホルモン値が100毎に測定してあるものや平均値の低値と高値で評価してあるものがあります。ただそれらでは明確な根拠がないため恣意的になるかと思われます。

投降後にロジスティック回帰を行いましたが有意差を認めませんでした。
低値群で妊娠率が低いため、統計的に評価したいのですが恣意的になってしまいますか。
宜しくお願い致します







Re: 低値・高値の基準について


　　投稿者：青木繁伸  2019/06/15(Sat) 13:15  No. 22753   


> 投稿後にロジスティック回帰を行いましたが有意差を認めませんでした
ということなら，仕方ないでしょう。

有意ではないが医学的には意味があるようなら，奥の手で，もう少しデータを増やすとか。







Re: 低値・高値の基準について


　　投稿者：H  2019/06/16(Sun) 11:50  No. 22755   


返信ありがとうございます。
有意差なしでデータを出そうと思います。
ありがとうございました。











サンプルについて








　　投稿者：林檎  2019/06/12(Wed) 17:25  No. 22744  HomePage 


青木先生

はじめまして。
サンプルについてわからないことがあり、質問をさせていただきました。

私は、10施設の企業を対象に質問紙調査による横断研究を実施し、解析をしています。
500名の有効回答があり、施設ごとと従属変数にて一元配置分散分析を行ったところ、1施設のみ合計得点平均が有意に高いことが分かり、指導教授より、その施設は特殊であるため抜いたほうがいいのではないかとご指摘を受けました。ちなみにこの施設の対象は20名でした。
世の中の論文をみると、このような有意差を見ている論文はさまざまであり、今回先生のご意見をお聞かせいただけましたらと思っております。

宜しくお願い致します。







Re: サンプルについて


　　投稿者：青木繁伸  2019/06/12(Wed) 22:03  No. 22745   


> その施設は特殊であるため抜いたほうがいいのではないかとご指摘

そもそも，調査の目的は何だったのですか？というか，調査をして，何が分かると思っていたのですか？

「どの施設も同じようである」ということを知りたかったのなら，「1施設のみ合計得点平均が有意に高い」というのは，ビックリするような知見ではないですか。その施設は何らかの問題があると。

対象者が20名ということが特殊なのか，それ以外の条件が特殊なのか。いずれにしても，サンプルサイズが20であっても，それでも有意な差であると言うことですよ。

「その施設は特殊であるため抜いたほうがいい」ということなら，最初からその施設は対象としない方が良かったでしょう。不都合な結果が出たからなかったことにしようというのは，はなはだ科学的ではないです。

その施設を除いて，「どの施設も同じようでした」という結果を出して，何の意味があるのでしょうか？

指導教授とよく話し合うことをお勧めします。







Re: サンプルについて


　　投稿者：林檎  2019/06/13(Thu) 08:02  No. 22747  HomePage 


青木先生

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

目的としては、重回帰分析を行い、企業の社員のストレスの要因を明らかにすることです。
重回帰分析に入れるための変数を探すこと、属性を明らかにするために一元配置分散分析を実施したら、有意差があったという結果でした。
調査をお願いするときには、似たような企業だと思い調査を行ったのですが他の施設よりストレスが有意に高かったことが分かりました。

指導教授と相談をしてみます。
ありがとうございました。







Re: サンプルについて


　　投稿者：青木繁伸  2019/06/13(Thu) 10:32  No. 22748   


> 他の施設よりストレスが有意に高かったことが分かり

それこそ，その企業が他の企業とどこが違うか（ストレスの原因）を知るために使える重要なサンプルではないですか。







Re: サンプルについて


　　投稿者：林檎  2019/06/13(Thu) 14:00  No. 22749  HomePage 


青木先生

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

この企業は大切なサンプルであるということを改めて実感致しましたので、先生からのお返事をもとに解析を進めていきたいと思います。

とても丁寧なご対応をいただき、ありがとうございました。











クラスタリングされた個体の割合の比較とライアンの方法について








　　投稿者：PV  2019/06/09(Sun) 19:11  No. 22739   


青木先生、

はじめまして、
ご多忙のところ失礼いたします。

クラスタリングされた個体を群間で比較し、割合が有意に異なるかどうか、またどの群間でことなるか、検定を行いたいのですが、なかなか方法が見つからず調べていたところ、青木先生のホームページにたどり付きました。
以下記載いたしました質問 (1), (2) ご教授いただけましたら幸いです。

下の表は対象生物170個体をクラスター解析し、
各群(A〜F）で各クラスター（V1〜V9)に分配された個体数です。
	V1	V2	V3	V4	V5	V6	V7	V8	V9
A群	6	6	4	3	4	9	3	8	6
B群	2	1	8	4	0	1	1	2	3
C群	9	2	1	1	1	0	0	4	1
D群	8	2	3	2	0	3	4	2	3
E群	1	0	1	0	1	0	0	0	0
F群	3	3	4	9	6	4	4	7	10
質問(1)「クラスターに含まれる個体数の傾向は群間で異なるか」検定を行いたいのですが、
その場合、カイ二乗検定→ライアンの方法で多重比較　を行う方法は正しいでしょうか？

また、青木先生の説明およびＲコードを基に、カイ二乗検定で全体としての差が認められたことを確認後、ライアン法をひとまず行ってみたのですが、
エラーメッセージ「Error:  n > 0 は TRUE ではありません 」と表示され、検定を行うことができませんでした。
質問(2)ライアン法は0が含まれていると行えないのでしょうか？

(1), (2) についてご解答のほどよろしくお願いいたします。







Re: クラスタリングされた個体の割合の比較とライアンの方法について


　　投稿者：青木繁伸  2019/06/10(Mon) 09:47  No. 22740   


ライアン法は，あなたのデータには使えません（使用例参照）

全ての2群の組み合わせ（A:B, A:C, ..., E:F の15組〕について 2 x 9 の行列で chisq.test を行い，p 値をBonferroni法で解釈する（0.05/15)

手作業でやるのは大変なので，以下のようにすればよい
x = matrix(c(
  6, 6, 4, 3, 4, 9, 3, 8, 6,
  2, 1, 8, 4, 0, 1, 1, 2, 3,
  9, 2, 1, 1, 1, 0, 0, 4, 1,
  8, 2, 3, 2, 0, 3, 4, 2, 3,
  1, 0, 1, 0, 1, 0, 0, 0, 0,
  3, 3, 4, 9, 6, 4, 4, 7, 10),
  nrow = 6, byrow=TRUE)  

g = function(ij) {
  i = ij[1]
  j = ij[2]
  y= rbind(x[i,], x[j,])
  z = colSums(y)
  y2 = y[, z != 0]
  p = chisq.test(y2)$p.value
  cat(i, j, p, "\n")
  return(c(i, j, p))
}

gn = nrow(x)
cn = ncol(x)
p0 = 0.05 / choose(gn, 2)
res = data.frame(t(combn(6, 2, g)))
colnames(res) = c("i", "j", "p.value")
res$sig = ifelse(res$p.value > p0, "n.s.", "*")
print(res)
結果は以下のようになる
   i j     p.value  sig
1  1 2 0.060309971 n.s.
2  1 3 0.095874639 n.s.
3  1 4 0.365146072 n.s.
4  1 5 0.558144611 n.s.
5  1 6 0.384082770 n.s.
6  2 3 0.034933166 n.s.
7  2 4 0.234596686 n.s.
8  2 5 0.246627056 n.s.
9  2 6 0.177712096 n.s.
10 3 4 0.267205617 n.s.
11 3 5 0.438290867 n.s.
12 3 6 0.005751997 n.s.
13 4 5 0.160816574 n.s.
14 4 6 0.077873765 n.s.
15 5 6 0.441102843 n.s.







Re: クラスタリングされた個体の割合の比較とライアンの方法について


　　投稿者：PV  2019/06/10(Mon) 12:05  No. 22741   


青木先生、

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

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

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

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

猫も杓子も構造化
http://nekomosyakushimo.hatenablog.com/entry/2017/08/26/222728

のカテゴリー数を増やしただけだと思ったのですが、この解釈は間違っているということでしょうか？
統計解析が素人のため愚問で大変恐縮です。

今後解析を行う上でぜひ参考にさせていただきたいためご解答いただけましたら幸いです。
どうぞよろしくお願いいたします。







Re: クラスタリングされた個体の割合の比較とライアンの方法について


　　投稿者：青木繁伸  2019/06/10(Mon) 12:15  No. 22742   


ライアン法，チューキー法は，群（グループ）の数を k としたとき，k x 2 の行列データを対象にするのです。例えば，正否とか失敗・成功のような2値データが対象です。一方（もう一方）の「比率」が対象になります。

一般的に k x m 行列は，群により「分布が違うか」を対象にすることになります。

「カテゴリー数を増やしただけ」ではありません。







Re: クラスタリングされた個体の割合の比較とライアンの方法について


　　投稿者：PV  2019/06/10(Mon) 18:42  No. 22743   


青木先生、

なるほど、大変勉強になりました。
お忙しいところ大変丁寧にご説明頂き誠にありがとうございました。











R textでのイタリック表示








　　投稿者：コロン  2019/05/26(Sun) 17:39  No. 22733   


お世話になっております。

R plotで描画をしておりますが，グラフ内に p<0.01と表示させたいと思っております。

text(1.5, 70, labels="p <- 0.01)

pだけをイタリックで表示させたいのですが，どのようにすればよろしいでしょうか。

ご教示頂けますでしょうか。







Re: R textでのイタリック表示


　　投稿者：青木繁伸  2019/05/26(Sun) 21:08  No. 22734   


実行可能な最小限のテストプログラムを添付してください







Re: R textでのイタリック表示


　　投稿者：コロン  2019/05/27(Mon) 08:38  No. 22735   


青木先生

申し訳ございませんでした。簡単なコードではございますが，よろしくお願いいたします。

m1 <- 30
m2 <- 40
x <- c(1,2)
allm <- c(m1, m2)
plot(x, allm, ylim=c(0, 100))
text(1.5, 60, labels="p < 0.05")







Re: R textでのイタリック表示


　　投稿者：青木繁伸  2019/05/27(Mon) 18:13  No. 22736   


できばえはちょっといまいちなのですが（字が薄い），以下のようにすれば可能
m1 <- 30
m2 <- 40
x <- c(1,2)
allm <- c(m1, m2)
plot(x, allm, ylim=c(0, 100))
text(1.5, 60, labels="p", vfont = c("sans serif", "italic"))
x2 = 1.58 # この値を locator で決めるか，試行錯誤で決める
text(x2, 60, labels=" < 0.05", vfont = c("sans serif", "plain"))
または，書き始めの位置指定を文字列の左端として（pos=4），2 番目の文字列の最初を p のための空白を一つ確保し，
text(1.5, 60, labels="p", vfont = c("sans serif", "italic"), pos=4)
text(1.5, 60, labels="  < 0.05", vfont = c("sans serif", "plain"), pos=4)
とする（この方が試行錯誤不要なのでよいかも）。
"sans serif"  の代わりに "serif" も試して，良いと思う方を使ってください。







Re: R textでのイタリック表示


　　投稿者：コロン  2019/05/27(Mon) 20:52  No. 22737   


青木先生

お忙しい中，ありがとうございました。今作成しておりますグラフにも使用することができました。











Brunner-Munzel検定








　　投稿者：青木繁伸  2019/05/07(Tue) 18:23  No. 22724   


https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
でも紹介されているが，R の
lawstat パッケージの brunner.munzel.test() も
brunnermunzel パッケージの brunnermunzel.test() も
alternative="greater", "less" の場合の信頼区間がおかしい（"two.sided" の場合とまったく同じになっている）
nparcomp パッケージの npar.t.test で method="t.app" を指定して得られるのが正しい信頼区間のようだ。







Re: Brunner-Munzel検定


　　投稿者：青木繁伸  2019/05/09(Thu) 11:38  No. 22725   


> nparcomp パッケージの npar.t.test で method="t.app" を指定して得られるのが正しい信頼区間のようだ。

と書いたが，そもそも npar.t.test は片側検定で "greater" と "less" を取り違えているようだ。p 値が t.test と比べて逆になっている。これは困ったことになる。
df = data.frame(
  d = c(1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 2, 4, 1, 1, 3, 3, 4, 3, 1, 2, 3, 1, 1, 5, 4),
  g = rep(1:2, c(14, 11)))
library(nparcomp)
cl = 0.95
alt = "greater"
a = npar.t.test(d ~ g, data=df, round=5, info=FALSE, method="t.app", alternative=alt, conf.level=cl)
a$Analysis
b = t.test(d ~ g, data=df, alternative=alt, conf.level=cl)
b
を実行すると，npar.test の結果は
  Effect Estimator  Lower Upper       T p.Value
1 p(1,2)   0.78896 0.6291     1 3.13747 0.00289
だが，t.test の結果は
t = -2.9486, df = 15.916, p-value = 0.9953
p 値が 0.00289 と 0.9953 と 180 度異なっている。
alt="less" では 0.99711 と 0.004739 である。
ちなみに，wilcox.test の p 値は t.test と同一方向だ。 







Re: Brunner-Munzel検定


　　投稿者：青木繁伸  2019/05/10(Fri) 18:00  No. 22727   


引き続き

並べ替え検定については，奥村先生が
https://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html
でも独自のプログラム例を示しているが，
brunnermunzel パッケージには brunnermunzel.permutation.test がある。
しかし，これは場合によって正しい答えを出さない。というか，奥村先生のプログラムと結果が異なることがある。
そのようなデータの例は，x = c(36, 45, 51, 49, 57)，y = c(48, 42, 49, 39, 45) の場合。

奥村先生のプログラムでは 0.365079365079365 になるのだが，brunnermunzel.permutation.test では 0.341269841269841 になる。

前者は 92/252，後者は 86/252 だ。92 の内訳は，データを並べ替えたデータにより計算される検定統計量の絶対値が，観察値における検定統計量の絶対値と同じかより極端な場合の数 45*2 と，検定統計量が +Inf, -Inf になる 2 通りの計 92 である。

観察値における検定統計量の絶対値と同じになるのは 4*2 = 8 通りあるが brunnermunzel.permutation.test では，計算誤差のためか 8 個全てが同じとなっていないのであろう。つまり，数え落としている。

また，Inf, -Inf となる場合の扱いもチェックする必要があるだろう。
前に書いて消したが，検定統計量が Inf になるのは R の場合だけで，Python の場合には 0 による割り算エラーで計算が中止される。

いろいろ問題ありで，それぞれのプログラムどれも信頼して使うわけにはいかな状態である。

自分で書くしかないか。







Re: Brunner-Munzel検定


　　投稿者：青木繁伸  2019/05/10(Fri) 21:15  No. 22728   


Brunner-Munzel検定 最強！！

なんて，書かれているページが多くて，ですね。ちょっと，立ち止まって検証してみませんか？という趣旨で書きました。

ほかの人の意見も聞きたいので，「拡散希望」（^_^;）







Re: Brunner-Munzel検定


　　投稿者：it may be that  2019/05/19(Sun) 19:00  No. 22729   


npar.t.test の greater/less の扱いですが、取り違えているわけではないと思いますよ。
summary(a)
に
#---------------------------Interpretation---------------------------------------------#
 p(a,b) > 1/2 : b tends to be larger than a 
と出ます。
したがって、群1と2 (g1とg2) の観測量をX1, X2として、対立仮説は
'greater'はP(X1 < X2) > 1/2, 
'less'はP(X1 < X2) < 1/2
です。
Estimate が 0.78896 > 1/2 であることから、P(1,2)はP(X2 - X1 > 0)であると思われます。
よって、書き直すと
'greater'はP(X2 - X1 > 0) > 1/2,
'less'はP(X2 - X1 > 0) < 1/2 あるいは P(X2 - X1 < 0) > 1/2
となります。

一方、g1とg2の位置母数をそれぞれm1とm2とすると、t.testの対立仮説は
'greater' :m1 - m2 > 0
'less' : m1 - m2 < 0
です。
つまり、npar.t.testの対立仮説は
m2 - m1 > 0
m2 - m1 < 0
に相当し、基準が逆になっています。
確かに紛らわしいですが、lm, glmなどでcontr.treatmentに馴染んでしまった身としては、むしろt.testやwilcox.testのfactorの扱いに違和感を覚えます。
普段contr.SASを使っていればt.test,wilcox.testが自然でしょう。







Re: Brunner-Munzel検定


　　投稿者：青木繁伸  2019/05/21(Tue) 11:01  No. 22730   


「npar.t.test の greater/less の扱い」についてのコメント，ありがとうございます。
確かに，説明の全体を見ればその通りですね。







Re: Brunner-Munzel検定


　　投稿者：青木繁伸  2019/05/21(Tue) 11:04  No. 22731   


追加の考察

http://oku.edu.mie-u.ac.jp/~okumura/stat/brunner-munzel.html

では，p = Pr(x＜y)+0.5*Pr(x=y) ではなく，t 分布による漸近近似を使っている。

n1, n2 が小さいときに行う並べ替え検定で，n1, n2 が大きいことを仮定する漸近近似を行うという，相矛盾する状況である。

p そのままを使って並べ替え検定を行うと，
foo2 = function(X) {
  brunner.munzel.test(xandy[X], xandy[-X])$estimate
}
z2 = combn(1:N, n1, foo2)
bm2 = brunner.munzel.test(x, y)$estimate
mean(abs(z2-0.5) >= abs(bm2-0.5))
0.006690223 が得られる。

なお，この値の正当性としては，マン・ホイットニーの U 検定統計量と Brunner-Munzel の p との間に p=1-U/n1/n2 という単純な関係があるので，マン・ホイットニーの U 検定統計量を計算して使う並べ替え検定と同じ値になることを示せば十分。

または，coin パッケージの wilcox_test で正確な p 値を求めると，同じ値 0.006690223 が得られる。結局は wilcox_test が並べ替え検定をやっているということだ。
> df = data.frame(d=c(x, y), g=factor(rep(1:2, c(length(x), length(y)))))
> wilcox_test(d ~ g, data=df, distribution="exact")

    Exact Wilcoxon-Mann-Whitney Test

data:  d by g (1, 2)
Z = -2.6934, p-value = 0.00669
alternative hypothesis: true mu is not equal to 0
ということは，マン・ホイットニーの U 検定と Brunner-Munzel 検定って，前者が正規分布，後者が t 分布（など）で漸近近似を行うだけで，本質は同じということだ...
> npar.t.test(d ~ g, data=df, info=FALSE, method="perm", nperm=400000, round=5)$Analysis
       Estimator Statistic   Lower   Upper p.value
id       0.78896   3.13747 0.59708 0.99057 0.00708
logit    0.78896   2.38394 0.57650 0.91561 0.00694
probit   0.78896   2.51949 0.58109 0.92446 0.00694
となれば，奥村先生のページの冒頭の説明，

> 二つの確率変数 X1，X2 が同じ分布に従うという帰無仮説を検定するには，有名なWMW検定（Wilcoxon-Mann-Whitney test）が使えます。しかし，分布が異なる（例えば分散が異なる）場合には，これでは正確に検定できません。

> そこで，分布が同じことは仮定せず，両群から一つずつ値を取り出したとき，どちらが大きい確率も等しいという帰無仮説を検定することにします。

という，Brunner-Munzel 検定の優位性というのもなくなるのか。（前提が異なるといえども，検定統計量が同じなんだからねぇ）







Re: Brunner-Munzel検定


　　投稿者：青木繁伸  2019/05/24(Fri) 18:57  No. 22732   


検定統計量ではなく，その近似（t分布に近似する）のときに「ベーレンス・フィッシャー問題を考慮した」ということですね。











クロンバックαの係数








　　投稿者：wtcs  2019/05/05(Sun) 08:24  No. 22721   


青木先生

はじめまして。お忙しい中 またお休み中に失礼いたします。

昨年尺度開発を行い、現在はその開発尺度を用いて縦断調査にて影響要因を探索しています。
尺度開発の論文を執筆し投稿しましたところ、査読者から信頼性の検討のところでご指摘をいただきました。
尺度は全15項目で、下位因子は5項目、4項目、6項目の3因子です。
クロンバックαが0.89,0.77,0.70で全体では0.78でした。
それぞれの下位因子のクロンバックαが全体のクロンバックαより高いのはおかしいとご指摘をいただいたのですが、どこを調べてもわかりませんでした。文献も同じようなデータのものもあります。
確かに項目数が多い全体の方が高くなりやすいとは思いますが、必ずしも全体の方が高い値でないとその尺度の信頼性は低いのでしょうか？
尺度開発は初めてのため、初歩的な質問だと思いますがどうぞよろしくお願いいたします。







Re: クロンバックαの係数


　　投稿者：青木繁伸  2019/05/05(Sun) 15:00  No. 22722   


ある数学的（統計学的）命題が誤りであることを示すには，反例を一つあげれば十分ですね。

妥当なデータを，信頼のおけるプログラムで分析した結果がその命題に反するものであった場合，その命題は誤っているということでしょう。

今の場合，命題は「下位尺度のクロンバックのアルファは全体のクロンバックのα信頼性係数より小さい」ということ。

信頼のおけるプログラムは，たとえば R とか SPSS でいいでしょう。
妥当なデータというのはちょっとやっかいで，確かにそのようなデータがあるということでは不十分なのです。その統計手法が想定している条件を満たしているデータでなければいけません。

クロンバックのα信頼性係数は各質問項目への回答の不偏分散の和と各質問項目への回答の合計値（尺度値）の不偏分散で決まります。合計値の計算の元になる各質問項目間の相関の符号は一致していますか？もし，符号が一致しない項目があると合計を取るときに項目間で多少の値の打ち消しが起こり合計値の分散は小さくなり，必然的に各項目の分散の合計との割合は大きくなります。そして，クロンバックのαは小さくなります。もし，下位尺度にこのような符号の不一致がない場合はその尺度のクロンバックのα信頼性係数は小さくはならないでしょう。そうすれば，全体のクロンバックのα信頼性係数のほうが小さくなってしまうということもありえるでしょう。

以下のようなn=15, 9変数（下位因子数3個ずつ）のデータで
x1  x2  x3  y1  y2  y3  z1  z2  z3
39  37  50  72  50  52  47  63  65
42  48  39  46  47  56  48  48  40
48  54  50  53  43  55  52  55  58
32  30  29  43  26  48  25  50  40
54  39  39  46  57  55  47  31  48
42  50  44  34  38  49  59  29  55
56  62  61  52  50  54  60  49  75
57  51  55  61  56  39  68  51  51
47  65  69  58  59  51  61  61  53
48  50  45  36  51  24  43  43  41
73  66  54  40  56  49  46  60  44
49  48  47  45  43  39  45  54  45
46  50  56  58  67  50  55  49  45
61  56  51  54  59  60  48  57  50
57  44  61  53  47  67  45  51  40
全体のクロンバックのα信頼性係数は 0.8022528 であり，下位尺度のクロンバックのα信頼性係数は 0.8212634, 0.5153157, 0.4307981 となっていることが，エクセルでちょっとした計算をするとわかります。
原因は，相関係数行列をみるとわかりますが，4つの項目間の相関係数がほんのちいさな負の相関を持っているからです。
通常は，尺度を作る際にこのような項目は選ばれないと思いますが，こんなこともあるよということを示すために掲示しておきます。







Re: クロンバックαの係数


　　投稿者：wtcs  2019/05/05(Sun) 17:21  No. 22723   


青木先生

早々のご回答をありがとうございました。
データの因子行列を確認しました。符号の一致しない（負の相関）因子構造はありませんでした。
分散の大きさなどはあまり意識していませんでした。とても勉強になりました。
青木先生のデータ例から、査読者へも回答できそうですし、尺度の信頼性を確認でき安心しました。

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











勉強中








　　投稿者：T-test  2019/05/01(Wed) 15:52  No. 22716   


青木先生

お世話になります。お忙しい中大変申し訳ないのですが、アドバイスを頂ければと思います。

現在、同じ人物が回答したpooled データの解析を行っており、pooled T-test ができないか方法を探しております。

例：

血液型     健康良い群　　　　悪い群
A　　　　　　30　　　　　　　　45　　　
B　　　　　　29　　　　　　　　39　　
O　　　　　　34　　　　　　　　50
AB　　　　　35　　　　　　　　　42

30などの数字は、医療費です。

この場合、各血液毎、例えばA型の列の30と45で、対応のないT-testは、統計としてダメでしょうか？

サンプルは、同じ人物が3か月ごとに回答した結果が3回分あり、その合計です。
Pooled データのため、どのように有意差を解析していいかわかりません。
お手数をお掛けし大変申し訳ございませんが、どうぞよろしくお願い致します。







Re: 勉強中


　　投稿者：青木繁伸  2019/05/02(Thu) 10:57  No. 22717   


実験でサイン的にはあまり芳しいものではないと思うが，繰り返しのない二元配置分散分析でしょうかね？
http://aoki2.si.gunma-u.ac.jp/lecture/TwoWayANOVA/TwoWay.html（繰り返しが1，すなわち繰り返しがない場合）







Re: 勉強中


　　投稿者：T-test  2019/05/02(Thu) 19:17  No. 22718   


青木先生

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

私の書き方が悪かったのですが、同じ回答者が1回から3回、繰り返して回答しています。

データは、
ID   血液　健康状態　医療費
1　　A        良い　　  30
1    A　　　  良い　 　 28　　
1    A　　　　良い　　　29
2    B　　　　悪い　　　35　
2    B　　　　悪い　　　40
2    B　　　　良い　　　32
3    O　　　　良い　　　35
3    O　　　　良い　　　32　

ただし、3回回答している者と2回や1回の者がおります。

このようなデータでも、繰り返しのない二元配置でよろしいでしょうか？

多重比較は、シェッフェでよろしいでしょうか？

分散分析は初めてで、勉強してみたいと思います。
お手数をお掛けし大変申し訳ございませんが、どうぞよろしくお願い致します。







Re: 勉強中


　　投稿者：青木繁伸  2019/05/03(Fri) 14:08  No. 22719   


> このようなデータでも、繰り返しのない二元配置でよろしいでしょうか？

よくわからないけど。無理だと思う。

複数の観察がある場合には（ランダムに）どれか一つに絞れば，対応のない二元配置分散分析にはなりますよね。

分析のことを考えずにデータを取ってしまうと，後で困るという典型例です。
実験デザイン（調査デザイン）をしっかり立ててからデータを取るべきでした。







Re: 勉強中


　　投稿者：T-test  2019/05/03(Fri) 21:25  No. 22720   


青木先生

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

はい、その通りで、データの二次利用のため、非常に困っております。
なんでこういう調査を行うのかと思うのですが、やってしまった後ではどうにもなりません。

対応のない二元配置分散分析、勉強してみます。

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











lda関数の判別係数の再利用について








　　投稿者：初学者  2019/04/14(Sun) 20:16  No. 22713   


青木先生

こんにちは。
いつも掲示板で勉強させていただいております。

このたびは正準判別分析の判別係数を別の研究で再利用できないかと思い、
ご質問させていただきました。
Rのlda関数で分析しております。
print()で出力される判別係数を使って、別の新しいデータの分類に利用できれば、
と考えているのですが、そのようなことは可能なのでしょうか。
R上ではpredict()等で簡単に実現できると思うのですが、
Rが使いこなせない人のために、例えばEXCEL上で判別係数から新データの分類が
簡便にできないかと考えました。

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







Re: lda関数の判別係数の再利用について


　　投稿者：青木繁伸  2019/04/15(Mon) 11:58  No. 22714   


library(MASS)
MASS:::predict.lda
の 2 行を入力すると 89 行のソースプログラムが表示されます。
まるまる全部が必須というわけではないですが，その程度の処理が必要ということです。
それを Excel のワークシート機能だけを使って R を知らない人のために用意するというのはあまりお勧めできません。たとえ用意できたとして，使い方もある程度は複雑になるでしょう。データを用意してスタートボタンを押すだけというようにすることもできるでしょうが，かなりたいへんでしょう。それくらいなら predict(obj, newdata) と入力するのが遙かに簡単では？







Re: lda関数の判別係数の再利用について


　　投稿者：初学者  2019/04/15(Mon) 12:28  No. 22715   


青木先生

早速のご回答、誠にありがとうございました。
たしかに、おっしゃる通りだと感じました。

他の人に使っていただくときには、
・判別分析から学習したobjectをRDSファイル等にセーブ
・RDSファイル、実行するRのscript、簡単な解説、などを渡す
・predict(obj,newdata)してもらう
等の方法で検討してみたいと思いました。

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











マンホイットニーのＵ検定〜同順位の補正〜








　　投稿者：確率は難しい  2019/04/03(Wed) 14:29  No. 22708   


青木先生　どうかよろしくお願いします。

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

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

しかし、満月の日の翌日は平均気温が低くなるということで考えると、小数点以下３桁ぐらいまで平均気温を計算して同順位がないとすると、
第１群 気温９.２５３、１０．８６５度・・・・など、９度〜１１度の日が１５個
第２群 気温８．０２５度のデータが１個（最も低い）
これで計算してみると、確率は0.1ということになりました。

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

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

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

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







Re: マンホイットニーのＵ検定〜同順位の補正〜


　　投稿者：青木繁伸  2019/04/03(Wed) 18:03  No. 22709   


stats の wilcox.test で，同順位があるために「タイがあるため、正確な p 値を計算することができません」となる場合は，coin パッケージの wilcox_test を使ってみましょう。
library(coin)
x <- c(rep(0, 15), 100)
g <- factor(c(rep("a", 15), "b"))
wilcox_test(x ~ g, distribution="exact")
以下が表示されます。
	Exact Wilcoxon-Mann-Whitney Test

data:  x by g (a, b)
Z = -3.873, p-value = 0.0625
alternative hypothesis: true mu is not equal to 0
Z は，同順位補正をしたときの正規近似の Z 値で，表示されている p 値の算出には使われていません。

同時に表示されている p-value = 0.0625 は正確な p 値です。

coin パッケージの出現前には exactRankTests パッケージの wilcox.exact が使われていました。
library(exactRankTests)
wilcox.exact(x ~ g)
以下が表示されます。
	Exact Wilcoxon rank sum test

data:  x by g
W = 0, p-value = 0.0625
alternative hypothesis: true mu is not equal to 0
W は wilcoxon 統計量です（どちらを1群にするかで，もう一方の極端な側の W = 15 が返されることもあります）。表示された p 値は，coin:::wilcox_test と同じ 0.0625 です。

もう少し他の道具でも試してみるには，
http://aoki2.si.gunma-u.ac.jp/R/index.html
の中の，
XV. exact 検定とモンテカルロ法による近似検定
	3. マン・ホイットニーの U 検定(exact test) 
に記載した R プログラム exact.mw を使ってみると以下のようになります。
exact.mw(x=rep(0, 15), y=1)
U = 0, Z = 3.87298, P 値 = 0.000107511
正確な P 値 = 0.0625
査察した分割表の個数は 2 個
示されている U = 0, Z = 3.87298 は，coin:::wilcox_test の Z と同じですね。

次の行に正確な p 値が示されています。0.0625 で，coin:::wilcox_test，exactRankTests:::wilcox.exact の p 値と同じです。

exact.mw は exact.mw(matrix(c(15, 0, 0, 1), 2))  のようにも書けますが，
http://aoki2.si.gunma-u.ac.jp/exact/utest/getpar.html
にオンライン版を用意してあります。
列数を 2 として，次の画面で
15 0
0 1
の2×2分割表を入力すれば， 0.0625 が表示されます。

ということで，正確な p 値は 0.0625 でよいようです。

====================

最後に，「乱数をプラスして同順位をなくして解決する方が、同点の場合の数式を使うよりも実態に合っているように思う」ということですが，乱数を加えて検定を行うという流儀は存在します。
しかし，今回の場合は
stats:::wilcox.test(1:15, 16)

	Wilcoxon rank sum test

data:  1:15 and 16
W = 0, p-value = 0.125
alternative hypothesis: true location shift is not equal to 0
とすることと同じなので，わざわざ乱数を加えて同順位をなくすということがどの程度適切か疑問ですね。







Re: マンホイットニーのＵ検定〜同順位の補正〜


　　投稿者：確率は難しい  2019/04/04(Thu) 10:16  No. 22710   


青木先生

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

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

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







Re: マンホイットニーのＵ検定〜同順位の補正〜


　　投稿者：青木繁伸  2019/04/04(Thu) 12:19  No. 22711   


それぞれのデータを wilcox_test(x ~ g, distribution="exact") で分析する数行のプログラムを書けばよいだけでしょう。

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







Re: マンホイットニーのＵ検定〜同順位の補正〜


　　投稿者：確率は難しい  2019/04/04(Thu) 15:28  No. 22712   


先生　ありがとうございます。

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

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











二項ロジスティック回帰分析








　　投稿者：春  2019/03/30(Sat) 02:00  No. 22701   


青木先生

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

仮説に基づいて、独立変数群を5モデル作り、それぞれの影響をみるため順番に強制投入しています。以下の２つの方法を試みましたが、投入したモデルは同じなのに結果が違っています（B,標準偏差、wald,有意差、EXP（B)）。なぜ違っているのか、また、どちらのやり方が適しているのかを教えていただけないでしょうか。

�‘販�変数投入の際にブロックごとに次々と指定していった場合
�▲屮蹈奪�ごとの指定せずに「総当たり法」を繰り返し、�，里箸�と同じブロックごとに変数を増やした場合

以上です。

お忙しいところ、申し訳ありませんが、どうか、よろしくお願いいたします。







Re: 二項ロジスティック回帰分析


　　投稿者：青木繁伸  2019/03/30(Sat) 09:16  No. 22702   


（2）の実際の分析過程があいまいなのではっきりしたことが言えませんが，一般的には，「分析に用いる変数セットにより，結果的に選択される変数セットが異なることはよくあること」

（1）の「独立変数投入の際にブロックごとに次々と指定」する場合は，投入順序により結果が異なることがある。増減法であろうと同じことが起きる。
増加法と減少法，増減法と減増法で結果が違うのもよくあること。

「ブロックごとの指定せずに「総当たり法」を繰り返し、（1）のときと同じブロックごとに変数を増やした場合」も同じ。ブロックごとに総当たり法で増やすべき変数を決めても，結局は（1）と同じように，投入順序で結果が違うこともある。

「分析に用いる変数セット」に依存しない方法は，全ての変数を対象に総当たり法でモデルを選ぶ。もっとも，理論的に解釈可能なモデルであるかはチェックする必要がある。

=======

以下のような思考実験をするとよく分かる。

以上のようにして「最良のモデル」が得られた後，ある独立変数を分析に加えたらどうなるだろうか？

その変数が橋にも棒にもかからないヘボ変数で，結果は変わらない。
その変数が，新たにモデルに加えられる。
その変数が，モデルに加えられる一方，今までにモデルに含まれていた変数が不適切とわかり除去される（変数1個とは限らない）
その変数が，モデルに加えられる一方，今までにモデルに含まれていた変数が不適切とわかり除去される（1個とは限らない），そして今までモデルに含まれていなかった変数がモデルに加えられる（1個とは限らない）。
その際，同じブロック内の変数に入れ替えが生じる。
また，ブロック内に適当な変数がなくなる。
さらに，ブロック内の複数の変数がモデルに含まれる。
などなど，全ての状況が考えられる。

だって，それが多変量解析だもの







Re: 二項ロジスティック回帰分析


　　投稿者：春  2019/03/30(Sat) 12:09  No. 22703   


青木先生

早々に明快なご説明をいただきまして、誠に、ありがとうございます。

�，皚△睚竸瑤療蠧�順序は同じで、結果的にも変数は除去されなかったのですが、値（B,標準偏差、wald,有意差、EXP（B)）が異なるのは、なぜだろうかと不思議に思っていました。
しかし、先生が明示してくださった「全ての変数を対象に総当たり法でモデルを選ぶ」方法で解析を進めようと思います。ありがとうございました。

今後も、どうぞよろしくお願い致します。







Re: 二項ロジスティック回帰分析


　　投稿者：青木繁伸  2019/03/30(Sat) 16:58  No. 22704   


文意がとりにくかったのですが，「最終的に採用された変数は同じ」ということですか？

データに欠損値があったりしますか？
2つの方法で，実際に使用されたサンプルサイズがいくつかわかりますか？
変数が同じでも，実際に分析に使われるサンプルサイズが異なれば，B,標準偏差、wald,有意差、EXP（B)）が異なるのは当たり前です。







Re: 二項ロジスティック回帰分析


　　投稿者：春  2019/03/31(Sun) 23:27  No. 22705   


青木先生

ご連絡ありがとうございます。
本当に、このように親身になって頂いてありがたいです。

ご質問して頂いてように「最終的に採用された変数は同じ」です。
」
そして、「実際に使用されたサンプルサイズ」」を確認すべく、もう一度、spssの結果のlog�，鉢△鮓�比べてみました。すると、青木先生の仰るようにサンプルサイズが違っていました！
�，諒�法だと投入した全ての変数の欠損値を含むので、�△離汽鵐廛襯汽ぅ困茲蠅眈�さくなっていました。
これが原因だったのですね。よくわかりました。

ちなみに、�，任呂匹離皀妊襪�Nの値が同じになりますが、�△世肇皀妊襪瓦箸�Nの値が違ってくるということになりますね。論文にモデル1、モデル2、モデル3、モデル4、モデル5と階層的に示したいときには、Nを揃えた�，諒�法が適切なのでしょうか？
何度も質問させていただき、恐縮ですが、ご教示いただけないでしょうか。

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







Re: 二項ロジスティック回帰分析


　　投稿者：青木繁伸  2019/04/01(Mon) 08:29  No. 22706   


モデルの評価を AIC などで行う場合は，サンプルサイズが同じでないと評価不能ですね。

ただし，対象者の属性により欠損値が生じやすくなるような状況がある場合は，結果の評価に注意が必要になるかもしれません。

なお，変数選択の結果を見て，モデルに残った変数だけを指定して再度分析を行うべきです。
以下のような状況を考察すれば分かるでしょう。
変数が a ~ e，対象者が 9 人。o は測定値がある（欠損値ではない），x は欠損値。
       a b c d e
case1  o o o o x
case2  o o o x x
case3  o o x o o
case4  o o o o o
case5  o o o o o
case6  o o o o o
case7  o o o o o
case8  o o o o o
case9  o o o o o
5変数全部を使って分析を始める場合，c, d, e が欠損値の case1, case2, case3 は分析対象から外されます。分析対象は 6 人です。
その後の変数選択で  a, b, d が残ったとします（c, e は採用されなかった）。
これで終わってはいけません。
a, b, d を使って分析することにすると，欠損値を持つとして分析から外されるのは case2 のみとなります。分析対象は 8 人です。 case1 も case3 も分析対象となり，2 人増えます。







Re: 二項ロジスティック回帰分析


　　投稿者：春  2019/04/01(Mon) 10:56  No. 22707   


青木先生

何度も丁寧にご教示いただきまして誠にありがとうございました。

ようやく分かってまいりました。
Nを揃えるときの注意点に気をつけながら、モデルに残った変数でやり直しをしてみます。

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











Split plot + RCBD with repeated measures








　　投稿者：森  2019/03/25(Mon) 18:32  No. 22700   


青木先生

自分で調べてみましたが、解決できなかったので質問させていただきます。
本当に申し訳ありませんが、よろしくお願いいたします。

10つの森林を対象に、それぞれの森林に種Aの樹木1個体と樹種Bの樹木1個体を植栽し、それらの個体の高さ（樹高）を2ヶ月おきに観察しています。
私が明らかにしたいのは、樹高に対する種の影響と経時的な変化（測定時期の影響）、それらの交互作用です。
森林をブロックとして捉え、その中に、処理区（樹種）を設置し、さらに、その経時変化を調べていますが、それぞれのブロックの処理区に個体が1つずつしか存在しません。
経時変化の測定はそれぞれの個体を繰り返し測定します（repeated measures）。
これを統計ソフトRで解析したいのですが、一般化線形混合モデルを用いる場合、変量効果をどのようにモデルに組み込めばいいのかが分かりません。
glmer(樹高~樹種＋(1|森林/樹種/測定時期)＋測定時期＋樹種×測定時期)
などの形を考えていますが、これでいいのかも分かりません。
よろしければご教示いただけないでしょうか。

また、2つ目の質問になってしまいますが（不躾で本当に申し訳ありません）、上記の実験デザインの中で、それぞれのブロックの処理区に2つの個体を配置した場合はどのようなモデルになりますでしょうか。
こちらにつきましても、
glmer(樹高~樹種＋(1|森林/樹種/個体/測定時期)＋測定時期＋樹種×測定時期)
などの形を考えていますが、よくわかっていません。
よろしければお教えいただけないでしょうか。

私の勉強不足なところをお尋ねしてしまい、大変申し訳ありません。
お手数をおかけいたしますが、どうぞよろしくお願いいたします。











最小選択肢にほとんど丸のつかない可能性のある尺度開発








　　投稿者：風  2019/03/11(Mon) 08:13  No. 22695   


青木先生
初めて質問させていただきます。
初学者ですよろしくお願い致します。

現在、自己評価尺度の開発を行っています。尺度原案の信頼性妥当性検証を行っている
ところです。
尺度原案のリッカートを6件法「1．全く当てはまらない」、「2．当てはまらない」、「3．どちらかというと当てはまらない」、「4．どちらかというと当てはまる」、「5．当てはまる」、「6．非常に当てはまる」として調査を行いました。
返送されてきた調査票は、「4.どちらかというと当てはまる」「5．当てはまる」に回答が集まる傾向がありました。

このように最小選択肢にほとんど〇のつかない可能性のある尺度はどうなのでしょうか…？

（リッカートを変更する予定ですが…）







Re: 最小選択肢にほとんど丸のつかない可能性のある尺度開発


　　投稿者：青木繁伸  2019/03/11(Mon) 09:20  No. 22696   


もっと極端な場合，つまり一つの選択肢にほとんど全てが回答するという場合は，そのような質問は不適切ということになりますね。
そのような場合と比べてどうすべきかということでしょう。明確な基準というものはないと思います。他の研究者とのコンセンサスが得られるかどうかも判断基準になるでしょう。







Re: 最小選択肢にほとんど丸のつかない可能性のある尺度開発


　　投稿者：風  2019/03/11(Mon) 12:20  No. 22697   


青木先生

早速お返事をいただきありがとうございました。
極端な場合は、質問項目として除外されることになると考えます。

他の研究者の意見も伺って検討していきたいと思います。ありがとうございました。











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





- J o y f u l 　N o t e -

Modified by　i s s o