logo标签列表

在R中，非线性回归(nls)存在的问题

rnls
3

我正在尝试在R中解决一个非线性回归问题，但是一直收到语法错误的提示信息。我已经尝试过调试，并且也招募了同事帮忙，但都没有成功。由于我刚接触R，我不确定是否正确设置了程序。

如果有任何建议或见解，将不胜感激。

谢谢。

## read in files from snotel sites to develop parameters
BASE1 <-read.csv("base_data.csv")

## work with the base file first
attach(BASE1)

## set parameters
S1 <- seq(-10,10,0.5)
S2 <- seq(0.1,2.5,0.2)
M1 <- seq(-10,10,0.5)
M2 <- seq(0.1,2.5,0.2)

## define the sub-functions used in the model
S_EXP <- (exp(-(PRISM_T + S1)/S2))
M_EXP <- (exp(-(PRISM_T + M1)/M2))
C_SNOW <- (1 - (1/(1 + S_EXP)))
C_MELT <-(1/(1 + M_EXP))

swe.mod <- nls(SNOTEL_SWE ~ SURPLUS * C_SNOW + PREV_SWE * (1 - C_MELT),data=BASE1,
+ start=list(S1 = -10, S2 = 0.1, M1 = -10, M2 = 0.1) , trace=TRUE)

玩具数据集：

BASE1<-structure(list(OBS = 1:61, SURPLUS = c(55.9, 124.5, 138.4, 124.1, 107.8, 102.9, 84.8, 74.4, 40.1, -23.1, -23.1, 20.7, 73.8, 267.9, 282.6, 244.2, 234.9, 199, 118.9, 55, -8.1, -59.4, -51.7, -10.4, 23.6, 111, 111.7, 107.7, 88.2, 88.1, 54.1, 18.8, -17.1, -65.4, -54.1, -16, 62.6, 178.5, 201, 192.8, 170.5, 173.5, 85.6, 30.4, -15.3, -69.7, -50.5, -4.8, 90, 243.2, 247.7, 234.6, 213.7, 194, 105.5, 16.6, -21.9, -73.8, -62.7, -6.8, 74.8), PRISM_T = c(1.9, -3.6, -6.3, -7, -5.7, -4.2, -1.3, 2.4, 5.8, 10.2, 10.4, 6.5, 5.9, 0.4, -2.3, -2.9, -2.1, -0.7, 1.4, 4.7, 8.6, 13.2, 13.3, 10.4, 5, -0.2, -3, -4.3, -2, 0, 2.1, 5.8, 9.8, 14.6, 14.4, 10.9, 7.1, 1.7, -0.3, -1.2, -0.6, 0.4, 2.2, 6.1, 9.7, 14.8, 14.8, 12, 8.5, 3, 1.3, 0.3, 1.3, 1.5, 2.8, 6.6, 10.6, 15.2, 15.7, 13.3, 7), SNOTEL_SWE = c(7.62, 50.8, 180.34, 317.5, 434.34, 562.61, 660.4, 622.3, 306.07, 36.83, 0, 0, 2.54, 49.53, 241.3, 488.95, 711.2, 895.35, 1093.47, 957.58, 372.11, 0, 0, 0, 1.27, 25.4, 137.16, 256.54, 379.73, 501.65, 549.91, 287.02, 8.89, 0, 0, 0, 2.54, 27.94, 177.8, 318.77, 459.74, 612.14, 730.25, 584.2, 142.24, 0, 0, 0, 0, 6.35, 91.44, 167.64, 256.54, 330.2, 267.97, 129.54, 0, 0, 0, 0, 10.16), PREV_SWE = c(0, 7.62, 50.8, 180.34, 317.5, 434.34, 562.61, 660.4, 622.3, 306.07, 36.83, 0, 0, 2.54, 49.53, 241.3, 488.95, 711.2, 895.35, 1093.47, 957.58, 372.11, 0, 0, 0, 1.27, 25.4, 137.16, 256.54, 379.73, 501.65, 549.91, 287.02, 8.89, 0, 0, 0, 2.54, 27.94, 177.8, 318.77, 459.74, 612.14, 730.25, 584.2, 142.24, 0, 0, 0, 0, 6.35, 91.44, 167.64, 256.54, 330.2, 267.97, 129.54, 0, 0, 0, 0)), .Names = c("OBS", "SURPLUS", "PRISM_T", "SNOTEL_SWE", "PREV_SWE"), class = "data.frame", row.names = c(NA, -61L))
4
请提供用 dput() 函数生成的玩具数据集，以便我们可以运行代码，而不是使用 BASE1 <- read.csv("base_data.csv") 读取真实数据集。注意不要改变原始数据集的含义或信息。 - Michael
我们需要查看BASE1中保存的一些数据细节，请包括?str、?head和/或?dput。（PS：attach以搞砸事情而闻名，因此请使用nls的数据参数代替） - Brandon Bertelsen
dput(toy_data_set)将返回一个结果，我们可以将其读入R中。 - Michael
1个回答
3

我知道你想做什么，但是有点难理解。你不能像这样定义“子函数”。你必须在模型中明确指定它们。例如，你的公式应该是:

swe.mod <- 
nls(
  SNOTEL_SWE ~ SURPLUS * (1-(1/(1 + (exp(-(PRISM_T + S1)/S2)))))
  + PREV_SWE * (1 - 1/(1 + (exp(-(PRISM_T + M1)/M2)))),
data=BASE1,
start=list(S1 = -10, S2 = 0.1, M1 = -10, M2 = 0.1))

如果您觉得这样有点丑，可以使用函数来定义您的“子函数”，如 ?nls 中所示的示例。请尝试以下操作：
S_EXP  <- function(PRISM_T,S1,S2) (exp(-(PRISM_T + S1)/S2))
C_SNOW <- function(PRISM_T,S1,S2) (1 - (1/(1 + S_EXP(PRISM_T,S1,S2) )))
M_EXP  <- function(PRISM_T,M1,M2) (exp(-(PRISM_T + M1)/M2))
C_MELT <- function(PRISM_T,M1,M2) (1/(1 + M_EXP(PRISM_T,M1,M2)))

swe.mod <- nls(SNOTEL_SWE ~ SURPLUS * C_SNOW(PRISM_T,S1,S2) 
+ PREV_SWE * (1 - C_MELT(PRISM_T,M1,M2)),
data=BASE1,start=list(S1 = -10, S2 = 0.1, M1 = -10, M2 = 0.1), trace=TRUE)

我无法理解您所说的是什么意思：

## set parameters
S1 <- seq(-10,10,0.5)
S2 <- seq(0.1,2.5,0.2)
M1 <- seq(-10,10,0.5)
M2 <- seq(0.1,2.5,0.2)

这些参数是您想要拟合的，对吗？您只需要指定初始参数，就像您所做的那样。此外，请避免使用attach。您不需要它，而且它会导致非常难以找到的错误。

以下是结果：

Formula: SNOTEL_SWE ~ SURPLUS * C_SNOW(PRISM_T, S1, S2) + PREV_SWE * (1 - 
    C_MELT(PRISM_T, M1, M2))

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
S1  -0.6546     0.4222  -1.550  0.12657    
S2   1.4182     0.4907   2.890  0.00543 ** 
M1  -7.4500     0.2335 -31.906  < 2e-16 ***
M2   1.6118     0.2046   7.879 1.09e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 51.01 on 57 degrees of freedom

Number of iterations to convergence: 20 
Achieved convergence tolerance: 6.617e-06
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接