This problem is not for the faint of heart. Doug Bates, author of
nls(...) has said that a general purpose implementaion of R code for
multiresponse nonlinear regression is unlikely in the near future. You
have a large set of issues to deal with here. First, you have a system
of differential equations that are nonlinear so they need to be solved
precisely. This is your starting point. Next you have three responses
- Arne's comment about minimizing the determinant of the covariance
model is a standard approach but not universal - it depends on the error
structure of the data. Third, optimization routines (there are several
in R) can be used but they are often persnickity (sp?) when used with
numerical solutions to the differential equations. I've used nlm and
optim both for multiresponse problems in R but with explicit solutions
rather than numerical methods. Fourth, there may be issues with
correlation among responses that can render the residual covariance
matrix (near) singular and the determinant can vanish. To get started I'd
1. Ask myself if simultaneous estimation is what you really want to do
or if you can get what you want from single response estimation using
nls(...). If yes, then
2. Read the book by Bates and Watts (there are others, but this one is
very concise and has examples). It's a John Wiley book called Nonlinear
regression analysis and its applications, 1988. It's very likely in the
UCSB library.
3. Start by estimating parameters for each response individually to see
if the model can fit the data. If it can't, you need to reformulate your
model and go back to 1. If so, then there's hope and proceed to trying
to use two responses, then three. Whenever you have multiple responses
you need to do the eigenvalue/eigenvector analysis described in that
book and several previous papers by George Box and colleagues.
4. Finally, if all goes well, calculate the Bayesian 95% joint
confidence regions for the parameter pairs to assess their uncertainty
and check the model residuals for compliance with normality,
independence, and constant variance.
5. Collect Nobel prize!
This is one approach I've used with success. There are others, as hinted
at by Arne. Good luck - keep us posted on how it goes.
Regards
David Stevens
On 11/15/2011 12:17 PM, Arne Henningsen wrote:
> Dear Louise
>
> On 15 November 2011 19:03, lstevenson
> wrote:
>> Hi all,
>>
>> I'm trying to estimate model parameters in R for a pretty simple system of
>> equations, but I'm having trouble. Here is the system of equations (all
>> derivatives):
>> eqAlgae<- (u_Amax * C_A) * (1 - (Q_Amin / Q_A))
>> eqQuota<- (p_max * R_V) / (K_p + R_V) - ((Q_A-Q_Amin)*u_Amax)
>> eqResource<- -C_A * (p_max * R_V) / (K_p + R_V)
>> eqSystem<- list(C_A = eqAlgae, Q_A = eqQuota, R_V = eqResource)
>>
>> I want to estimate u_Amax, Q_Amin, p_max and Q_Amin with the data I've
>> collected using least squares. I've tried using systemfit but I'm not sure
>> how to write out the equations (my attempt is above but that doesn't work
>> since I haven't given values to the parameters I'm trying to estimate -
>> should I give those parameters initial values?). I've looked into the other
>> functions to get least squares estimates (e.g. lm() ) but I'm not sure how
>> to use that for a system of equations. I have some experience with R but I'm
>> a novice when it comes to parameter estimation, so any help would be much
>> appreciated! Thank you!
> Your system of equations is non-linear in parameters. As lm() and
> systemfit() can only estimate models that are linear in parameters,
> you cannot use these commands to estimate your model. The "systemfit"
> package includes the function nlsystemfit() that is intended to
> estimate systems of non-linear equations. However, nlsystemfit() is
> still under development and often has convergence problems. Therefore,
> I wouldn't use it for "serious" applications. You can estimate your
> non-linear equations separately with nls(). If you want to estimate
> your equations jointly, I am afraid that you either have to switch to
> another software or have to implement the estimation yourself. You
> could, e.g., minimize the determinant of the residual covariance
> matrix with optim(), nlm(), nlminb(), or another optimizer or you
> could maximize the likelihood function of the FIML model using
> maxLik(). Sorry that I (and R) cannot present you a simple solution!
>
> Best wishes from Copenhagen,
> Arne
>
[[alternative HTML version deleted]]
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
