[R] ddply and nlminb

[Robert Clement]

Hello

I'm new to R (one week) so please excuse any obvious mistakes in my code 
or posting.

I am attempting to fit a non linear function defining the relationship 
between dependent variable A and  the variables PAR and T grouped by the 
condition Di.

The following steps are taken in the Rcode below:
1) load the data (not shown)
2) define the function to be fit
3) define the starting values of the fit parameters and their upper and 
lower limits
4) fit the function using all data using nlminb (Di groupings not 
considered)
5) estimate the parameter standard errors
6) fit the function with the data grouped by Di using ddply

The first 5 steps appear to be working alright and produce reasonable 
fit parameters and parameter errors.

However, when attempting to fit the function with the data grouped by Di 
the parameters are returned with the same value for each grouping:

 > R3Dpar
         Di         V1                       V2                   
V3             V4            V5          V6
1 0.033461 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
2 0.082682 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
3 0.133670 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
4 0.195940 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
5 0.272430 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
6 0.368960 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
7 0.500150 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518
8 0.748380 0.06081973 -2.119109e-05 31.75188 -2.631321 17.32565 0.821518


Can someone point out the error in my ways, and also the correct way to 
return both the parameter estimates and the parameter errors.

Thank you

Robert




# ------------ END OF CODE 
------------------------------------------------------------------------------------------


# ------------ Data loaded in data frame M2 and attached - not shown.

# ------------ Define function
fnonRecHypT <- function(x){
qeo   = x[1]
dqe   = x[2]
Fmo   = x[3]
TL    = x[4]
To    = x[5]
TH    = 45
theta = x[6]

qe = qeo + dqe*T
theta =
phi = (TH-To)/(To-TL)
Fm = Fmo*((T-TL)*(TH-T)^ phi) / ((To-TL)*(TH-To)^ phi)
Aest = (qe*PAR + Fm - ((qe*PAR + Fm)^2 - 4*theta*qe*PAR*Fm)^0.5)/(2.*theta)
result = sum((Aest-A)^2 )
}

# ------------ Define parameter starting values and limits
startval =  c(0.05,-0.01,20, -5,15,0.5)
lowval   =  c(0.01,-0.05, 5,-15, 7,0.1)
uppval   =  c(0.2,  0.02,50,  0,25,0.99)

# ------------ Fit using entire data set
R3 <-  nlminb(startval, fnonRecHypT, control=list(trace=1), 
lower=lowval, upper = uppval)

# ------------ estimate fit parameter standard errors
x = R3$par
D3 = hessian(fnonRecHypT,x)
e3 = sqrt(diag(solve(D3)))

# ------------ attempt to fit by grouping variable, Di   (Di is a real 
vector with 8 possible values)
R3Dpar= ddply(M2,
             .(Di),
             function(x) {res <- nlminb(startval, fnonRecHypT, 
control=list(trace=1), lower=lowval, upper = uppval)
             return(res$par)
}
)


# ------------------------- END OF CODE 
-----------------------------------------------------------------------------

-- 
School of Geosciences
University of Edinburgh
202 Crew Building
West Mains Road
Edinburgh
Scotland
EH9 3JN

Ph  +44 (0)131 650 7732
Fx  +44 (0)131 662 0478
email  robert.clement at ed.ac.uk


The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.