このドキュメントは、続・『続・心理統計学の基礎』（南風原朝和 著）読書会について作成した第3章3節以降のメモという位置づけです。主にRで簡単に実行できるコードを中心に構成しています。ただこのドキュメントのみでは意味が通らないかと思いますので、続・心理統計学の基礎と共にご覧ください。
共通の分散を持つ2つの正規母集団\(N(\mu_1, \sigma^2),N(\mu_2, \sigma^2)\)から、それぞれ大きさ\(n_1,n2\)の標本を独立にとるとき、2群間の平均値差の検定のための検定統計量は以下により算出されます: \[ \begin{eqnarray} t=d\times\sqrt{\frac{n_1 n_2}{n_1 + n_2}} \end{eqnarray} \]
なお、この\(d\)は標本標準化平均値差で、以下の式から求めます: \[ \begin{eqnarray} d=\frac{\bar{y_1}-\bar{y_2}}{s^*} \end{eqnarray} \]
この\(\bar{y_1},\bar{y_2}\)は各群の標本から算出した平均値なので、分子は標本2群の平均値間の差となります。分母は母集団標準偏差の推定量で、以下のような式で求めます:
\[ \begin{eqnarray} s^* = \sqrt{\frac{n_1 s_1^2 + n_2 s_2^2}{n_1 + n_2 - 2}} \end{eqnarray} \]
この\(s^*\)は要するに、標本2群の標準偏差\(s_1,s_2\)をプールしています。ようするに、この統計量は標本における2群間の平均値差が、(母集団)標準偏差の何倍くらいかを示すのが\(d\)です。効果量っぽいですね。でもこの算出された\(t\)値には2つの前提があります:
まず後者について、これはよく対応のないt検定で出てくるお話ですのでここでは置いておきます。前者の帰無仮説についてですが、この\(t\)は自由度\(n_1 + n_2 -2\)の\(t\)分布に従います。そうでなければ前章でお話したように分布が\(t\)分布から歪んでしまいます。そう、非心\(t\)分布に従います。なお非心度は以下のように算出されます:
\[ \begin{eqnarray} \lambda = \delta \times \sqrt{\frac{n_1 n_2}{n_1 + n_2}} \end{eqnarray} \]
なお、この\(\delta\)は母集団標準化平均値差となり、以下の式となります:
\[ \begin{eqnarray} \delta = \frac{\mu_1 - \mu_2}{\sigma} \end{eqnarray} \]
この\(\mu_1,\mu2\)は母集団における各群の平均値です。帰無仮説\(H_0 : \mu_1 = \mu_2\)が偽なので、この分子は0でない何らかの値をとります。\(\sigma\)は母集団における標準偏差(添字が入らないのは2群で共通だから)です。これを\(t\)に組み込んで指標を算出するお話は前章で触れているので省略します。
さて、この\(\delta\)は「母集団における2群の(標準化された)ズレ具合(を推定した値)」です。つまり母集団における効果量ということもできます。ではこの値はどうやって求めるのかですが、\(d\)と共通(不偏推定量)となります。なので効果量報告では\(d\)が報告されます。
この非心度\(\lambda\)と自由度\(n_1 + n_2 -2\)から、検定力を算出できますね。サンプルサイズ、決めれます。なおテキストの注16にあるとおり、この\(d\)はHedges`\(g\)とよばれることがあるようで、さらに色々面倒なことがあるようです。テキストでは「標準化平均値差\(d=***\)」という表現が望ましいとあります。
基本は本文を読んでください。以下ポイントを箇条書きで。
delta <- 1 #母集団効果量を設定
U3 <- pnorm(delta,mean=0,sd=1)
U3## [1] 0.8413447delta <- 1 #母集団効果量を設定
PoD <- pnorm(delta/sqrt(2),mean=0,sd=1)
PoD## [1] 0.7602499テキストの例をRで算出してみます。なおMBESSパッケージを使用するところは修正しています:
# 各種パラメータの設定
# dは平均値差と各群の標準偏差とサンプルサイズにより計算できます
alphaaa <- 0.05 #α
n1 <- 20
n2 <- 20
d <- 0.8 #標本より算出された母集団効果量の推定値
# 検定統計量tの算出
t_value <- d*sqrt((n1*n2)/(n1+n2))
t_value## [1] 2.529822# 非心度の下限値・上限値の算出
# テキストで紹介されているMBESSパッケージを利用します
# install.package("MBESS") # 未インストール時は左のコードを実行
library("MBESS")
lambdaaaa <- conf.limits.nct(t.value = t_value, df=(n1+n2-2),conf.level=(1-alphaaa)) #テキストからconf.levelに変更
lambdaaaa## $Lower.Limit
## [1] 0.4741762
##
## $Prob.Less.Lower
## [1] 0.025
##
## $Upper.Limit
## [1] 4.55476
##
## $Prob.Greater.Upper
## [1] 0.025# 非心度の下限値・上限値を、後述の式が見やすいよう変数に格納
lambda_L <- lambdaaaa$Lower.Limit
lambda_U <- lambdaaaa$Upper.Limit
c(lambda_L,lambda_U)## [1] 0.4741762 4.5547602# 効果量の信頼区間の限界値を算出
delta_L <- lambda_L*sqrt((n1+n2)/(n1*n2))
delta_U <- lambda_U*sqrt((n1+n2)/(n1*n2))
c(delta_L,delta_U)## [1] 0.1499477 1.4403416基本省略しますが、ブートストラップ法が有効では? という話です。
詳細は省略します。以下気をつけるポイントです:
詳細は省略します。以下気をつけるポイントです:
省略します。
省略します。なお計算は対応がない場合よりも簡単となります。
2つのカテゴリ変数間の連関を検定する場合、\(\chi^2\)を算出します:
\[ \begin{eqnarray} \chi^2 = \sum_{i=1}^a\sum_{j=1}^b\frac{(n_{ij}-e_{ij})^2}{e_{ij}} \end{eqnarray} \]
これは\(\chi^2\)検定でよく見る形なので説明は省略します。ここで算出された\(\chi^2\)値は、帰無仮説が真であれば自由度が\((a-1)(b-1)\)の\(\chi^2\)分布に従います。でも帰無仮説が偽であれば非心\(\chi^2\)分布に従い、その非心度は以下のようになります:
\[ \begin{eqnarray} \lambda = \sum_{i=1}^a\sum_{j=1}^b\frac{(\pi_{ij}-\pi_{0ij})^2}{\pi_{0ij}} \times N \end{eqnarray} \]
よって、母集団効果量は非心度からNの要素を抜いた、以下のようになります:
\[ \begin{eqnarray} \gamma = \sum_{i=1}^a\sum_{j=1}^b\frac{(\pi_{ij}-\pi_{0ij})^2}{\pi_{0ij}} \end{eqnarray} \]
母集団効果量\(\gamma\)に対応する標本統計量は\(\chi^2/N\)ですが、行数・列数に依存します(要するに影響をうけてしまう)。なので効果の大きさを見たいと思っても、行列サイズでその意味が変わるのでは困ります。そこでテキストではクラメルの連関係数\(V\)を紹介しています:
\[ \begin{eqnarray} V = \sqrt{\frac{\chi^2}{(min(a,b)-1)N}} \end{eqnarray} \]
省略です。テキストのままです。
テキストにある数値例を実際にRで計算してみます
# パラメータの設定
n <- 50
a <- 2 #行数
b <- 3 #列数
chi <- 7.8 #かいにじょーち
alphaaa <- .05
# クラメルの連関係数の算出
V <- sqrt(chi/((min(a,b)-1)*n))
V## [1] 0.3949684# このchi^2.valueのときの、非心度らむだぁの95%信頼区間の下限値・上限値を算出
# MBESSパッケージを利用します
library("MBESS")
lambdaaaa <- conf.limits.nc.chisq(Chi.Square=chi,df=(a-1)*(b-1),conf.level=1-alphaaa)
lambdaaaa## $Lower.Limit
## [1] 0.1173012
##
## $Prob.Less.Lower
## [1] 0.025
##
## $Upper.Limit
## [1] 21.27973
##
## $Prob.Greater.Upper
## [1] 0.025# イメージしやすいように値を変数に格納しときます
lambda_L <- lambdaaaa$Lower.Limit
lambda_U <- lambdaaaa$Upper.Limit
c(lambda_L,lambda_U)## [1] 0.1173012 21.2797331# 上記のラムダを使って効果量をそのまま算出もできますが、テキストにあわせてクラメルの連環計数の95%信頼区間の上限値・下限値を算出します
gamma_L <- sqrt(((a-1)*(b-1)+lambda_L)/((min(a,b)-1)*n))
gamma_U <- sqrt(((a-1)*(b-1)+lambda_U)/((min(a,b)-1)*n))
c(gamma_L,gamma_U)## [1] 0.2057815 0.6823450はい、一致しましたね。ちなみに注27にあるとおり、\(\chi^2\)分布は非負ですので、場合によってはエラーを出すことがあります:
# パラメータの設定
n <- 50
a <- 2 #行数
b <- 3 #列数
chi_dame <- 0.23 #適当にいれてみた
alphaaa <- .05
# not run(走りません).のでコメントアウトしてます
#lambdameee <- conf.limits.nc.chisq(Chi.Square=chi_dame,df=(a-1)*(b-1),conf.level=1-alphaaa)
#lambdameeeこのような場合、conf.level=NULLと水準として設定せず、上側確率と下側確率をそれぞれ個別に設定します。なお下側は(どうせ)マイナスなので0にしときます:
lambdaiyaaaa <- conf.limits.nc.chisq(Chi.Square=chi_dame,df=(a-1)*(b-1),conf.level=NULL,alpha.lower=0,alpha.upper=alphaaa/2)
lambdaiyaaaa## $Conf.Interval.type
## [1] "one-sided"
##
## $Lower.Limit
## [1] 0
##
## $Upper.Limit
## [1] 3.111237
##
## $Prob.Greater.Upper
## [1] 0.025以上です。