Как строятся регрессионные модели, доверительные интервалы и данные? - PullRequest
0 голосов
/ 05 ноября 2018

У меня есть следующая модель. Файл данных находится в: https://drive.google.com/open?id=1_H6YZbdesK7pk5H23mZtp5KhVRKz0Ozl

library(nlme)
library(lme4)
library(car)
library(carData)
library(emmeans)
library(ggplot2)
library(Matrix)
library(multcompView)
datos_weight <- read.csv2("D:/investigacion/publicaciones/articulos-escribiendo/pennisetum/pennisetum-agronomicas/data_weight.csv",header=T, sep = ";", dec = ",")

parte_fija_3 <- formula(weight_DM 
                    ~ Genotypes 
                    + Age 
                    + I(Age^2) 
                    + Genotypes*Age 
                    + Genotypes*I(Age^2))
heterocedasticidad_5 <- varComb(varExp(form = ~fitted(.)))
correlacion_4 <- corCompSymm(form = ~ 1|Block/Genotypes)

modelo_43 <- gls(parte_fija_3, 
             weights = heterocedasticidad_5, 
             correlation = correlacion_4, 
             na.action = na.omit, 
             data = datos_weight)
anova(modelo_43)

#response
Denom. DF: 48 
                   numDF  F-value p-value
(Intercept)            1 597.3828  <.0001
Genotypes              3   2.9416  0.0424
Age                    1 471.6933  <.0001
I(Age^2)               1  22.7748  <.0001
Genotypes:Age          3   5.9425  0.0016
Genotypes:I(Age^2)     3   0.7544  0.5253

Теперь я хочу построить регрессионные модели с доверительными интервалами и данными, разделенными для каждого генотипа. Я использовал ggplot2, и у меня есть данные, я не смог добавить регрессионные модели с доверительными интервалами.

library(ggplot2)
rango_X <- c(30,90) #x axis
rango_Y <- c(0,175) #y axis
ggplot(datos_weight, aes(x = Age, y = weight_DM)) +
  geom_point() + 
  xlab("Age") + 
  ylab("Dry matter") +
  xlim(rango_X) +
  ylim(rango_Y) +
  facet_wrap(~ Genotypes, ncol = 2)

График выглядит следующим образом:

enter image description here

Для следующего анализа тех же данных, где нет взаимодействия с квадратичным возрастом: Genotypes*I(Age^2), как бы вы добавили регрессионные модели с доверительными интервалами к графику?

parte_fija_3 <- formula(weight_DM 
                        ~ Genotypes 
                        + Age
                        + I(Age^2)
                        + Genotypes*Age) 
                        #+ Genotypes*I(Age^2))
> anova(modelo_44)
Denom. DF: 51 
              numDF  F-value p-value
(Intercept)       1 609.3684  <.0001
Genotypes         3   3.7264  0.0169
Age               1 479.0973  <.0001
I(Age^2)          1  21.9232  <.0001
Genotypes:Age     3   6.4184  0.0009

Линейные склоны от modelo_44:

(tendencias_em_lin <- emtrends(modelo_44,
                                "Genotypes",
                                var = "Age"))
Genotypes Age.trend        SE df lower.CL upper.CL
C          1.613619 0.1723451 51 1.267622 1.959616
E          1.665132 0.2024104 51 1.258776 2.071488
K          1.888587 0.2001627 51 1.486744 2.290430
M          1.059897 0.1205392 51 0.817905 1.301890

Квадратичные склоны есть?

(tendencias_em_quad <- emtrends(modelo_44,
                                "Genotypes",
                                var = "I(Age^2)"))
 Genotypes I(Age^2).trend           SE df    lower.CL   upper.CL
 C            0.013379926 0.0014290639 51 0.010510961 0.01624889
 E            0.013807066 0.0016783618 51 0.010437614 0.01717652
 K            0.015659927 0.0016597235 51 0.012327893 0.01899196
 M            0.008788536 0.0009994958 51 0.006781965 0.01079511

Confidence level used: 0.95 

Или сумасшедший из резюме: I(Age^2) = -0.01511? Я считаю, что наклон постоянен для всех генотипов, потому что взаимодействие Genotypes*I(Age^2) не было проверено в modelo_44:

summary(modelo_44)
Generalized least squares fit by REML
Model: parte_fija_3 
....
Coefficients:
                   Value Std.Error   t-value p-value
(Intercept)    -73.32555 11.236777 -6.525496  0.0000
GenotypesE       7.22267  9.581979  0.753776  0.4544
GenotypesK      -9.83285  9.165962 -1.072757  0.2884
GenotypesM      17.87000  8.085229  2.210203  0.0316
Age              3.43593  0.450041  7.634687  0.0000
I(Age^2)        -0.01511  0.004065 -3.717475  0.0005
GenotypesE:Age   0.05151  0.246724  0.208788  0.8354
GenotypesK:Age   0.27497  0.241923  1.136595  0.2610
GenotypesM:Age  -0.55372  0.195398 -2.833808  0.0066
...

Вопросы

  1. Как добавить регрессионные модели с доверительными интервалами и данными для каждого генотипа в отдельные графики, такие как представленные, с ggplot2 или другим вариантом, если мне пришлось строить графики для моделей: modelo_43 и modelo_44
  2. Правильно ли я рассчитал оценку квадратичного уклона с emtrends для modelo_44, как это правильно?

Большое спасибо за ответ

1 Ответ

0 голосов
/ 05 ноября 2018

Этот вопрос выглядит смутно знакомым - вы уже публиковали его раньше? Может быть, я думаю о каком-то подобном вопросе от кого-то еще.

Кажется, вы пытаетесь нанести яблоки, апельсины и бананы в одном масштабе. Я не уверен, в каких единицах реагирует переменная (сухое вещество находится в); скажем, его в кг. Затем результаты в tendencies_em_lin приведены в кг в год, а результаты в tendencies_em_quad - в кг в год ^ 2. Это три разные шкалы, и нет смысла «показывать данные» на графиках тех.

То, что я думаю, имеет смысл сделать примерно так:

emm <- emmeans(modelo_44, ~ Genotype*Age,
    at = list(Age = seq(from = 40, to = 80, by = 5)))

Это позволит получить прогнозы для данного возраста с каждым генотипом. Теперь вы можете построить их следующим образом:

plotobj <- emmip(emm, Genotype ~ Age, CIs = TRUE)
plotobj

Возвращенный plotobj является объектом ggplot, к которому вы можете добавить данные, используя методы, как показано в примере в конце графического раздела в https://cran.r -project.org / web / packages / emmeans / виньетки / basics.html # участки .

Или, вы можете использовать dat = as.data.frame(emm) в качестве фрейма данных, содержащего результаты, которые вам нужны, и строить их так, как вам нравится. Опять же, вы можете использовать методы ggplot2 , чтобы добавить данные наблюдений к этим графикам.

В любом случае, линейные тренды будут видны как увеличение или уменьшение в построенных EMM, а квадратичные тренды будут видны как кривизна на этих трассах.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...