[R] profile
KATARINA,DOMIJAN
kd2 at waikato.ac.nz
Mon Feb 11 02:18:04 CET 2002
I am running 1.3.1 on a Windows (NT 4.0) machine. I've fit a nonlinear
model intended to predict crop yield from nutrient information, and want to
use the profile function. If I type say,
profile(simparj.fm)
I get the following error message:
"Error in prof$getProfile(): number of iterations exceeded maximum of
5.25515e-308"
I used the profiler function to profile simparj,fm step by step, profiling
one parameter at the time and calculating the profile t statistics. I've
played around with different values for the parameters and found that
profiler doesn't work for certain ranges of the parameter values.
So now I want to reduce the range of the profiling. I've been looking at the
help page for profile.nls and also I hacked the codes for profile.nls (just
by typing profile.nls), but haven't really found an efficient way to do it.
For say, my parameter eta1, I found that the profiler function works between
-1.6<tau<1.6. I adjusted the code in profile.nls to profile within that
range only (this is a very crude adjustment, see the code below). The new
profile works, and I can increase the number of points profiled (by changing
delta.t), but for some reason I cannot plot it. I get the following error
message:
Loading required package: splines
Error in interpSpline.default(obj[[i]]$par.vals[,i], obj[[i]]$tau
Thank you for your info
Katarina
My code in full is:
(Note: I'm just using simulated data in this example, and the profile will
sometimes work with this)
# parameter assignments
Nmin <- 0.8852
Nopt1 <- 16.78
gN <- 0.5511
E.nfert1 <- 0.3271
E.nfert2 <- 0.6132
Beta <- 0.8902
Dls <- 0.5378
eta1 <- 0.3791
eta2 <- 0.6332
PopStd <- 90468
beta <- Beta
DIs <- Dls
MnmN <- Nmin
OptN <- Nopt1
# define model function
Y.model <- function(gN, MnmN, OptN, DIs, beta, eta1, eta2,
Popn, Dmax, AWC, SumEp, PotYield3, Nsupply)
{
Ymax<- 1-ifelse(Popn<=PopStd, eta1, eta2)*log(Popn/PopStd)
Ymax <- Ymax*PotYield3*Popn/1000
Ymax <- Ymax*ifelse(Dmax<=DIs*AWC, 1, 1 - beta*(Dmax
-DIs*AWC)/SumEp)
Nstar <- (Nsupply- MnmN*Ymax) / (OptN*Ymax - MnmN*Ymax)
Nstar<-pmax ( 0,Nstar)
Ystar<-ifelse(Nstar<1, (1 + gN*(1 - Nstar))* Nstar^(gN), 1)
Ystar<-pmax ( 0, Ystar)
Y.model <- Ystar*Ymax
Y.model
}
# simulate experimental data for predictors
nsim <- 300
Popn <- rnorm(nsim,PopStd,0.1*PopStd)
Dmax <- rnorm(nsim,140.68,47.45)
AWC <- rnorm(nsim,186.86,47.41)
SumEp <- rnorm(nsim,318.54,32.53)
PotYield3 <- rnorm(nsim,0.16180,0.01167)
Nsoil <- rnorm(nsim,94.07,34.06)
Bdfield <- rnorm(nsim,1.0590,0.1420)
Bdlab <- rnorm(nsim,0.7876,0.1169)
Nfert.broad <- runif(nsim,95.3,576.5)
Nfert.band <- runif(nsim,122,250)
broad <- rbinom(nsim,1,0.2)
Nfert.broad <- Nfert.broad*broad
Nfert.band <- Nfert.band*(1 - broad)
Nsupply<-Nsoil*Bdfield/Bdlab + Nfert.broad*E.nfert1 +
Nfert.band*E.nfert2
# generate response variable from model
Y.m <- Y.model(gN, MnmN, OptN, DIs, beta, eta1, eta2,
Popn, Dmax, AWC, SumEp, PotYield3, Nsupply)
Y.simulated <- Y.m + 0.1*rnorm(nsim,0,1)
library(nls)
# attempt to fit starting from true parameters
simparj.st <- c(gN, MnmN, OptN, DIs, beta, eta1, eta2)
names(simparj.st) <- c('gN', 'MnmN', 'OptN', 'DIs', 'beta', 'eta1',
'eta2')
simparj.fm <- nls(Y.simulated ~ Y.model(gN, MnmN, OptN, DIs, beta,
eta1, eta2, Popn, Dmax, AWC, SumEp, PotYield3, Nsupply), start =
simparj.st, trace =T)
pr1<-profile(simparj.fm)
# define a new profile function
pr2<- function (fitted, which = 1:npar, maxpts = 100, alphamax = 0.01,
delta.t = cutoff/5)
{
f.summary <- summary(fitted)
std.err <- f.summary$parameters[, "Std. Error"]
npar <- length(std.err)
nobs <- length(resid(fitted))
cutoff <- sqrt(npar * qf(1 - alphamax, npar, nobs - npar))
outmat <- array(0, c(nobs, npar))
out <- list()
prof <- profiler(fitted)
on.exit(prof$setDefault())
for (par in which) {
pars <- prof$getFittedPars()
prof$setDefault(varying = par)
sgn <- -1
count <- 1
varying <- rep(TRUE, npar)
varying[par] <- FALSE
tau <- double(2 * maxpts)
par.vals <- array(0, c(2 * maxpts, npar), list(NULL,
names(pars)))
tau[1] <- 0
par.vals[1, ] <- pars
base <- pars[par]
profile.par.inc <- delta.t * std.err[par]
pars[par] <- pars[par] - profile.par.inc
while (count <= maxpts) {
if (abs(pars[par] - base)/std.err[par] > 10 * cutoff)
break
count <- count + 1
prof$setDefault(params = pars)
ans <- prof$getProfile()
tau[count] <- sgn * sqrt(ans$fstat)
par.vals[count, ] <- pars <- ans$parameters
if (abs(tau[count]) > 1.6) #my adjustment
break
pars <- pars + ((pars - par.vals[count - 1, ]) *
delta.t)/abs(tau[count] - tau[count - 1])
}
tau[1:count] <- tau[count:1]
par.vals[1:count, ] <- par.vals[count:1, ]
sgn <- 1
newmax <- count + maxpts
pars <- par.vals[count, ]
while (count <= newmax) {
pars <- pars + ((pars - par.vals[count - 1, ]) *
delta.t)/abs(tau[count] - tau[count - 1])
if (abs(pars[par] - base)/std.err[par] > 10 * cutoff)
break
count <- count + 1
prof$setDefault(params = pars)
ans <- prof$getProfile()
tau[count] <- sgn * sqrt(ans$fstat)
par.vals[count, ] <- pars <- ans$parameters
if (abs(tau[count]) > 1.6) # my adjustment
break
}
out[[par]] <- structure(list(tau = tau[1:count], par.vals =
par.vals[1:count,
]), class = "data.frame", row.names = as.character(1:count),
parameters = list(par = par, std.err = std.err[par]))
prof$setDefault()
}
names(out)[which] <- names(coef(fitted))[which]
attr(out, "original.fit") <- fitted
attr(out, "summary") <- f.summary
class(out) <- c("profile.nls", "profile")
out
}
pr1<-pr2(simparj.fm, which=5)
# plot pr1
opar<-par(mfrow=c(1,1), oma=c(1.1,0,1.1,0), las=1)
plot(pr1, conf=c(95,90,80,50)/100)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list