crosshyou

主にクロス表(分割表)分析をしようかなと思いはじめましたが、あまりクロス表の分析はできず。R言語の練習ブログになっています。

RのHistDataパッケージのChestSizes

f:id:cross_hyou:20220205175146j:plain

Photo by Robert Thiemann on Unsplash 

RのHistDataパッケージの中のChestSizesです。

Queteletという人が調べたスコットランドの軍人のチェストサイズ、胸囲ですかね?のデータセットです。

このデータをつかって自然のデータが正規分布になっていることをデモンストレーションしたそうです。

データを読み込みます。

f:id:cross_hyou:20220205175558p:plain

str()関数で確認すると、16行2列のデータフレームでしたので全部表示してみました。

chestがサイズでその人数がcountですね。

では、ヘルプのコードを実行してみます。

f:id:cross_hyou:20220205175805p:plain

f:id:cross_hyou:20220205175815p:plain

plot()関数をつかって横軸がchest, 縦軸がcountのグラフを描きました。type = "b" としているので、線と丸印の両方が描かれています。2列のデータフレームだと1列目が横軸で2列目が縦軸として描かれるんですね。

試しに適当なデータでやってみます。

f:id:cross_hyou:20220205180308p:plain

f:id:cross_hyou:20220205180317p:plain

予想どおり2列のデータフレームだと、1列目が横軸、2列目が縦軸になりました。

では、3列のデータフレームではどうでしょうか?

f:id:cross_hyou:20220205180754p:plain

f:id:cross_hyou:20220205180811p:plain

3列だと散布図マトリックスになるんですね。

脱線してしまいました。ヘルプのコードに戻ります。

f:id:cross_hyou:20220205181004p:plain

f:id:cross_hyou:20220205181017p:plain

barplot()関数で棒グラフを描いています。正規分布っぽいかたちですよね。

次のヘルプのコードは、

f:id:cross_hyou:20220205181151p:plain

です。cestが正規分布するとしたらどういう分布になるか、期待値を計算しています。

ここまででヘルプのコードは終わっています。

この期待値と実際の値をくらべましょう。

f:id:cross_hyou:20220205181500p:plain

f:id:cross_hyou:20220205181510p:plain

黒の実線が実際の値で赤の点線が期待値です。非常によく似た形状ですよね。

もう一つ、barplotに期待値を赤の点線で重ねてみましょう。

f:id:cross_hyou:20220205181748p:plain

f:id:cross_hyou:20220205181758p:plain

barplot()関数で生成された棒グラフは、横軸の座標がよくわからないので、いったんbplotという名前のオブジェクトとして保存しています。

このbplotというオブジェクトをstr()関数、class()関数で調べると、16行1列のマトリックスであることがわかりました。そしてこの値が横軸の座標でした。

そこで、lines()関数で横軸をbplot, 縦軸をexpectedにして棒グラフに赤い点線を載せました。

これでおしまいです。

以下に今回のコードを記載しておきます。

#
library(HistData)
data("ChestSizes")
#
str(ChestSizes)
ChestSizes
#
# frequency polygon
plot(ChestSizes, type = "b")
#
# 適当なデータで実験
plot(data.frame(number = 1:16, count = rnorm(16)), type = "b")
#
# 3列で実験
plot(data.frame(r1 = 1:16, r2 = rnorm(16), r3 = rep(1:4, 4)))
#
# barplot
barplot(ChestSizes[ , 2], ylab = "Frequency", xlab = "Cest size")
#
# calculate expected frequencies under normality, 
# chest ~ N(xbar, std)
n_obs <- sum(ChestSizes$count)
n_obs
xbar <- with(ChestSizes, weighted.mean(chest, count))
xbar
std <- with(ChestSizes, sd(rep(chest, count)))
std
expected <- with(ChestSizes,
                 diff(pnorm(c(32, chest) + 0.5, xbar, std))*sum(count))
expected
#
# plot actual data and expected data
plot(ChestSizes, type = "b", lwd = 1)
lines(ChestSizes$chest, expected, lty = 2, col = 2)
#
# barplot objects
bplot <- barplot(ChestSizes[ , 2], xlab = "chest", ylab = "frequency")
str(bplot)
class(bplot)
bplot
lines(bplot, expected, col = "red", lwd = 3, lty = 2)