最初にテーブルをアップロードします。表には 9 行が含まれており、そのうちの 6 行は要因であり、残りの 3 行は 152 人の個人 (n01、n02、n03) の成長率の個別の測定値です。次に、要因を指定します。
`r$feed <- factor (r$feed)`
`r$ph <- factor (r$ph)`
`r$aq <- factor (r$aq)`
`r$ind <- factor (r$ind)`
`r$wc <- factor (r$wc)`
`r$p0<- factor (r$p0)`
次に、興味のある要素を含む新しいテーブル "r2" にデータフレームを溶かし、na.omit 関数で NA 値を削除します。
`r2 <- data.table::melt(r,id.vars=c("feed","ph","aq","wc"),
measure=c("n01","n12","n23"),
variable.name="time",value.name="G")`
`r2<-na.omit(r2)`
r2 は次のようになります。
data.frame(
G = c(0.184, 0.087, 1.747, 0.11, 0.39, 0.062, 0.08, 0.189, 0.068,
0.262, 0.048, 0.029, 0, 0.229, 0.175),
feed = as.factor(c("HF", "HF", "HF", "HF", "HF", "HF", "HF", "HF",
"HF", "HF", "HF", "HF", "HF", "HF", "HF")),
ph = as.factor(c("8.1", "8.1", "8.1", "8.1", "8.1", "8.1", "8.1",
"8.1", "8.1", "8.1", "8.1", "8.1", "8.1", "8.1",
"8.1")),
aq = as.factor(c("1", "1", "1", "1", "1", "1", "2", "2", "2", "2",
"2", "2", "2", "3", "3")),
wc = as.factor(c("3", "3", "2", "3", "2", "4", "3", "4", "2", "2",
"3", "3", "1", "4", "3")),
time = as.factor(c("n01", "n01", "n01", "n01", "n01", "n01", "n01",
"n01", "n01", "n01", "n01", "n01", "n01", "n01",
"n01"))
)
その後、固定分散を設定し、次のように 2 つの gls モデルを適用して実行します。
`vfix3 <- varIdent(form=~1|time*factor(aq))
mix1 <- gls(G ~ ph+feed, weights=vfix3,data=r2)
mix3 <- gls(G ~ ph+feed+wc+time, weights=vfix3,data=r2)`
モデルの要約と anova を取得できるため、モデルは正常に機能しているようです。次に、次のように、パッケージ emmeans の lsmeans 関数を使用して事後ペアワイズ比較を実行しようとします。
print(lsmeans(mix1, list(pairwise~ph|feed), adjust="tukey"))
lsmeans は、2 因子モデル mix1 で正常に機能するようです。ただし、モデル mix3 で lsmeans を実行すると、次のエラーが表示されます。
crossprod(x, y) のエラー: 数値/複素行列/ベクトルの引数が必要です
モデルを行列に変換しようとしましたが、lsmeans関数に対して正しいオブジェクトではありません。また、係数を設定せずに列を数値のままにしてみましたが、同じエラーが表示されます。lsmeans 関数について読んでいるときに、それに関連する crossprod 関数が見つかりません。