using non-linear model (nlm) as smooth as star_smooth function

[Peter Tittmann]

All,

I have a dataset that i'm using to create a nonlinear least-squares regression. I'd like to use the stat_smooth gemoetry to represent the line and confidence interval. The formula for the non-linear model is: dependent_variable =  a * independent_variable^b. 

I'm using the sort of lengthy approach proposed here http://perfdynamics.blogspot.com/2010/09/confidence-bands-for-universal.html

but wondered if anyone had suggestions for something simpler?

Thanks!

Peter

-- 
Peter Tittmann

[Dennis Murphy]

Hi:

One can fit your model in nls() as follows, per the following toy example:

x <- 1:20
y <- 2 * sqrt(x) + rnorm(20, s = 0.3) 
dd <- data.frame(x, y)
rm(x, y)

# True f = 2 * sqrt(x):
m <- nls(y ~ 0 + A * x^B, data = dd, start = list(A = 1, B = 1))
> summary(m)

Formula: y ~ 0 + A * x^B

Parameters:
  Estimate Std. Error t value Pr(>|t|)   
A  1.99163    0.12631   15.77 5.57e-12 ***
B  0.50860    0.02465   20.63 5.62e-14 ***

Residual standard error: 0.355 on 18 degrees of freedom

Number of iterations to convergence: 5
Achieved convergence tolerance: 1.010e-06

predict.nls() will produce predicted values from the fitted model, but does not generate standard errors of the predictions nor does it produce confidence limits, as documented on its help page.

Here's one way to get them and then produce a plot in ggplot2, but the key is understanding that the standard error of the predictions depends upon the derivative matrix of the model, which means that, unlike linear models, the standard errors of predictions are *model-dependent*, which is why they cannot be 'easily' implemented in a structured graphics engine like ggplot2. [Your reference blogger could probably stand to learn this little fact as well, given the carping at the end of his post.] You can probably work around this a little bit by generating your own selfStart() function in nls, but I'm not that sophisticated yet, so I'll just brute force it for illustration.

# Plotting confidence bands for a 'simple' nonlinear model in ggplot2
# Add predictions and margin of error (moe) to the original data frame dd

dd$yhat <- predict(m)   # the easy part...

# Trick is to get the margin of error for the CIs...

# In this example, f(x; A, B) = A * x^B
# The partial derivatives of f with respect to the parameters are
#  df/dA = x^B,  df/dB = A * x^B * log(x)
# The derivative matrix required has the derivatives as column functions,
# => each row consists of the two derivatives evaluated at the x value
# and the estimated values of A and B

# Estimated derivative matrix at the LS estimates (note that u = a * x^b,
# where a and b are the LSEs) [needed for the estimated SE below]

v0 <- with(dd, cbind(yhat/coef(m)[1], yhat * log(x)))

# Estimated vector of margins of error for the estimated expectation
# function f-hat: 

dd$moe <- qt(.975, 18) * sqrt(diag(v0 %*% vcov(m) %*% t(v0)))

# Now do the ggplot:

g <- ggplot(dd)
g + geom_point(aes(x = x, y = y), size = 2.5) +
    geom_ribbon(aes(x = x, y = yhat, ymin = yhat - moe, ymax = yhat + moe),
                fill = alpha('gray20', 0.1)) +
    geom_line(aes(x = x, y = yhat), color = 'red', size = 2) +
    geom_line(aes(x = x, y = yhat - moe), color = 'blue', linetype = 'dashed', size = 1) +
    geom_line(aes(x = x, y = yhat + moe), color = 'blue', linetype = 'dashed', size = 1)

It looks about right to me, but please double check. Not the most gorgeous graph I've ever done, but...

It may be possible to write this up into a geom, but one would need the gradient of the nonlinear function and, more likely, a selfStart() function to feed to nls() to get something useful out of it.

Hope this helps to explain why there is little about nonlinear regression with mechanistic models that is 'simple', including graphics.

Dennis

but wondered if anyone had suggestions for something simpler.

Thanks!

Peter

-- 
Peter Tittmann

--
You received this message because you are subscribed to the ggplot2 mailing list.
Please provide a reproducible example: http://gist.github.com/270442
 
To post: email ggp...@googlegroups.com
To unsubscribe: email ggplot2+u...@googlegroups.com
More options: http://groups.google.com/group/ggplot2

[Peter Tittmann]

Hi,

Thanks very much for taking the time on this. I realize that this is sort of a pure R question but I'll ask anyway:

can you tell me what exactly the %*% operator does to be able to multiply the vcov() and v0 matrices since the vcov is a 2x2 and the v0 is a 2x20?  Its sort of a weird question but Im trying to replicate the behavior in python. 

Thanks again!

-- 
Peter Tittmann

Attachments:

- nlfit.png

[Ben Bolker]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 10-12-01 04:08 PM, Peter Tittmann wrote:
>
>
> Hi,
>
> Thanks very much for taking the time on this. I realize that this is
> sort of a pure R question but I'll ask anyway:
>
> can you tell me what exactly the %*% operator does to be able to
> multiply the vcov() and v0 matrices since the vcov is a 2x2 and the v0
> is a 2x20? Its sort of a weird question but Im trying to replicate the
> behavior in python.
>
> Thanks again!

what's the problem? %*% is standard matrix multiplication --
multiplying an (m,n) by an (n,p) matrix should Just Work ...

