我正在使用R进行非线性回归分析,遇到了误差传播的问题.具体地说,我使用predictNLS来估计我的预测的误差,但在某些情况下,误差范围超过了生物学上合理的阈值(例如,预测间隔高于1或低于0,这在我的生物学背景下是没有意义的).

下面是我代码的一个可复制的例子.它模拟指数关系,拟合非线性最小二乘(NLS)模型,然后使用predictNLS来估计预测误差.我的挑战是处理错误传播导致生物学上不可能的值的实例.

Reproducible Example:

# simulate exponential relationship
set.seed(123)
# generate random x values between 0 and 60
x <- runif(100, 0, 60)
y <- 1 - exp(-0.075 * x) * rnorm(100, 0.7, 0.1)

data = data.frame(sr= x, fipar = y)

# create a nls model to fit the data
model <- nls(fipar ~ 1 - exp(-a * sr), data = data, start = list(a = 0.001))

# create an observed and predicted dataframe
data$predicted <- predict(model, data)

library(ggplot2)
data %>%
  ggplot(aes(x = sr, y = fipar)) + 
  geom_point() +
  geom_line(aes(y = predicted), color = "red")

# estimate the errors using predictNLS
newdat = data.frame(sr = seq(1, 60, 1))

prediction_se <- predictNLS(model, newdata = newdat, interval = "prediction", type = 'response')

prediction_se$summary$sr <- newdat$sr

prediction_se$summary %>%
  ggplot(aes(x = sr, y = Prop.Mean.1)) + 
  ylim(0, 1.2) +
  geom_point() +
  geom_ribbon(aes(ymin = `Prop.2.5%`, ymax = `Prop.97.5%`), alpha = 0.2) +
  geom_hline(yintercept = 1)

A scatter plot in ggplot2 showing relationship between model input and output, with 95% confidence interval band. A horizontal line represents the point at which any predicted value or error, would be going outside what is biologically possible

我对此的自然react 是简单地将任何错误值1替换为1,但我猜这忽略了一堆统计假设.>是否有任何推荐的方法来限制预测区间或调整误差估计以反映生物学约束?也许我的问题会更好地在统计堆栈交换,但我想我会寻找一个答案在这里太.

推荐答案

使用logit link函数的glm工作:

# Fit the GLM model 
model <- glm(fipar ~ sr, 
             data = data, 
             family = binomial(link = "logit"))

# Create an observed and predicted dataframe
data$predicted <- predict(model, data, type = "response",se.fit = T)$fit
data$se <- predict(model, data, type = "response",se.fit = T)$se.fit

# Plot the data and the model's predicted probabilities
library(ggplot2)
ggplot(data, aes(x = sr, y = fipar)) + 
  geom_point() +
  geom_line(aes(y = predicted), color = "red")+
  geom_ribbon(aes(ymin = predicted - se, ymax = predicted + se), alpha = 0.2)

enter image description here

R相关问答推荐

用单个表达匹配多个替代模式

使用ggplot将平滑线添加到条形图

Highcharter多次钻取不起作用,使用不同方法

如何在R中合并和合并多个rabrame?

迭代通过1个长度的字符串长字符R

是否可以创建一个ggplot与整洁判断的交互作用

如何自定义3D散点图的图例顺序?

合并DFS列表并将索引提取为新列

多个模拟序列间的一种预测回归关系

LOF中的插图短文字幕

有没有办法使用ggText,<;Sub>;&;<;sup>;将上标和下标添加到同一元素?

如何识别倒排的行并在R中删除它们?

按时间顺序对不同事件进行分组

调换行/列并将第一行(原始数据帧的第一列)提升为标题的Tidyr类似功能?

将全局环境变量的名称分配给列表中的所有元素

创建列并对大型数据集中的特定条件进行成对比较的更高效程序

如何将这个小列表转换为数据帧?

将字符变量出现次数不相等的字符框整形为pivot_wider,而不删除重复名称或嵌套字符变量

在具有条件的循环中添加行

从data.table列表中提取特定组值,并在R中作为向量返回