оценочные предельные средние с моделью смешанных эффектов: выбор модели и интерпретация результатов

Question

Сначала я пишу свои вопросы, чтобы вы могли задавать их по ходу дела:

1. Правильно ли выбрать лучшую модель со значением p и AIC так же, какЯ использовал?

2. Являются ли степени свободы проблемой в моей модели, особенно в отношении моделей с взаимодействиями (...~ treatment * level...)?

3. почему результаты emmeans (см. Таблицу cld.mixed_C или рисунок внизу) показывают, что нет существенной разницы между уровнем $ up и уровнем $ mid, когда, по-видимому, существует резкая разница между stack$level (почему CL такие большие в cld.mixed_C?)?если бы я использовал model_lmer_Clog10 вместо model_lmer_Cinter, эти различия были бы учтены emmeans.Это потому, что предположения модели не выполнены?

4. Есть ли способ убрать отрицательный нижний CL на графике?

Я хочу сравнить три stack$treatment (W; C; F), содержащие три разных stack$level (вверх, середина, низ) для различий в stack$continuous.

stack <- structure(list(treatment = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L), .Label = c("W", "C", "F"), class = "factor"), 
    block = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
    1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L
    ), level = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
    1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 
    3L, 3L, 3L), .Label = c("up", "mid", "low"), class = "factor"), 
    continuous = c(2.71097672953743, 3.65651341057938, 2.54600229936132, 
    63.2884147056461, 47.6120381265922, 35.7709752585528, 138.618364821712, 
    138.364742211133, 248.629954375939, 3.91496531831404, 4.14426270796186, 
    3.90287641449911, 51.619396872644, 56.2744358488691, 40.8303737226699, 
    174.673343321081, 185.181589166907, 260.523764718931, 4.03192255652607, 
    3.09513270445542, 2.45237238277809, 61.645621226552, 75.8775130347202, 
    47.1098329758748, 347.373110867901, 209.39376817372, 254.575864308321
    )), row.names = c(NA, -27L), class = "data.frame")

Вот boxplot(continuous~treatment+level,data=stack, xlab="treatment", ylab="continuous"):

Я подгоняю следующую модель к данным.Здесь у меня есть случайная переменная, состоящая из treatment, вложенных в block (level вложена в treatment, это не указано, но может быть опущено).

Сначала я проверил, если взаимодействиебыла лучшей моделью, сравнив ее с более простой моделью без взаимодействия (также преобразования log10 и sqrt), затем я запустил таблицу anova:

library(lme4)
# no interactions:
model_lmer_Cbase <- lmer((continuous) ~  treatment + level + 
                                 (1|block/treatment),data= stack, REML = FALSE)
#interactions:
model_lmer_Cinter <- lmer((continuous) ~  treatment * level + 
                                  (1|block/treatment),data= stack, REML = FALSE)


anova_models_compare <- anova(model_lmer_Cbase, model_lmer_Cinter)

anova_models_compare:

                  Df      AIC      BIC    logLik deviance    Chisq Chi Df Pr(>Chisq)
model_lmer_Cbase   8 274.2591 284.6258 -129.1296 258.2591       NA     NA         NA
model_lmer_Cinter 12 271.6778 287.2279 -123.8389 247.6778 10.58129      4 0.03169548

Лучшая модель - model_lmer_Cbase (model_lmer_Cinter имеет значение p <.05, что означает значительно меньшую объяснительную способность, даже если его AIC ниже, чем <code>model_lmer_Cbase):

model_lmer_C <- lmer((continuous) ~  treatment + level + (1|block/treatment), 
                     data= stack, REML = TRUE)

Затем я создаю график с приблизительным средним значением:

library(multcompView)
library(emmeans)

lsm.mixed_C<-emmeans::emmeans(model_lmer_C,pairwise ~  treatment * level, type="response")

cld.mixed_C<-data.frame(CLD(lsm.mixed_C,alpha=0.05,Letters=letters,
                                adjust="tukey"))

Таблица cld.mixed_C выглядит следующим образом:

  treatment level   response       SE       df   lower.CL  upper.CL .group
1         W    up  10.176389 13.77472 14.28571 -34.590414  54.94319      a
2         W   mid  53.379809 13.77472 14.28571   8.613006  98.14661      a
3         W   low 276.898486 13.77472 14.28571 232.131682 321.66529      b
4         C    up   8.625707 13.77472 14.28571 -36.141097  53.39251      a
5         C   mid  51.829127 13.77472 14.28571   7.062323  96.59593      a
6         C   low 275.347804 13.77472 14.28571 230.581000 320.11461      b
7         F    up  -9.446419 13.77472 14.28571 -54.213222  35.32038      a
8         F   mid  33.757001 13.77472 14.28571 -11.009802  78.52380      a
9         F   low 257.275678 13.77472 14.28571 212.508875 302.04248      b

Затем я запускаю (довольно грязный, извиняюсь) кодчтобы отобразить приведенное выше парное сравнение treatment и level:

trim <- function (x) gsub("^\\s+|\\s+$", "", x) 
cld.mixed_C$.group <- trim(cld.mixed_C$.group)
# I want to rename some variables:
library(dplyr)

## sort factors in a specific order:
cld.mixed_C$level = factor(cld.mixed_C$level ,levels=c("up", "mid", "low"))
cld.mixed_C$treatment = factor(cld.mixed_C$treatment,levels=c("W", "C", "F"))
# order table according to specific order:
cld.mixed_C <- cld.mixed_C[with(cld.mixed_C, order(treatment, level)), ]

library(ggplot2)
cld.mixed_C$treatment = factor(cld.mixed_C$treatment, 
                               levels=c("W", "C", "F")) ### Order the levels for printing
cld.mixed_C$level = factor(cld.mixed_C$level,levels=c("up", "mid", "low"))
cld.mixed_C<- with(cld.mixed_C, cld.mixed_C[order(treatment, level),])
rownames(cld.mixed_C) <- NULL 
pd = position_dodge(0.6)    ### How much to jitter the points on the plot
plot.mixed.lme<-ggplot(cld.mixed_C,aes(x = treatment,y=response, color= level,
                                       label=.group))+ 
        theme_bw()+
        geom_point(shape  = 15, size   = 4, position = pd) +
        geom_errorbar(aes(ymin  =  lower.CL,ymax  =  upper.CL),width =0.2,
                      size = 1, position = pd) +
        theme(axis.title   = element_text(face = "bold"),
              axis.text    = element_text(face = "bold"),
              plot.caption = element_text(hjust = 0)) +
        #following lines adds legend inside the boxplot. 
        theme(legend.justification=c(1,0), legend.position=c(0.1,0.7),
              legend.direction="vertical",
              legend.box="vertical",
              legend.box.just = c("top"), 
              legend.background = element_rect(fill=alpha('white', 0.4)))+
        theme(plot.title = element_blank(), axis.title.x = element_blank()) +
        #remove grid lines (I only remove managem):
        theme(panel.grid.major.x = element_blank() ,
              panel.grid.major.y = element_blank(), 
              panel.grid.minor.x = element_blank(),
              panel.grid.minor.y = element_blank()) +
        ylab("Estimated marginal mean\ncontinuous") +
        theme(axis.text = element_text(size = 15)) +labs(hjust=0.5) +
        geom_text(nudge_x = c(-0.3, 0, 0.3, -0.3, 0, 0.3,-0.3, 0, 0.3),
                  nudge_y = c(40, 40, 40, 40, 40, 40, 40, 40, 40),
                  color   = "black", size= 6) +
        scale_color_manual(values = c("#333333", "#669966", "slategray3"))

plot.mixed.lme

Заранее большое спасибо.