CODEWITHSUNDEEP
rnls
user1553041
user1553041

Reputation: 63

Problems with nonlinear regression (nls) in R

I am trying to solve a non-linear regression in R and am getting syntax error messages. I have tried to debug and also recruited a co-worker to no avail. I am not sure that the procedure is set up properly as I am new to R.

Any insight would be great-

Thanks

## 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)

toy_data_set:

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))

Upvotes: 3

Views: 929

Answers (1)

nograpes
nograpes

Reputation: 18323

It was a little hard to understand what you were trying to do, but I think I can guess. You can't define the "subfunctions" as you are doing. You must specify them explicitly in the model. For example, your formula should be:

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))

If you find that a little ugly, you can use functions to define your "subfunctions" as shown in the examples in ?nls. Try this:

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)

What I couldn't figure out is what you meant by:

## 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)

Those parameters are what you are trying to fit, right? You only need to specify initial parameters, as you did. Also, avoid attach. You don't need it, and it will lead to very hard to find bugs.

Here are the results:

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

Upvotes: 3

Related Questions