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.
