Bing Image Creator で生成: Close up of a pink Rhododendron flower, background is blue clear sky, photo
www.crosshyou.info
の続きです。前回はおおまかな分析で2009年と2022年の体力テストの結果を比較して、2009年のほうがテスト結果が良いようだとわかりました。
今回は、もう少しちゃんと分析してみます。
前回は、学年の違いや男女の別を考慮せずに、2009年の平均値と2022年の平均値を計算しました。
今回は、学年の違いと男女の別考慮して分析してみます。
infer という統計的推論のパッケージを使用します。
Tidy Statistical Inference • infer <<< このサイトに詳しい使い方があります。
まず、infer パッケージを library() 関数で読み込みましょう。
year は 2009 と 2022 の二つの値しかないので、ファクター型に変換しておきます。
そうしたら、spceify() 関数と fit() 関数を使い、observed fit を求めます。
year2022 の係数は -0.0154 です。つまり、学年や男女差を考慮しても、2022年のほうが2009年よりも体力テストの結果は 0.0154 低いということです。
これが、統計的に有意に 0 と違うのかどうかが問題です。
ブートストラップ法で、いまの作業を 1000 回繰り返します。
year の係数だけ抜き出してみましょう。
1 回目の推計では、係数がプラスです。2 回目から 10回目はマイナスで推計されました。
95%信頼区間を求めます。
year2022 の 95% 信頼区間は -0.0529 ~~ 0.0218 です。0 を含んでいますから、学年と男女別を考慮した後では、2009 年と 2022 年では統計的に有意な違いは認められない、ということですね。
ということで、学年、男女の別を考慮した状況で 2009 年と 2022 年の体力測定の結果は、平均すると、2022 年のほうが悪いけど、統計的に有意な違いとは言えないということがわかりました。
もちろん、これは、1 学年、男子、女子 一人だけの結果とみなしたものです。
実際は、多数の児童・生徒の平均値がこのデータになっているので、2022年のほうが体力テストの結果は悪くなっているのでしょうね。
なので、仮想の多数の児童・生徒のデータセットを作って同じことをやってみます。
まず、仮想のデータセットを作ります。
生徒の数は 1,0000 人にしてみました。東京都の生徒・児童数は小学生は 60 万人、中学生は 20万人 ぐらいだそうですが、データ数を多いと計算時間も長くなりそうなので、とりあえず、このくらいでやってみます。
あとは、同じ作業を繰り返すだけです。
observed fit を求めて、
ブートストラップ法で繰り返して、
信頼区間を求めます。
year2022 の係数の 95% 信頼区間は、-0.0256 ~~ -0.00739 と 0 を含んでいません。つまり、2022年 は有意に体力テストの結果は 2009年と比較して悪い、ということです。
今回は以上です。
初めから読むには、
です。
今回のコードは以下になります。
#
# infer パッケージを読み込む
library(infer)
#
# year をファクター型の変数に変換しておく
df <- df |>
mutate(year = factor(year))
df
#
# データから観測されたモデル
obs_fit <- df |>
specify(scaled ~ year + grade + gender) |>
fit()
obs_fit
#
# ブートストラップ法
set.seed(113)
boot_dist <- df |>
specify(scaled ~ year + grade + gender) |>
generate(reps = 1000, type = "bootstrap") |>
fit()
#
# year の係数だけ表示
boot_dist |>
filter(term == "year2022")
#
# 95%信頼区間
conf_ints <- get_confidence_interval(
boot_dist,
level = 0.95,
type = "percentile",
point_estimate = obs_fit
)
conf_ints
#
# 仮想の個人個人のデータを作る
set.seed(137)
index <- sample(1:nrow(df), size = 10000, replace = TRUE)
df2 <- df[index, ]
df2
#
# observed fit
obs_fit2 <- df2 |>
specify(scaled ~ year + grade + gender) |>
fit()
obs_fit2
#
# ブートストラップ法
set.seed(149)
boot_dist2 <- df2 |>
specify(scaled ~ year + grade + gender) |>
generate(reps = 1000, type = "bootstrap") |>
fit()
#
# 95% 信頼区間
conf_ints2 <- get_confidence_interval(
boot_dist2,
level = 0.95,
type = "percentile",
point_estimate = obs_fit2
)
conf_ints2
#
