CODEWITHSUNDEEP
r
Pinus2Quercus
Pinus2Quercus

Reputation: 91

Creating a fitted logarithmic model

I'm trying to plot out a flood frequency chart with some data I have. Here's the type of data I'm working with:

#Set up maximum flow data
flow=sample(seq(10,1000,20),100,replace=TRUE)
flow=as.data.frame(flow[order(flow, decreasing=TRUE)])
names(flow)="max"
#rank flows from largest to smallest
flow$"rank"=seq(1,nrow(flow),1)
#Calculate the return interval in years
flow$"RI"=(nrow(flow)+1)/flow$"rank"
plot(flow$"max"~flow$"RI",type='p',log='xy', xlab='Return Interval', ylab='Max flow')

Now we have the maximum recorded flow per year and an estimate of the recurrence interval. Now what I'd like to do is find the logarithmic line of best fit. I've been playing around a bit with the nls function but keep coming up with this error.

Error in parse(text = x) : <text>:2:0: unexpected end of input
1: ~ 
   ^

Here's a sample of what I'm doing with the nls function:

logMod = nls((flow$"max"~(a*log10(flow$"RI")+b)),start = list(a = 0, b = 0))

Can someone help me out and let me know where I've gone astray?

Upvotes: 3

Views: 92

Answers (1)

MrFlick
MrFlick

Reputation: 206242

You should not use $ when specifying data from a data.frame in a nls formula. You should use the data= parameter to specify which data.frame to use to look up variable values. Thus you should change your call to 

logMod = nls( max ~ (a*log10(RI)+b), data=flow, 
    start = list(a = 0, b = 0))

The problem seems to be from using a$"b" rather than a$b as is more common. nls() extracts variable names using all.vars(). And observe

all.vars(flow$"max" ~ (a * log10(flow$"RI") + b))
# [1] "flow" "a"    "b"   
all.vars(flow$max ~ (a * log10(flow$RI) + b))
# [1] "flow" "max"  "a"    "RI"   "b"

This is because when you use the quotes, you no longer are specifying your columns as symbol/names for all.vars() to find, but rather you are passing them as character value which are not extracted. So in this case a$b is slightly different than a$"b"

Also, as @Gregor pointed out, if you're just just doing a log transform on one of your predictors, this is basically still a linear model. You can do

lm( max ~ log10(RI), data=flow)

Upvotes: 2

Related Questions