2

[背景として私が持っている実験の詳細を説明していますlmer-sの方法については明確ですが、必要な値を抽出する/手動で計算する方法が不明であるため、これをCVではなくSOに投稿しました。これが正しい投稿場所であったことを願っています!]

データはこちら。

私の実験には、ブロック/プロット/サブプロットのレベルを持つ分割プロット デザインがあります。

6ブロックあります。各ブロックには 2 つのプロットがあり、各プロットには 2 つのサブプロットがあります。処理 1 には 2 つのレベル (A と B) があり、プロット レベルで適用されます。各ブロックには、処理 1 レベル A を受け取る 1 つのプロットと、処理 1 レベル B を受け取る 1 つのプロットがあります。

処理 2 はサブプロット レベルで適用され、さらに 2 つのレベル (C および D) があります。各プロットには、処理 2 レベル A を受け取る 1 つのサブプロットと、処理 2 レベル B を受け取る 1 つのサブプロットがあります。

実験は 3 年間行われた。2 つの治療法の各組み合わせが従属変数 (DV) にどのように影響するかに興味があります。

そのため、4つの治療の組み合わせがあります。

TMT1A:TMT2C

TMT1B:TMT2C

TMT1A:TMT2D

TMT1b:TMT2D

分割プロットの設計を考慮して、モデルに lmer を使用しています。私はクロスイヤーモデルを実行していますが、各年のモデルも順番に実行しています（実験での複製では、クロスイヤーモデルでの年の効果のテストが許可されていないため、モデルは過剰にパラメータ化されていました）。

各年のlmers は次のようになります。

m2011<- lmer (DV2011~ TMT1*TMT2 + (1|Block/TMT1))
m2012<- lmer (DV2012~ TMT1*TMT2 + (1|Block/TMT1))
m2013<- lmer (DV2013~ TMT1*TMT2 + (1|Block/TMT1))

時間の経過に伴うこれらの治療手段の変化をグラフで表現するために、毎年の各治療の各レベル (上記の 4 つのレベルを参照) の治療平均を抽出し、これらを実験の年ごとにプロットします。この投稿の例

オブジェクトから 4 つの異なる処理の組み合わせ (上記のリストなど) の処理手段を抽出することは可能lmerですか? それとも、手で計算する必要がありますか?

私が考えた 1 つの方法は、4 つの治療の組み合わせを表す別の因子を実際に作成することです (貼り付けたデータの列 "TMT1x2" を参照)。次に、毎年次のモデルを実行できます。

m2011<- lmer (DV2011~ TMT1x2 + (1|Block/TMT1))

4 つのレベルのそれぞれの処理手段をそのように抽出します。ただし、この新しい 4 レベル ファクターは、それを構成するレベルのネストされた性質を無視するため (ランダム効果はそれを無視しませんが)、この方法が分割プロット デザインを適切に制御するかどうかはわかりません...

さらに、手動で処理手段を計算する必要がある場合、私の実験で入れ子のレベルを考慮してこれを行う方法を知っている人はいますか?

また、これらの各処理手段のエラーバーを計算したいと思います...

誰かがこれについて何か洞察を持っているなら、それは大歓迎です!

4

3 に答える 3

2

あなたが求めているのは、クラスのデフォルトメソッドがない（少なくともCRANのバージョン）の何らかの形式だと思いpredict()ます。ただし、使用できます。merlme4ez::ezPredict

library(ez)
library(ggplot2)
to_predict <- expand.grid(TMT1=c("A","B"), TMT2=c("C","D"))
t_means <- rbind(ezPredict(m2011, to_predict=to_predict, boot=F), ezPredict(m2012, to_predict=to_predict, boot=F), ezPredict(m2013, to_predict=to_predict, boot=F) )
t_means$YEAR = rep(2011:2013, each = 4)
ggplot(t_means, aes(x=YEAR, y=value, color=TMT1:TMT2)) + geom_point() + geom_line()

この関数には、ブートストラップされた値を提供するなど、役立つ可能性のある追加機能がいくつかあります。

