时间戳

首页 > 解决方案 > 回归模型、置信区间和数据是如何绘制的？

r - 回归模型、置信区间和数据是如何绘制的？

问题描述

我有以下模型。数据文件在：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)

图表如下：

在此处输入图像描述

对于相同数据的下一次分析，与二次年龄没有交互：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. 我是否用emtrendsfor正确计算了二次斜率的估计值modelo_44，它是如何正确的？

非常感谢你的回复

标签： rggplot2plotregressionemmeans

解决方案

这个问题看起来有点眼熟——你以前发过吗？也许我正在考虑其他人的一些类似问题。

您似乎正在尝试以相同的比例绘制苹果、橙子和香蕉。我不确定响应变量（干物质的单位）是什么单位；假设它以公斤为单位。然后结果以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

返回的是一个 ggplot 对象，您可以使用https://cran.r-project.org/web/packages/emmeans/vignettes/basicsplotobj中图形部分末尾的示例中所示的技术向其中添加数据.html#plots。

或者，您可以将dat = as.data.frame(emm)其用作包含所需结果的数据框，并根据需要绘制它们。同样，您可以使用ggplot2技术将观察到的数据添加到这些图中。

无论哪种方式，线性趋势将随着绘制的 EMM 的增加或减少而可见，而二次趋势将作为这些路径中的曲率可见。

推荐阅读