【问题标题】:Plotting a list of non linear regressions with ggplot使用 ggplot 绘制非线性回归列表
【发布时间】:2016-09-04 12:21:17
【问题描述】:

作为此链接的非线性回归分析的输出图

https://stats.stackexchange.com/questions/209087/non-linear-regression-mixed-model

有了这个数据集:

zz <-(" iso temp diam
 Itiquira   22  5.0
 Itiquira   22  4.7
 Itiquira   22  5.4
 Itiquira   25  5.8
 Itiquira   25  5.4
 Itiquira   25  5.0
 Itiquira   28  4.9
 Itiquira   28  5.2
 Itiquira   28  5.2
 Itiquira   31  4.2
 Itiquira   31  4.0
 Itiquira   31  4.1
 Londrina   22  4.5
 Londrina   22  5.0
 Londrina   22  4.4
 Londrina   25  5.0
 Londrina   25  5.5
 Londrina   25  5.3
 Londrina   28  4.6
 Londrina   28  4.3
 Londrina   28  4.9
 Londrina   31  4.4
 Londrina   31  4.1
 Londrina   31  4.4
    Sinop   22  4.5
    Sinop   22  5.2
    Sinop   22  4.6
    Sinop   25  5.7
    Sinop   25  5.9
    Sinop   25  5.8
    Sinop   28  6.0
    Sinop   28  5.5
    Sinop   28  5.8
    Sinop   31  4.5
    Sinop   31  4.6
    Sinop   31  4.3"
)
df <- read.table(text=zz, header = TRUE)

而这个拟合模型,有四个参数:

thx:最佳温度

thy:最佳直径

thq:曲率

thc:偏度

library(nlme) 

df <- groupedData(diam ~ temp | iso, data = df, order = FALSE) 

n0 <- nlsList(diam ~ thy * exp(thq * (temp - thx)^2 + thc * (temp - thx)^3),               
      data = df, 
      start = c(thy = 5.5, thq = -0.01, thx = 25, thc = -0.001))

> n0
# Call:
#  Model: diam ~ thy * exp(thq * (temp - thx)^2 + thc * (temp - thx)^3) | iso 

# Coefficients:
              thy          thq      thx           thc
# Itiquira 5.403118 -0.007258245 25.28318 -0.0002075323
# Londrina 5.298662 -0.018291649 24.40439  0.0020454476
# Sinop    5.949080 -0.012501783 26.44975 -0.0002945292

# Degrees of freedom: 36 total; 24 residual
# Residual standard error: 0.2661453

有没有办法在 ggplot 中绘制拟合值,比如 smooth() 的特定函数?

我想我找到了...(基于http://rforbiochemists.blogspot.com.br/2015/06/plotting-two-enzyme-plots-with-ggplot.html)

ip <- ggplot(data=daf,  aes(x=temp, y=diam, colour = iso)) +  
  geom_point() + facet_wrap(~iso)

ip + geom_smooth(method = "nls", 
                method.args = list(formula = y ~ thy * exp(thq * (x-thx)^2 + thc * (x - thx)^3), 
                                   start = list(thy=5.4, thq=-0.01, thx=25, thc=0.0008)),
                se = F, size = 0.5, data = subset(daf, iso=="Itiquira")) +

  geom_smooth(method = "nls", 
              method.args = list(formula = y ~ thy * exp(thq * (x-thx)^2 + thc * (x - thx)^3), 
                                 start = list(thy=5.4, thq=-0.01, thx=25, thc=0.0008)),
              se = F, size = 0.5, data = subset(daf, iso=="Londrina")) +

  geom_smooth(method = "nls", 
              method.args = list(formula = y ~ thy * exp(thq * (x-thx)^2 + thc * (x - thx)^3), 
                                 start = list(thy=5.4, thq=-0.01, thx=25, thc=0.0008)),
              se = F, size = 0.5, data = subset(daf, iso=="Sinop")) 

【问题讨论】:

  • 您可以使用predict 和geom_line 连接变量的最小值和最大值之间的大量点吗?这是否会给您所需的结果?
  • 您能提供一些链接作为示例吗?
  • 应该不难,但你能提供一些数据来运行它吗?没有它就不是一个完整的问题。您可以使用 dput(df) 并将输出添加到您的问题中。
  • @dww 数据在上面的链接中
  • 数据的正确位置在问题中。如果您希望人们帮助您,那么您必须努力使他们尽可能轻松,这才是公平的。另外,如果链接的问题被删除会怎样?

标签: r list ggplot2 nlme


【解决方案1】:

用更原则的ggplot 方法来回答这个问题(将输出组合成一个结构与原始数据匹配的单个数据帧)。不幸的是,找到nls 预测的置信区间并不容易(搜索涉及引导或增量方法的解决方案):

tempvec <- seq(22,30,length.out=51)
pp <- predict(n0,newdata=data.frame(temp=tempvec))
## combine predictions with info about species, temp
pdf <- data.frame(iso=names(pp),
                  temp=rep(tempvec,3),
                  diam=pp)

创建图表:

library(ggplot2)
ggplot(df,aes(temp,diam,colour=iso))+
  stat_sum()+
  geom_line(data=pdf)+
  facet_wrap(~iso)+
  theme_bw()+
  scale_size(range=c(1,4))+
  scale_colour_brewer(palette="Dark2")+
  theme(legend.position="none",
        panel.spacing=grid::unit(0,"lines"))

【讨论】:

    猜你喜欢
    • 2021-05-31
    • 2019-03-06
    • 2019-10-20
    • 2017-01-21
    • 2016-07-09
    • 1970-01-01
    • 1970-01-01
    • 2021-10-12
    • 2017-02-21
    相关资源
    最近更新 更多