Photo by Alexander Schimmeck on Unsplash
RのHistDataパッケージのCushnyPeeblesのデータは、Cushnyという方とPeeblesという方が睡眠に関する研究をしたときのデータのようです。
まずはデータを呼び出します。
ヘルプにないですが、このCushnyPeeblesのデータを表示してみます。
11行 x 4列のデータです。11行目のControlだけがNAです。
ヘルプのコードではplot()関数で散布図を描いています。
Control以外の3つの変数は結構な相関にみますね。
cor()関数で確認してみます。
やっぱり大きな相関係数ですね。
ヘルプのコードに戻ります。boxplot()関数で箱ひげ図を作っています。
Controlが一番、Hours of Sleepが短いですね。このデータの数値は睡眠時間なのですね。
ヘルプには無いですが、summary()関数でそれぞれの平均値をみてみます。
ヘルプに戻ります。白状しますと、以下のコードは何をしているのかよくわからないです。
lm()関数で線形モデルを作っていいると思うのですが、y ~ a + b + c ... の形式でなくて、yの部分がcbind関数でマトリックスになっていて、 ~ の右側が1なので切片だけ、つまり平均値を計算しているのですかね。。。ヘルプは無いのですが、結果をみてみます。
なんと、lm()関数で、y ~ a + b1x + b2x + ... の線形モデルのyのほうをマトリックスにすると、それぞれの列をyにして同時に複数のモデルを分析してくれるのですね。。
ヘルプのコードに戻ります。
carパッケージのAnova()関数でANOVA分析をしているのですね。
普通にanova()関数を使ったらどうなるでしょうか？
ヘルプのコードに戻ります。
次は、heplotsというパッケージを使っています。
そして最後はpairs()関数をつかっています。
これもよくわからないです。
これでCushnyPeeblesのヘルプのコードは終了です。
今回の発見は、lm()関数のモデルで y ~ x のyが複数の変数でも大丈夫ということでした。
以下が今回のコードです。
# HistDara_CushuyPeebles
library(HistData)
#
data(CushnyPeebles)
# CushnyPeeblesを表示
CushnyPeebles
# quick looks at the data
plot(CushnyPeebles)
# 相関係数
cor(CushnyPeebles[-11, ])
boxplot(CushnyPeebles, ylab = "Hours of Sleep", xlab = "Treatment")
# サマリー関数
summary(CushnyPeebles)
##########################
# Repeated measures MANOVA
CPmod <- lm(cbind(Control, L_hyoscyamine, L_hyoscine, DL_hyoscine) ~ 1,
data = CushnyPeebles)
# CPmodの結果を表示
summary(CPmod)
# Assign within-S factor and contrasts
Treatment <- factor(colnames(CushnyPeebles), levels = colnames(CushnyPeebles))
contrasts(Treatment) <- matrix(
c(-3, 1, 1, 1,
0,-2, 1, 1,
0, 0,-1, 1), ncol = 3)
colnames(contrasts(Treatment)) <- c("Control.Drug", "L.DL", "L_hy.DL_hy")
Treats <- data.frame(Treatment)
if (require(car)) {
(CPaov <- Anova(CPmod, idata = Treats, idesign = ~Treatment))
}
# ただのanova()だとどうなるか？
anova(CPmod)
#
summary(CPaov, univariate = FALSE)
library(heplots)
heplot(CPmod, idata = Treats, idesign = ~Treatment, iterm = "Treatment",
xlab = "Control vs Drugs", ylab = "L vs DL drug")
pairs(CPmod, idata = Treats, idesign = ~Treatment, iterm = "Treatment")
################################
# reshape to long format, add Ns
CPlong <- stack(CushnyPeebles)[,2:1]
colnames(CPlong) <- c("treatment", "sleep")
CPN <- stack(CushnyPeeblesN)
CPlong <- data.frame(patient = rep(1:11,4), CPlong, n = CPN$values)
str(CPlong)
