[R] sas retain statement in R or fitting differene equations in NLS

[Berend Hasselman]

On 08-03-2012, at 14:50, Journals wrote:

> I wish to fit a dynamical model in R and I am running in a problem that 
> requires some of your wisdom to solve. For SAS users I am searching for 
> the equivalent of the */retain/ *statement.
> 
> For people that want to read complicated explanations to help me:
> 
> I have a system of two equations written as difference equations here.
> 
> To boil it down. I have a dataframe with three variables y, X1, X2 which 
> are measured. Now, I want to estimate a model which says y= 
> yhat+epsillon as in standard nonlinear regression. Where yhat is the 
> predicted value of y.
> 
> yhat can be calculated as follows:
> 
> # This code illustrates how I could simulate the expected values of y if 
> I knew the values of the parameters tau and b.
> # but in reality I would like to estimate them.
> # code is for illustration of the principles and is not meant to be 
> functional!!
> yaux[1]<-0
> b<- a_number # b would have to become estimated by nls or nlme
> tau<- another_number # tau would also be estimated in nls or nlme
> for (t in 2:1000)
> {
> yaux[t+1] <- yaux[t] + (X1-yaux[t])/tau
> yhat[t+1] <- yaux[t+1]*X2[t+1]/(X2[t+1]+b)
> }
> 

This can't be correct. When t==2 you are referencing yaux[2] which has no value.
Secondly, you haven't indexed X1 in the equation for yaux[t+1]. Is it X1[t]?

Besides that you don't have to calculate yaux that way.
You can use filter. As in

# Initial value yaux is 0
yaux <- as.vector(filter(a*X1,c(1-a),method="recursive"))

where a = 1/tau and assuming you meant X1[t]. 

And then you don't need to calculate yhat recursively.
You can do it with

yhat <- yaux * X2/(X2+b)

> Now, my problem is that I do not know the values of /tau/ and /b/ and I 
> would like to estimated them by non-linear regression. This is easy in 
> the case of /b/ but for /tau /nls (or nlme etc) would have to remember 
> the value of /yaux /for the previous observation and I did not find any 
> syntactical mean to do that.
> 

You could do a non linear regression with

y.rhs <- function(a,b,X1,X2) as.vector(filter(a*X1,c(1-a),method="recursive")) * X2/(X2+b)
nls(y ~ y.rhs(a,b,X1,X2), start=list(a=a,b=b)) 

where tau=1/a

Here is an example with random data

N <- 500
a <- .4
b <- 1

# generate some random data
set.seed(1)
X1 <- 5*round(runif(N),2)
X2 <- round(runif(N),2)

yaux <- as.vector(filter(a*X1,c(1-a),method="recursive"))
yhat <- yaux * X2/(X2+b)
eps <- runif(N)/2
eps <- eps - mean(eps)
y <- yhat + eps 
sum((y-yhat)^2)

# starting values
a <- .6
b <- 2

y.rhs <- function(a,b,X1,X2) as.vector(filter(a*X1,c(1-a),method="recursive")) * X2/(X2+b)
nls(y ~ y.rhs(a,b,X1,X2), start=list(a=a,b=b)) 

Berend