生物統計学の質問

以下の問題を解いてほしいです。協力お願いします。
dataA = c(18.3, 16.9, 18.1, 17.7, 17.9, 17.3, 18.7, 18.3, 16.3, 15.0, 17.9, 18.2, 15.3, 18.5, 16.9, 18.4, 17.4, 17.9, 17.8, 16.3, 17.4, 17.2, 16.7, 17.4, 17.7, 15.5, 17.3, 16.3, 15.1, 16.4)
dataB = c(17.8, 19.1, 18.6, 17.6, 18.7, 19.2, 15.9, 18.9, 19.1, 18.7, 13.1, 16.2, 16.8, 18.5, 18.0, 17.5, 18.4, 16.1, 18.1, 16.3, 16.5, 17.2, 18.2, 15.9, 16.0, 17.7, 18.7, 18.8, 18.2, 18.6)

① それぞれのデータにおける「平均値（算術平均）」，「自由度」，「不偏分散」を求めよ．

平均値
mean(dataA)
mean(dataB)

自由度
length(dataA)-1
length(dataB)-1

計算式で出した不偏分散
sum((dataA -mean(dataA))^2)/(length(dataA)-1)
sum((dataB -mean(dataB))^2)/(length(dataB)-1)

② それぞれのデータにおける「標準誤差」および「95%信頼区間」を求めよ．

標準誤差 標準偏差 / sqrt(n)

sqrt(var(dataA)) / sqrt(length(dataA))
sqrt(var(dataB)) / sqrt(length(dataB))

95％信頼区間

mean(dataA)+(qnorm(0.95) * sd(dataA))
mean(dataA)-(qnorm(0.95) * sd(dataA))

mean(dataB)+(qnorm(0.95) * sd(dataB))
mean(dataB)-(qnorm(0.95) * sd(dataB))

③ それぞれのデータを平均0，標準偏差1のデータにZ変換せよ．

④ 5％水準において，dataAとdataBの平均値に差があるかどうか判定せよ．なお，回答には対立仮説および帰無仮説を明記すること．また「入力したコマンド」や「出力結果」を示すとともに，「判定するために読み取った場所」および「なぜそのように判定したのか」を説明すること．

t.test(x=dataA,y=dataB,conf.level=0.95,var.equal=F)

⑤ 解析の結果から，品種AおよびBの糖度にはどのようなことを言うことができると考えられるか，考察せよ

No.7 回答者： kamiyasiro 回答日時： 2022/05/18 08:57

#6です。

ご質問者の書かれたＲのスクリプトは、先生が示されたものですか。う～ん、ですね。

ちゃんとした教科書は、検定・推定においては、『母分散既知の場合』と『母分散未知の場合』で区別して、正規分布とｔ分布を使い分けています。

ご質問者も、ちゃんとした教科書で学んだ方が良いと思います。

No.6 回答者： kamiyasiro 回答日時： 2022/05/18 08:28

訂正というか、他の方のご意見も伺いたいということで・・・。

データの95％信頼区間ですが、初学者はｔ分布 qt() でなく正規分布 qnorm() を使えば良いかも。
ただ、数値0.95は0.975に修正して下さい。

mean(dataA)+(qnorm(0.975) * sd(dataA))
mean(dataA)-(qnorm(0.975) * sd(dataA))

mean(dataB)+(qnorm(0.975) * sd(dataB))
mean(dataB)-(qnorm(0.975) * sd(dataB))

ただ、矛盾はあります。なぜなら、sd() は不偏分散の平方根で、そもそもサンプリングデータを前提としているんだ、と言っているようなもんです。なのに、想定している分布は全数データを前提としているものを使うんですから。

そもそも、各30個のデータがサンプリングデータだとは、どこにも書いてないですが、問題には「『不偏分散』を求めよ」とあります。だから、これはサンプリングデータかもしれないのです。

で、追い打ちを掛けるようで申し訳ないですが、もしこれがサンプリングデータでなく全数データなら、次のようにしなければなりません。ただし、n はデータ数とします。

mean(dataA)+(qnorm(0.975) * sqrt(var(dataA) * (n-1) / n))
mean(dataA)-(qnorm(0.975) * sqrt(var(dataA) * (n-1) / n))

mean(dataB)+(qnorm(0.975) * sqrt(var(dataB) * (n-1) / n))
mean(dataB)-(qnorm(0.975) * sqrt(var(dataB) * (n-1) / n))

なにをやっているかと言うと、不偏分散 var() は(n-1) で割った値なので、一旦それを掛けてｎで割り直します。これによって標本分散になります。標準偏差はその平方根を使います。
こうすれば、qnorm() を使うことと矛盾はしません。

いったい、正解はどっちだと言うんでしょうか。

なお、ウチの会社では、信頼区間、予測区間はｔ分布で求めるよう教育しています。市販の統計ソフトもそうなっていますしね。

No.5 回答者： kamiyasiro 回答日時： 2022/05/17 23:20

#2,3です。

Ｒのｔ検定はデフォルトで両側検定を行いますが、念のため、以下のように両側検定であることを明示しておいた方が良いです。

と言うのも、var.equal = F を入れなくても、Rはデフォルトで「ウェルチの検定」を行います。これを書くくらいなら、両側検定であることも明示した方が良いと思った次第です。


t.test(x = dataA, y = dataB, alternative = c("two.sided"), conf.level = 0.95, var.equal = F)

No.4 回答者： kamiyasiro 回答日時： 2022/05/17 22:57

#2,3です。

前の投稿のヒストグラムを見て頂くと、dataBには外れ値があることが分かります。
このような外れ値は測定ミスや異常品と考えられるので、サンプルから除外することが必要です。

すると検定結果は、「両者の差は、5%有意水準で有意である」と変化します。この点を考察に加えて頂ければと思います。