処理平均の点推定だけが必要な場合は、特に 3 つのモデルすべてが同じ計画行列を持っているため、手動で計算するのと同じくらい簡単です。

mm = unique(model.matrix(m2011))
Y_bar <- c(mm%*%fixef(m2011), mm%*%fixef(m2012), mm%*%fixef(m2013))
ggplot(t_means, aes(x=YEAR, y=Y_bar, color=TMT1:TMT2)) + geom_point() + geom_line()

「処理手段を計算する...私の実験での入れ子のレベルを説明する」という意味が正確にはわかりません。混合モデルの変量効果は、母集団レベルの効果 (固定効果) からの構造化された正規分布偏差ranef(m2011)ですm2011@Zt。

したがって、プロットしたいのが母集団レベルの処理平均だけである場合は、固定効果の推定値とfixef(m2011)固定効果計画行列model.matrix(m2011)を上記のように操作するだけで済みます。母集団レベルの予測に何らかの不確実性の尺度を含めたい場合、または各ブロック/プロット/サブプロットの予測が必要な場合は、ランダム効果と固定効果の両方を使用する必要があります。http://glmm.wikidot.com/faqの「予測および/または予測の信頼 (または予測) 間隔」という見出しの下を見ることから始めることをお勧めします。

編集 2013 年 8 月 26 日:

bootMer()for (パラメトリック ブートストラップされた) 信頼区間の開発バージョンでは、ランダム効果の分散に不確実性を組み込む必要lme4があり、GLMM で機能する必要があります (たとえば、このスレッドを参照してください)。

アイデアは、対象のモデルからシミュレートし、シミュレートされた値で再適合し、再適合モデルから対象の統計を計算することです。simulate()と を使用して、自分で手順を実行できますrefit()。

t_sim <- apply(simulate(model, 999), 2, function(x) combn(unique(model.matrix(model))%*%fixef(refit(model, x)), 2, diff) )

これは、処理手段間のペアごとの差の 999quantile()個のブートストラップ担当者を生成します。

apply(t_sim, 1, function(.) quantile(., c(0.975, 0.025)))
于 2013-08-25T20:51:23.523 に答える
2

package の関数を使用する代替手段languageR。私はあなたのデータセットを呼び出しますdf。

library(lme4)
library(languageR)
library(ggplot2)

# fit model
# n.b. I don't claim that this is a sensible model
# It is just used to demonstrate the plot
mod <- lmer(DV ~ TMT1 * TMT2 + (1|Block), data = df)

# create MCMC matrix
mcmc <- pvals.fnc(mod, nsim = 1000, withMCMC = TRUE)
# pval.fnc also calculates MCMC-based p-values and HPD confidence intervals,
# and plot the posterior distributions of the parameters

# plot using plotLMER.fnc 
# in addition, set withList = TRUE to create a list of data frames with plot data
# which can be used for a (possibly prettier) plot in ggplot
ll <- plotLMER.fnc(mod, withList = TRUE, pred = "TMT1", 
               intr = list(
                 "TMT2",
                 c("C", "D"),
                 "end",
                 list(c("red",  "blue"), rep(1, 2))),
               addlines = TRUE,
               mcmcMat = mcmc$mcmc)

 # here follows additional steps to plot using ggplot 

 # convert list to data frame
 df <- do.call(rbind, ll$TMT1)

 # rename 
 names(df)[names(df) == "Levels"] <- "TMT1"

 # add TMT2
 df$TMT2 <- rep(c("C", "D"), each = 2)

# plot using ggplot
dodge <- position_dodge(width = 0.1)
ggplot(data = df, aes(x = TMT1, y = Y, col = TMT2, group = TMT2)) +
   geom_point(position = dodge, size = 3) +
   geom_errorbar(aes(ymax = upper, ymin = lower, width = 0.1), position = dodge) +
   geom_line(position = dodge) +
   ylab("DV") +
   theme_classic()
于 2013-08-26T17:54:20.210 に答える