计算并绘制广义非线性模型的 95% 置信区间
Calculate and plot 95% confidence intervals of a generalised nonlinear model
我用 R 包 nlme 和包含的 gnls() 函数构建了几个广义非线性最小二乘模型（指数衰减）。我不简单地使用基础 nls() 函数构建非线性最小二乘模型的原因是因为我希望能够对异方差建模以避免转换。我的模型看起来像这样：
model <- gnls(Response ~ C * exp(k * Explanatory1) + A,
start = list(C = c(C1,C1), k = c(k1,k1), A = c(A1,A1)),
params = list(C ~ Explanatory2, k ~ Explanatory2,
A ~ Explanatory2),
weights = varPower(),
data = Data)
与简单 nls() 模型的主要区别在于 weights 参数，它可以通过解释变量对异方差进行建模。 gnls() 的线性等效是广义最小二乘法，即 运行 与 nlme.
的 gls() 函数
现在我想计算 R 中的置信区间，并将它们与我的模型一起绘制在 ggplot() 中（ggplot2 包）。我为 gls() 对象执行此操作的方式是：
NewData <- data.frame(Explanatory1 = c(...), Explanatory2 = c(...))
NewData$fit <- predict(model, newdata = NewData)
到目前为止，一切正常，我的模型也很合适。
modmat <- model.matrix(formula(model)[-2], NewData)
int <- diag(modmat %*% vcov(model) %*% t(modmat))
NewData$lo <- with(NewData, fit - 1.96*sqrt(int))
NewData$hi <- with(NewData, fit + 1.96*sqrt(int))
这部分不适用于 gnls()，所以我无法获得我的上限和下限模型预测。
由于这似乎不适用于 gnls() 对象，我查阅了教科书以及以前提出的问题，但 none 似乎符合我的需要。我发现的唯一类似问题是 。在最佳答案中，建议使用 investr::predFit() 或使用 drc::drm() 构建模型，然后使用常规 predict() 函数。 None 这些解决方案帮助我 gnls()。
我目前的最佳解决方案是使用 confint() 函数计算所有三个参数（C、k、A）的 95% 置信区间，然后为置信上限和置信下限编写两个单独的函数，即一种使用 Cmin、kmin 和 Amin，一种使用 Cmax、kmax 和 Amax。然后我使用这些函数来预测值，然后用 ggplot() 绘制这些值。但是，我对结果并不完全满意，不确定这种方法是否最优。
这是一个最小的可重现示例，为简单起见，忽略了第二个分类解释变量：
# generate data
set.seed(10)
x <- rep(1:100,2)
r <- rnorm(x, mean = 10, sd = sqrt(x^-1.3))
y <- exp(-0.05*x) + r
df <- data.frame(x = x, y = y)
# find starting values
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 9.98071, C = 10.85413, logk = -3.14108
plot(m) # clear heteroskedasticity
# fit generalised nonlinear least squares
require(nlme)
mgnls <- gnls(y ~ C * exp(k * x) + A,
start = list(C = 10.85413, k = -exp(-3.14108), A = 9.98071),
weights = varExp(),
data = df)
plot(mgnls) # more homogenous
# plot predicted values
df$fit <- predict(mgnls)
require(ggplot2)
ggplot(df) +
geom_point(aes(x, y)) +
geom_line(aes(x, fit)) +
theme_minimal()
根据 Ben Bolker 的回答进行编辑
应用于第二个模拟数据集的标准非参数自举解决方案，它更接近我的原始数据并包含第二个分类解释变量：
# generate data
set.seed(2)
x <- rep(sample(1:100, 9), 12)
set.seed(15)
r <- rnorm(x, mean = 0, sd = 200*x^-0.8)
y <- c(200, 300) * exp(c(-0.08, -0.05)*x) + c(120, 100) + r
df <- data.frame(x = x, y = y,
group = rep(letters[1:2], length.out = length(x)))
# find starting values
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 108.9860, C = 356.6851, k = -2.9356
plot(m) # clear heteroskedasticity
# fit generalised nonlinear least squares
require(nlme)
mgnls <- gnls(y ~ C * exp(k * x) + A,
start = list(C = c(356.6851,356.6851),
k = c(-exp(-2.9356),-exp(-2.9356)),
A = c(108.9860,108.9860)),
params = list(C ~ group, k ~ group, A ~ group),
weights = varExp(),
data = df)
plot(mgnls) # more homogenous
# calculate predicted values
new <- data.frame(x = c(1:100, 1:100),
group = rep(letters[1:2], each = 100))
new$fit <- predict(mgnls, newdata = new)
# calculate bootstrap confidence intervals
bootfun <- function(newdata) {
start <- coef(mgnls)
dfboot <- df[sample(nrow(df), size = nrow(df), replace = TRUE),]
bootfit <- try(update(mgnls,
start = start,
data = dfboot),
silent = TRUE)
if(inherits(bootfit, "try-error")) return(rep(NA, nrow(newdata)))
predict(bootfit, newdata)
}
set.seed(10)
bmat <- replicate(500, bootfun(new))
new$lwr <- apply(bmat, 1, quantile, 0.025, na.rm = TRUE)
new$upr <- apply(bmat, 1, quantile, 0.975, na.rm = TRUE)
# plot data and predictions
require(ggplot2)
ggplot() +
geom_point(data = df, aes(x, y, colour = group)) +
geom_ribbon(data = new, aes(x = x, ymin = lwr, ymax = upr, fill = group),
alpha = 0.3) +
geom_line(data = new, aes(x, fit, colour = group)) +
theme_minimal()
这是生成的图，看起来很整洁！
我实施了 bootstrapping 解决方案。我最初做了标准的非参数 bootstrapping，它对 observations 进行了重新采样，但这会产生 95% 的置信区间，看起来范围很广——我 认为 这是因为 bootstrapping 的形式无法保持 x 分布的平衡（例如，通过重采样，您最终可能无法观察到较小的 x 值）。 （也有可能是我的代码中存在错误。）
作为第二个镜头，我切换到从初始拟合中重新采样 残差 并将它们添加到预测值；这是一个相当标准的方法，例如在 bootstrapping 时间序列中（尽管我忽略了残差自相关的可能性，这需要 block bootstrapping）。
这是基本的 bootstrap 重采样器。
df$res <- df$y-df$fit
bootfun <- function(newdata=df, perturb=0, boot_res=FALSE) {
start <- coef(mgnls)
## if we start exactly from the previously fitted coefficients we end
## up getting all-identical answers? Not sure what's going on here, but
## we can fix it by perturbing the starting conditions slightly
if (perturb>0) {
start <- start * runif(length(start), 1-perturb, 1+perturb)
}
if (!boot_res) {
## bootstrap raw data
dfboot <- df[sample(nrow(df),size=nrow(df), replace=TRUE),]
} else {
## bootstrap residuals
dfboot <- transform(df,
y=fit+sample(res, size=nrow(df), replace=TRUE))
}
bootfit <- try(update(mgnls,
start = start,
data=dfboot),
silent=TRUE)
if (inherits(bootfit, "try-error")) return(rep(NA,nrow(newdata)))
predict(bootfit,newdata=newdata)
}
set.seed(101)
bmat <- replicate(500,bootfun(perturb=0.1,boot_res=TRUE)) ## resample residuals
bmat2 <- replicate(500,bootfun(perturb=0.1,boot_res=FALSE)) ## resample observations
## construct envelopes (pointwise percentile bootstrap CIs)
df$lwr <- apply(bmat, 1, quantile, 0.025, na.rm=TRUE)
df$upr <- apply(bmat, 1, quantile, 0.975, na.rm=TRUE)
df$lwr2 <- apply(bmat2, 1, quantile, 0.025, na.rm=TRUE)
df$upr2 <- apply(bmat2, 1, quantile, 0.975, na.rm=TRUE)
现在画图：
ggplot(df, aes(x,y)) +
geom_point() +
geom_ribbon(aes(ymin=lwr, ymax=upr), colour=NA, alpha=0.3) +
geom_ribbon(aes(ymin=lwr2, ymax=upr2), fill="red", colour=NA, alpha=0.3) +
geom_line(aes(y=fit)) +
theme_minimal()
pink/light-red区域是观察级bootstrapCI（可疑）；灰色区域是残差 bootstrap CIs.
最好也尝试 delta 方法，但是 (1) 它比 bootstrapping 更强大 assumptions/approximations 并且 (2) 我没时间了。
我用 R 包 nlme 和包含的 gnls() 函数构建了几个广义非线性最小二乘模型（指数衰减）。我不简单地使用基础 nls() 函数构建非线性最小二乘模型的原因是因为我希望能够对异方差建模以避免转换。我的模型看起来像这样：
model <- gnls(Response ~ C * exp(k * Explanatory1) + A,
start = list(C = c(C1,C1), k = c(k1,k1), A = c(A1,A1)),
params = list(C ~ Explanatory2, k ~ Explanatory2,
A ~ Explanatory2),
weights = varPower(),
data = Data)
与简单 nls() 模型的主要区别在于 weights 参数，它可以通过解释变量对异方差进行建模。 gnls() 的线性等效是广义最小二乘法，即 运行 与 nlme.
gls() 函数
现在我想计算 R 中的置信区间，并将它们与我的模型一起绘制在 ggplot() 中（ggplot2 包）。我为 gls() 对象执行此操作的方式是：
NewData <- data.frame(Explanatory1 = c(...), Explanatory2 = c(...))
NewData$fit <- predict(model, newdata = NewData)
到目前为止，一切正常，我的模型也很合适。
modmat <- model.matrix(formula(model)[-2], NewData)
int <- diag(modmat %*% vcov(model) %*% t(modmat))
NewData$lo <- with(NewData, fit - 1.96*sqrt(int))
NewData$hi <- with(NewData, fit + 1.96*sqrt(int))
这部分不适用于 gnls()，所以我无法获得我的上限和下限模型预测。
由于这似乎不适用于 gnls() 对象，我查阅了教科书以及以前提出的问题，但 none 似乎符合我的需要。我发现的唯一类似问题是 investr::predFit() 或使用 drc::drm() 构建模型，然后使用常规 predict() 函数。 None 这些解决方案帮助我 gnls()。
我目前的最佳解决方案是使用 confint() 函数计算所有三个参数（C、k、A）的 95% 置信区间，然后为置信上限和置信下限编写两个单独的函数，即一种使用 Cmin、kmin 和 Amin，一种使用 Cmax、kmax 和 Amax。然后我使用这些函数来预测值，然后用 ggplot() 绘制这些值。但是，我对结果并不完全满意，不确定这种方法是否最优。
这是一个最小的可重现示例，为简单起见，忽略了第二个分类解释变量：
# generate data
set.seed(10)
x <- rep(1:100,2)
r <- rnorm(x, mean = 10, sd = sqrt(x^-1.3))
y <- exp(-0.05*x) + r
df <- data.frame(x = x, y = y)
# find starting values
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 9.98071, C = 10.85413, logk = -3.14108
plot(m) # clear heteroskedasticity
# fit generalised nonlinear least squares
require(nlme)
mgnls <- gnls(y ~ C * exp(k * x) + A,
start = list(C = 10.85413, k = -exp(-3.14108), A = 9.98071),
weights = varExp(),
data = df)
plot(mgnls) # more homogenous
# plot predicted values
df$fit <- predict(mgnls)
require(ggplot2)
ggplot(df) +
geom_point(aes(x, y)) +
geom_line(aes(x, fit)) +
theme_minimal()
根据 Ben Bolker 的回答进行编辑
应用于第二个模拟数据集的标准非参数自举解决方案，它更接近我的原始数据并包含第二个分类解释变量：
# generate data
set.seed(2)
x <- rep(sample(1:100, 9), 12)
set.seed(15)
r <- rnorm(x, mean = 0, sd = 200*x^-0.8)
y <- c(200, 300) * exp(c(-0.08, -0.05)*x) + c(120, 100) + r
df <- data.frame(x = x, y = y,
group = rep(letters[1:2], length.out = length(x)))
# find starting values
m <- nls(y ~ SSasymp(x, A, C, logk))
summary(m) # A = 108.9860, C = 356.6851, k = -2.9356
plot(m) # clear heteroskedasticity
# fit generalised nonlinear least squares
require(nlme)
mgnls <- gnls(y ~ C * exp(k * x) + A,
start = list(C = c(356.6851,356.6851),
k = c(-exp(-2.9356),-exp(-2.9356)),
A = c(108.9860,108.9860)),
params = list(C ~ group, k ~ group, A ~ group),
weights = varExp(),
data = df)
plot(mgnls) # more homogenous
# calculate predicted values
new <- data.frame(x = c(1:100, 1:100),
group = rep(letters[1:2], each = 100))
new$fit <- predict(mgnls, newdata = new)
# calculate bootstrap confidence intervals
bootfun <- function(newdata) {
start <- coef(mgnls)
dfboot <- df[sample(nrow(df), size = nrow(df), replace = TRUE),]
bootfit <- try(update(mgnls,
start = start,
data = dfboot),
silent = TRUE)
if(inherits(bootfit, "try-error")) return(rep(NA, nrow(newdata)))
predict(bootfit, newdata)
}
set.seed(10)
bmat <- replicate(500, bootfun(new))
new$lwr <- apply(bmat, 1, quantile, 0.025, na.rm = TRUE)
new$upr <- apply(bmat, 1, quantile, 0.975, na.rm = TRUE)
# plot data and predictions
require(ggplot2)
ggplot() +
geom_point(data = df, aes(x, y, colour = group)) +
geom_ribbon(data = new, aes(x = x, ymin = lwr, ymax = upr, fill = group),
alpha = 0.3) +
geom_line(data = new, aes(x, fit, colour = group)) +
theme_minimal()
这是生成的图，看起来很整洁！
我实施了 bootstrapping 解决方案。我最初做了标准的非参数 bootstrapping，它对 observations 进行了重新采样，但这会产生 95% 的置信区间，看起来范围很广——我 认为 这是因为 bootstrapping 的形式无法保持 x 分布的平衡（例如，通过重采样，您最终可能无法观察到较小的 x 值）。 （也有可能是我的代码中存在错误。）
作为第二个镜头，我切换到从初始拟合中重新采样 残差 并将它们添加到预测值；这是一个相当标准的方法，例如在 bootstrapping 时间序列中（尽管我忽略了残差自相关的可能性，这需要 block bootstrapping）。
这是基本的 bootstrap 重采样器。
df$res <- df$y-df$fit
bootfun <- function(newdata=df, perturb=0, boot_res=FALSE) {
start <- coef(mgnls)
## if we start exactly from the previously fitted coefficients we end
## up getting all-identical answers? Not sure what's going on here, but
## we can fix it by perturbing the starting conditions slightly
if (perturb>0) {
start <- start * runif(length(start), 1-perturb, 1+perturb)
}
if (!boot_res) {
## bootstrap raw data
dfboot <- df[sample(nrow(df),size=nrow(df), replace=TRUE),]
} else {
## bootstrap residuals
dfboot <- transform(df,
y=fit+sample(res, size=nrow(df), replace=TRUE))
}
bootfit <- try(update(mgnls,
start = start,
data=dfboot),
silent=TRUE)
if (inherits(bootfit, "try-error")) return(rep(NA,nrow(newdata)))
predict(bootfit,newdata=newdata)
}
set.seed(101)
bmat <- replicate(500,bootfun(perturb=0.1,boot_res=TRUE)) ## resample residuals
bmat2 <- replicate(500,bootfun(perturb=0.1,boot_res=FALSE)) ## resample observations
## construct envelopes (pointwise percentile bootstrap CIs)
df$lwr <- apply(bmat, 1, quantile, 0.025, na.rm=TRUE)
df$upr <- apply(bmat, 1, quantile, 0.975, na.rm=TRUE)
df$lwr2 <- apply(bmat2, 1, quantile, 0.025, na.rm=TRUE)
df$upr2 <- apply(bmat2, 1, quantile, 0.975, na.rm=TRUE)
现在画图：
ggplot(df, aes(x,y)) +
geom_point() +
geom_ribbon(aes(ymin=lwr, ymax=upr), colour=NA, alpha=0.3) +
geom_ribbon(aes(ymin=lwr2, ymax=upr2), fill="red", colour=NA, alpha=0.3) +
geom_line(aes(y=fit)) +
theme_minimal()
pink/light-red区域是观察级bootstrapCI（可疑）；灰色区域是残差 bootstrap CIs.
最好也尝试 delta 方法，但是 (1) 它比 bootstrapping 更强大 assumptions/approximations 并且 (2) 我没时间了。