(hope I'm not missing something obvious myself)

Ben Bolker

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.10 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iEYEARECAAYFAkz2uqMACgkQc5UpGjwzenPazwCfcUqC2gDoqIyvbVJew7V0y9F4
HYQAn1xDibYNAQEj+Y15QYxJoKcT0aOX
=VgHg
-----END PGP SIGNATURE-----

[PeterT]

You are absolutely right. It was my mistake, it was a problem with my
python data structures not letting me do matrix multiplication.

Thanks,

Peter

On Dec 1, 1:14 pm, Ben Bolker <bbol...@gmail.com> wrote:
> -----BEGIN PGP SIGNED MESSAGE-----
> Hash: SHA1
>
> On 10-12-01 04:08 PM, Peter Tittmann wrote:
>
>
>
> > Hi,
>
> > Thanks very much for taking the time on this. I realize that this is
> > sort of a pure R question but I'll ask anyway:
>
> > can you tell me what exactly the %*% operator does to be able to
> > multiply the vcov() and v0 matrices since the vcov is a 2x2 and the v0
> > is a 2x20?  Its sort of a weird question but Im trying to replicate the
> > behavior in python.
>
> > Thanks again!
>
>   what's the problem? %*% is standard matrix multiplication --
> multiplying an (m,n) by an (n,p) matrix should Just Work ...
>
>   (hope I'm not missing something obvious myself)
>
>   Ben Bolker
>
> -----BEGIN PGP SIGNATURE-----
> Version: GnuPG v1.4.10 (GNU/Linux)

> Comment: Using GnuPG with Mozilla -http://enigmail.mozdev.org/

[Peter Tittmann]

Dennis and Ben,

Thanks again for your assistance. I've successfully generated the plot using ggplot2 with python via rpy2. Your effort is very much appreciated. I've included the fruits of this labor which represents the relationship between diameter and height for +/- 1300 trees. If anyone is really interested in the python code they can contact me, the essence is in Dennis' response to my initial question..

gratefully,

-- 
Peter Tittmann

On Wednesday, December 1, 2010 at 2:41 PM, Dennis Murphy wrote:

Hi:

Ben's got it right - it's simply a matrix multiplication operator. In the example, the multiplication produces a 20 x 20 matrix, the only relevant part of which is the diagonal.

Essentially, the formula works on an individual basis (Bates and Watts, 1988, section 2.3.2, pp. 58-59). The term under the square root is a quadratic form v'Av, where v is a row of the derivative matrix evaluated at x and the LSEs and A is the covariance matrix of the estimated coefficients. All I did was to try to 'vectorize' the operation by computing the covariance matrix of all the y's and then just pull out the diagonals (the variances).

Ref:  Bates and Watts (1988). Nonlinear Regression Analysis and Its Applications. NY: Wiley.

HTH,
Dennis

[Peter Tittmann]

I used the following, BTW this is python with Rpy (http://rpy2.bitbucket.org/)

`ggplot2.theme_bw()+ggplot2.scale_x_continuous('height (m)')+ggplot2.scale_y_continuous('DBH (cm)’)`

Peter

On Aug 30, 2016, at 5:45 AM, Mrunmayee Sahoo <mrunmay...@gmail.com> wrote:

hi

the code helps me a lot. can you please tell me how to add legend for the lines and dots here.

thank you in advance

[Mrunmayee Sahoo]

hi

the code helps me a lot. can you please tell me how to add legend for the lines and dots here.

thank you in advance

On Wednesday, December 1, 2010 at 6:27:31 AM UTC+5:30, Dennis Murphy wrote:

[aidepl...@gmail.com]

You are has the search for financing for either revive your activities either for the realisation of a project, either to buy you a apartment but alas the bank poses to you has conditions of which you are unable to fill. Most concern me I am a particular I provides loans ranging from 1000 to 2000 000 dollars has all people able to meet its commitments.moreover, the interest rate is 2% per year.

Either you need money for other reasons ; do not hesitate to contact me for more information. Mail:aidepl...@gmail.com

[cab...@dse.ufpb.br]

Hi Denis,

I have carefully read your thread bellow. It has been some time since you posted this response, so it is copied bellow. I am sorry for resurrecting this topic from the dead, but this was the best response I found within days of search.

I wish to build a graph from a non-linear model through nils function. It seems not to be a way of doing that in nls, so I got here, for I wish to include error in the graph.

My data is composed of two variables (x, y) which I measured in the field, and I wish to use to estimate two parameters (A, B) based on the model y ~ A - 20*log10(x) - B*x

TO make it easier, I put an example of data and the nls model: 

#Input data (measured variables)

d <- matrix(c(70,10,64,20,58,40,52,79), ncol=2)

colnames(d) <- c('x', 'y')

trial.table <- as.table(d)

#model

m <- nls(y ~ A - 20*log10(x) - B*x, data = dd, start = list(A = 1, B = 1))

As I understood from you response, to get to the SE you have to:

# Estimate the derivative matrix at the LS estimates (note that u = a * x^b,
# where a and b are the LSEs). 

However, at this point I am not sure how to construct the code bellow from function f(x) partial derivatives.
#v0 <- with(dd, cbind(yhat/coef(m)[1], yhat * log(x)))

From the model, we can easily calculate the partial derivatives as following df/dA = 1 and df/dB = x. 

And as I understood (and I am sorry if I lost something), the SE calculations depends on the model, and we should use the partial derivatives to calculate v0. 

But how do I do that? 

Best wishes,

Carlos.

PS: This is the first time I write in this forum, so I sorry if the message is somehow unclear or out of the format.