HistDataパッケージのDrinksWagesのデータセットは、Pearsonが両親のアルコール飲酒と子どもの賃金に関するデータセットです。Pearsonは両親がいっぱいアルコール飲酒している子どもの賃金は低いということを実証しようとしていたようですが、実際は失敗の終わり、両親のアルコール摂取と子どもの賃金は関係ないようです。
それではヘルプのコードをやってみましょう。
まずはデータを呼び出します。
ここでヘルプには無いですが、str()関数とhead()関数でデータの様子を確認しておきます。
70行、6列のデータフレームです。classというのは賃金クラスとのことです。
tradeは職業のようです。soberは飲酒しない人の人数、drinksは飲酒する人の人数、wageは賃金、totalは全体の人数です。sober + drinks = total という関係性が成り立ちます。
ヘルプのコードはplot()関数でデータ全体をプロットしています。
次は、plot()関数でsoberとwage, classをプロットしています。
横軸がwage, 賃金で、縦軸がsober/n、全体の人数なので全体の人数に対する飲酒しない人の比率ですね。青がclass A, 赤が class B, 緑が class C です。Cクラスのほうが賃金が高いようです。
Pearsonは恐らく、飲酒しない人の比率が高いほうが賃金が高い、ということを言いたかったのでしょうが、このグラフを見る限りではそんなことはないようです。
ヘルプでは続いてglm()関数でロジスティクス回帰分析をしています。
飲酒しない人の比率を賃金で回帰分析しているようですが、wageの係数のp値は0.88443となっていて、wageは関係ない、ということがわかります。
そして最後はplot()関数でモデルのグラフを描いています。
ヘルプには無いですが、wageをsober/nとclassで回帰分析してみます。
I(sober/n)の係数が2.0520でp値が0.0339なので5%水準で有意な値です。
soberの比率が0.01上昇すると、賃金は0.02上昇するといます。Personの狙い通りの結果ですよね。
今回は以上です。
下には今回のコードを記載しておきます。
# 2022-03-05
# HistData DrinksWages
library(HistData)
data(DrinksWages)
#
# data structure
str(DrinksWages)
#
# view a few rows of data
head(DrinksWages)
#
# plot data
plot(DrinksWages)
#
# plot proportion sober vs. wage | class
with(DrinksWages, plot(wage, sober/n,
col = c("blue", "red", "green")[class]))
#
# fit logistic regression model of sober on wage
mod.sober <- glm(cbind(sober, n) ~ wage, family = binomial,
data = DrinksWages)
summary(mod.sober)
#
# plot model
par(mfrow = c(2, 2))
plot(mod.sober)
par(mfrow = c(1, 1))
#
# wage = beta0 + beta1*(sober/n) + class
summary(lm(wage ~ I(sober/n) + class, data = DrinksWages))
