[R] ddply and nlminb

From: Robert Clement <robert.clement_at_ed.ac.uk>
Date: Wed, 13 Apr 2011 13:16:16 +0100


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.


______________________________________________ R-help_at_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.

Received on Wed 13 Apr 2011 - 14:34:43 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 13 Apr 2011 - 15:00:29 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive