在nls中使用“predict”

Question

我从USGS Nation Water Data网站获取了数据。目前我正在尝试绘制并适应该数据的曲线，以用于预测数据集中不同测量值（溶解氧、pH、高程和温度）与流量的关系。我使用了“nls”命令，并使用了一本方程书来查找要使用的曲线……对于这个例子，我特别使用了Schumacher的方程（书上第48页）。

数据链接如下：

曲线书籍：http://www.for.gov.bc.ca/hfd/pubs/docs/bio/bio04.htm 我使用的数据：http://waterdata.usgs.gov/mi/nwis/uv?referred_module=qw&search_station_nm=River%20Rouge%20at%20Detroit%20MI&search_station_nm_match_type=anywhere&index_pmcode_00065=1&index_pmcode_00060=1&index_pmcode_00300=1&index_pmcode_00400=1&index_pmcode_00095=1&index_pmcode_00010=1&group_key=NONE&sitefile_output_format=html_table&column_name=agency_cd&column_name=site_no&column_name=station_nm&range_selection=date_range&begin_date=2013-11-18&end_date=2013-12-18&format=html_table&date_format=YYYY-MM-DD&rdb_compression=file&list_of_search_criteria=search_station_nm,realtime_parameter_selection

我的问题是，一旦我选择了一条曲线编码它，我就无法使用nls来预测新值...我也无法弄清楚如何绘制它...我猜这可以通过残差来解决吗？在代码中，我使用“聚合”提取列出测量值和相应的放电速率的平均值，现在我只需要让R为我预测即可。我已经得到了我认为是拟合值...但我不确定并被“?nls”挡住了。

##Create new dataframes with means given date for each constituent
ph <- aggregate(Discharge~pH, data=River.Data, mean)

##pH models
pH <- ph$pH
disch <- ph$Discharge
phm <- nls(disch~exp(a+(b/pH)), data=ph, trace=T, start=list(a=-47.06 ,b=400.2))
newph<- data.frame(ph=c(3.0,4.0,5.0,6.0,7.0,8.0,9.0))
predict(phm, newdata=newph)

- Jarrod

---

有点棘手：您必须在predict中使用与原始模型完全相同的变量名称。否则，就像您发现的那样，predict将返回拟合值而没有任何警告消息。我认为predict(phm,newdata=list(ph=newph))会起作用。 - Carl Witthoft

Carl，我确实尝试过这个方法...但它仍然返回拟合值。最让人沮丧的是，我无法理解为什么它只返回这些拟合值：predict(phm, newdata=newph) [1] 663.69857 460.76412 322.92607 228.39464 162.95840 117.25539 85.05862 62.18766 - Jarrod

Answer 1

似乎你已经得到了答案(??)，但是:

ph    <- aggregate(Discharge~pH, data=River.Data, mean)
phm   <- nls(Discharge~exp(a+(b/pH)), data=ph, trace=T, start=list(a=-47.06 ,b=400.2))
newph <- data.frame(pH=seq(3,9,by=0.1))
Discharge.pred <- predict(phm, newdata=newph)

plot(ph$pH, ph$Discharge, xlim=c(3,9), ylim=c(0,1000))
par(new=t)
plot(newph$pH,Discharge.pred, xlab="", ylab="", axes=F, xlim=c(3,9), ylim=c(0,1000), type="l")

问题是您的数据是在pH为[7.5,8.2]范围内，但您正在尝试预测pH为[3,9]。您选择的模型在范围外远离的pH值不稳定。

- jlhoward

---

是的，看起来是这样！谢谢您的回复，它似乎有效，但也许我需要选择一个更好的模型？ - Jarrod

Answer 2

Jarrod，试试这个。干杯，Robert。

#Try this
#pH <- ph$pH # you don't need this
#disch <- ph$Discharge # you don't need this
phm <- nls(Discharge~exp(a+(b/pH)), data=ph, trace=T, start=list(a=-47.06 ,b=400.2))
newph<- data.frame(pH=seq(3,9,0.1)) # it'll be smoother with a sequence in increments of 0.1
plot(newph,predict(phm, data=newph,type="l"))

- wobetfwoze

---

它好像也不喜欢那个... 它在绘图时出现了一个错误: > plot(newph,predict(phm, data=newph,type="l")

• plot(newph,predict(phm, data=newph,type="l") 错误: 在以下内容中出现了意外符号: "plot(newph,predict(phm, data=newph,type="l") plot"

- Jarrod

如果你将这个函数从plot中取出并仅使用predict，它仍然只会返回拟合值...也许这是模型的错误？我应该使用“glm”吗？ - Jarrod

我忘记了在 'plot()' 语句中关闭圆括号。我已经对其进行了编辑，如上所示。 - wobetfwoze

Answer 3

这个答案基于@Carl Witthoft的评论。我将其突出显示为单独的答案，因为我只是浏览了一下评论，并没有真正理解Carl提出的建议。希望它能帮助其他也发现自己在这里的人。

从?predict.nls文件中可以看出，新数据需要是一个“命名列表或数据框”。这就是我过去失败的地方；列表或数据框的列需要有一个名称。正如Carl所说，这个名称需要完全与模型中使用的名称相同。我认为OP的问题可能是因为他们在模型中使用了pH，但在创建新数据时使用了ph。

- EcologyTom