www.crosshyou.info

政府統計の総合窓口のデータや、OECDやUCIやのデータを使って、Rの練習をしています。ときどき、読書記録も載せています。

RのHistDataパッケージのCushnyPeebles

f:id:cross_hyou:20220220210205j:plain

Photo by Alexander Schimmeck on Unsplash 

RのHistDataパッケージのCushnyPeeblesのデータは、Cushnyという方とPeeblesという方が睡眠に関する研究をしたときのデータのようです。

まずはデータを呼び出します。

f:id:cross_hyou:20220220210605p:plain

ヘルプにないですが、このCushnyPeeblesのデータを表示してみます。

f:id:cross_hyou:20220220210742p:plain

11行 x 4列のデータです。11行目のControlだけがNAです。

ヘルプのコードではplot()関数で散布図を描いています。

f:id:cross_hyou:20220220210943p:plain

f:id:cross_hyou:20220220210955p:plain

Control以外の3つの変数は結構な相関にみますね。

cor()関数で確認してみます。

f:id:cross_hyou:20220220211211p:plain

やっぱり大きな相関係数ですね。

ヘルプのコードに戻ります。boxplot()関数で箱ひげ図を作っています。

f:id:cross_hyou:20220220211354p:plain

f:id:cross_hyou:20220220211406p:plain

Controlが一番、Hours of Sleepが短いですね。このデータの数値は睡眠時間なのですね。

ヘルプには無いですが、summary()関数でそれぞれの平均値をみてみます。

f:id:cross_hyou:20220220211626p:plain

ヘルプに戻ります。白状しますと、以下のコードは何をしているのかよくわからないです。

f:id:cross_hyou:20220220211756p:plain

lm()関数で線形モデルを作っていいると思うのですが、y ~ a + b + c ... の形式でなくて、yの部分がcbind関数でマトリックスになっていて、 ~ の右側が1なので切片だけ、つまり平均値を計算しているのですかね。。。ヘルプは無いのですが、結果をみてみます。

f:id:cross_hyou:20220220212234p:plain

f:id:cross_hyou:20220220212257p:plain

f:id:cross_hyou:20220220212307p:plain

なんと、lm()関数で、y ~ a + b1x + b2x + ... の線形モデルのyのほうをマトリックスにすると、それぞれの列をyにして同時に複数のモデルを分析してくれるのですね。。

ヘルプのコードに戻ります。

f:id:cross_hyou:20220220212602p:plain

carパッケージのAnova()関数でANOVA分析をしているのですね。

普通にanova()関数を使ったらどうなるでしょうか?

f:id:cross_hyou:20220220213324p:plain

ヘルプのコードに戻ります。

f:id:cross_hyou:20220220213518p:plain

f:id:cross_hyou:20220220213530p:plain

次は、heplotsというパッケージを使っています。

f:id:cross_hyou:20220220213736p:plain

f:id:cross_hyou:20220220213753p:plain

そして最後はpairs()関数をつかっています。

f:id:cross_hyou:20220220213920p:plain

f:id:cross_hyou:20220220213931p:plain

これもよくわからないです。

これで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)