# [R] regression questions

Paul, David A paulda at BATTELLE.ORG
Wed Sep 10 23:11:57 CEST 2003

```I have been puzzling over how to fit some fixed effects models
to a set of data.  My response function is

response <- function(a, b, c, alpha1, alpha2, indicator, t, t2)
{
z = a +
b * (t) * exp(-alpha1 * t) +
indicator *c * (t2) * exp(-alpha2 * t2)
}

where t2 = t - 4 and "indicator" is a 0-1 vector denoting
when t > 4.  Each test subject receives equal doses at t = 0 and
t = 4.  The dose can vary from subject to subject.

Also note the following:
1.  Var(e[it]) = sigma1^2 for t<=4; Var(e[it]) = sigma2^2 for t>4.
This is motivated by my data exploration.
2.  b,c > 0 for biological interpretability
3.  t varies over {0,2,4,6,8,10}.
4.  For a variety of reasons, "a", "alpha1", and "alpha2" must
be held constant over all of the test subjects.

The function nlsList( ) is not appropriate because it assumes
that all of the parameters are allowed to vary with each level of
a specified grouping variable (in this case, "subject.id").
I have been able to fit nls( ) models using the following
syntax:

model.nls1 <-
nls(y ~ response(10, b[subject.id], c[subject.id], alpha1, alpha2,
indicator, t, t2),
data = foo.frame,
start = list(b = rep(25,12), c = rep(100,12),
alpha1 = 0.5, alpha2 = 0.5),
trace = T)

The start values were motivated by some data exploration, and
the results appear to be stable.  The value "a=10" was fixed also
as a result of the initial data exploration, and appears necessary
in order for the model to be stable.

Unfortunately, the estimated b- and c-values for several subjects are
negative.  Also, nls( ) does not allow a "weights = " statement
like gnls( ) does.  When I try

model.nls1 <-
gnls(y ~ response(10, b[subject.id], c[subject.id], alpha1, alpha2,
indicator, t, t2),
data = foo.frame,
start = list(b = rep(25,12), c = rep(100,12),
alpha1 = 0.5, alpha2 = 0.5),
trace = T)

I get the message "Error in eval(expr, envir, enclos) : Object "b" not
found"
This surprises me, since my understanding is that gnls( ) is essentially
nls( ) but with "weights = " and "correlation = " options.  I suppose
that separate fixed effects for each subject could be estimated from
gnls( ) if I created a separate indicator variable for each subject and
added them to the data frame (I have not yet done this); however, this
does not address the need for the b,c parameters to be constrained
greater than zero.

I would gratefully welcome suggestions.