R help - Dec 2009 - need a solution to an R-problem: consultant available?

I am trying to get confidence bands for a non-linear power function
(y=mx^b).   I thought I should be able to figure it out, but can't.  Are
there any R consultant?  I would be willing to pay some amount of money, but
not sure such consultants exist.  I fit power functions to lots of data, and
this would be very useful.  I would ideally like to have confidence bands
for the mean function and a single prediction.  

I have posted my current Rscript below, with the location of the error
noted, in addition, below the Rscript, is some example data. 

Let me know what you think.

Cheers,

Andre

###RSCIPT STARTS

 

require(graphics)

rm(list=ls()) # clears everything in working directory

data = read.delim("F:/R_code/exampledata.txt")

#head(data) # shows header

 

attach(data) #attaches data to current run

 

datasort=data[order(Stage), ] #sorts data so predicted values follow in
systematic order

 

#datasort # shows sorted data

 

starts = list(const =0.00001, alpha = 20)

 

fm <- nls(Discharge ~ const*Stage ^alpha, start=starts) #runs a non linear
least squares regression model 

 

#predict(fm)              # fitted values at observed times

 

#summary(fm)

coef=coefficients(fm)

const_val=round(coef["const"],4) #makes a variable for the equation in
the
plot and rounds it to 4 decimal places

exp_val=round(coef["alpha"],4)

plot(Discharge~Stage, 

                #xlim=c(min(Stage),max(Stage)), #sets limits of x and y axis
based on max and min in data set

                #ylim=c(min(Discharge),max(Discharge)),

                xlab="Stage in East Creek", #label of X axis

                ylab="Discharge in East Creek Plunge Pool",  # label
of y
axis

                #xaxs="i",  #fits axis to exact tick mark location,
not
standard + 4 %

                #yaxs="i", 

                log = "xy",

                pch=15

)

#adds a legend, exponent should really be up higher, but not sure how to do
this.

legend("bottomright",  ,bty="n",
substitute(paste("Discharge", " = ",
const_val, , "Stage^",exp_val),

 list(const_val = const_val, exp_val=exp_val)))

 

tt <- seq(min(Stage),max(Stage), length = 101)

 

points(Stage,fitted.values(fm),type="l")

 

se.fit <- sqrt(apply(attr(predict(fm,list(Stage =
tt)),"gradient"),1,

                  function(x) sum(vcov(fm)*outer(x,x))))

###Above is where the error arrives. dim(X) must have a positive length.
"gradient' data does not get created with the predict function, which
is
likely the source of the problem 

matplot(tt, predict(fm,list(Stage = tt))+

               outer(se.fit,qnorm(c(.5, .025,.975))),type="l")  

###example data below


Discharge

Stage


0.069

1.9


0.099

1.93


0.001

1.74


0.001

1.76


0.005

1.71


1.558

2.77


0.44

2.26


0.467

2.25


0.086

1.84


0.012

1.86


0.045

1.93


0.043

1.93


0.05

1.91


0.016

1.84


0.142

2.07


0.656

2.37


0.576

2.33


0.329

2.22


0.321

2.21


0.378

2.24


0.375

2.24


0.132

2.07


0.125

2.07


0.465

2.27

 

 


	[[alternative HTML version deleted]]