r - R 中的非线性回归 (nls) 问题

Question

我正在尝试解决 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))

Answer 1

有点难以理解您要做什么，但我想我可以猜到。您不能像现在这样定义“子功能”。您必须在模型中明确指定它们。例如，您的公式应该是:

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