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.
Much thanks in advance,
david paul