> # 外れ値を除外して検定すると
> dataB <- dataB[dataB > 14]
> t.test(x = dataA, y = dataB, conf.level = 0.95, var.equal = F)

Welch Two Sample t-test

data: dataA and dataB
t = -2.0287, df = 56.531, p-value = 0.04721
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.124049963 -0.007214405
sample estimates:
mean of x mean of y
17.20333 17.76897

No.3 回答者： kamiyasiro 回答日時： 2022/05/17 22:48

このサイトは空白を詰めて表示するので、非常に見辛いですが、結果は以下のとおりです。

検定結果ですが、最後の出力の
p-value = 0.1985
を見れば良く、結論は「品種AおよびBの糖度には差があるとは言えない」です。

> # Rで実行
> rm(list = ls()) # ご破算で願いましては
>
> dataA <- c(18.3, 16.9, 18.1, 17.7, 17.9, 17.3, 18.7, 18.3, 16.3, 15.0, 17.9, 18.2, 15.3, 18.5, 16.9, 18.4, 17.4, 17.9, 17.8, 16.3, 17.4, 17.2, 16.7, 17.4, 17.7, 15.5, 17.3, 16.3, 15.1, 16.4)
> dataB <- c(17.8, 19.1, 18.6, 17.6, 18.7, 19.2, 15.9, 18.9, 19.1, 18.7, 13.1, 16.2, 16.8, 18.5, 18.0, 17.5, 18.4, 16.1, 18.1, 16.3, 16.5, 17.2, 18.2, 15.9, 16.0, 17.7, 18.7, 18.8, 18.2, 18.6)
>
> par(mfrow = c(2, 1))
> bins = seq(12, 20, length = 21)
>
> hist(dataA, xlim = c(12, 20), breaks = bins)
> abline(v = mean(dataA), col = 2)
> hist(dataB, xlim = c(12, 20), breaks = bins)
> abline(v = mean(dataB), col = 2)
>
> # 平均値
> mean(dataA)
[1] 17.20333
> mean(dataB)
[1] 17.61333
>
> # 自由度
> length(dataA) - 1
[1] 29
> length(dataB) - 1
[1] 29
>
> # 計算式で出した不偏分散
> sum((dataA - mean(dataA))^2) / (length(dataA) - 1)
[1] 1.080333
> var(dataA)
[1] 1.080333
>
> sum((dataB - mean(dataB))^2) / (length(dataB) - 1)
[1] 1.894989
> var(dataB)
[1] 1.894989
>
> # 標準誤差 標準偏差 / sqrt(n)
> sqrt(var(dataA)) / sqrt(length(dataA))
[1] 0.1897659
> sqrt(var(dataB)) / sqrt(length(dataB))
[1] 0.251329
>
> # データの95％信頼区間（ここは修正してあります）
> df <- length(dataA)-1
> mean(dataA)+(qt(0.975, df) * sd(dataA))
[1] 19.32913
> mean(dataA)-(qt(0.975, df) * sd(dataA))
[1] 15.07754
>
> df <- length(dataB)-1
> mean(dataB)+(qt(0.975, df) * sd(dataB))
[1] 20.42877
> mean(dataB)-(qt(0.975, df) * sd(dataB))
[1] 14.7979
>
> # 平均値の95%信頼区間（ここは追加しました）
> fit <- lm(dataA ~ 1)
> confint(fit, level = 0.95)
2.5 % 97.5 %
(Intercept) 16.81522 17.59145
>
> fit <- lm(dataB ~ 1)
> confint(fit, level = 0.95)
2.5 % 97.5 %
(Intercept) 17.09931 18.12736
>
> # それぞれのデータを平均0，標準偏差1のデータにZ変換せよ．
> # scale(dataA) 結果は割愛
> # scale(dataB) 結果は割愛
>
> # 平均値の差の検定
> t.test(x = dataA, y = dataB, conf.level = 0.95, var.equal = F)

Welch Two Sample t-test

data: dataA and dataB
t = -1.3019, df = 53.955, p-value = 0.1985
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.0413975 0.2213975
sample estimates:
mean of x mean of y
17.20333 17.61333

>

No.2 回答者： kamiyasiro 回答日時： 2022/05/17 22:41

このまま、Rに突っ込んでやるだけ、と思いますよ。

でも、間違いが多々ありますね。95%信頼区間を求めるときは、0.975として片側2.5%ずつにしないといけません。

それに、データはサンプリングデータだと思いますが、そうであるなら信頼区間はなぜ正規分布で計算するのですか。ｔ分布を使うべきです。でないと、最後のｔ検定と矛盾が生じますよ。

それらを直すと、

データの95％信頼区間
df <- length(dataA)-1
mean(dataA)+(qt(0.975, df) * sd(dataA))
mean(dataA)-(qt(0.975, df) * sd(dataA))

df <- length(dataB)-1
mean(dataB)+(qt(0.975, df) * sd(dataB))
mean(dataB)-(qt(0.975, df) * sd(dataB))

③ それぞれのデータを平均0，標準偏差1のデータにZ変換せよ．
scale(dataA)
scale(dataB)

④検定せよ．
帰無仮説　H0：μA = μB
対立仮説　H0：μA ≠ μB

次の投稿では、解析結果をアップします。

No.1 回答者： yhr2 回答日時： 2022/05/17 14:07

なんか、「統計学」の「序論」である「記述統計」の基本式の復習のための演習問題ですよね。

「統計学」の本論は、それを前提にした「推測統計」ですから、そのための準備としてしっかりマスターしておいてください。

単に「そこに書かれている定義式どおりにガシガシと数値計算」するだけの話です。
というより、パソコンか何かの統計ツールを使いこなすための演習なんじゃないの？「データを入力する」のが大変なだけ。