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<louise.stevenson@lifesci.ucsb.edu>  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]]

