[R-sig-ME] nlme question...

[Andrew Robinson]

Hi Jeffrey,  

in the first instance, I think that the problem is that "if" is not
vectorised.  Try ifelse() instead.

Cheers

Andrew

On Sat, Aug 28, 2010 at 07:35:42PM -0400, Jeffrey Harring wrote:
> I would like to fit a nonlinear model with "nlme."
> 
> The model is a nonlinear growth model with five time points: y = X*b + 
> e, where design matrix X is defined as
> 
> X= |  1   0  |
>   |  1   1  |
>   |  1   k  |
>   |  1  2k  |
>   |  1  3k  |
> 
> and parameter vector b = (b0, b1). And where "k" is a parameter to be 
> estimated. Of course I also want to estimate the intercept (b0), slope 
> (b1). Error variances (e) with 5 free parameters and random effects 
> covariance matrix  (2x2: for b0 and b1).
> 
> I am not certain how to get the design matrix in the code so that the 
> "nlme" procedure will be able to
> estimate parameter "k" as well as the intercept and slope.
> 
> I have tried the following:
> 
> meanfunc <- function(x,b0,b1,k){
> meangrad <- array(0,c(length(x),2),list(NULL,c("b0","b1")))
> 
> for(i in 1:length(x)){
>    if(x[i]==0){
>     err[i] <- b0 + b1*x[i]
>       meangrad[i,"b0"] <- 1
>       meangrad[i,"b1"] <- x[i]}
>    if(x[i]==1){
>     err[i] <- b0 + b1*x[i]
>       meangrad[i,"b0"] <- 1
>       meangrad[i,"b1"] <- x[i]}
>    if(x[i]==2){
>     err[i] <- b0 + b1*(x[i]-1)*k
>     meangrad[i,"b0"] <- 1
>       meangrad[i,"b1"] <- (x[i]-1)*k}
>    if(x[i]==3){
>       err[i] <- b0 + b1*2*(x[i]-1)*k
>     meangrad[i,"b0"] <- 1
>       meangrad[i,"b1"] <- (x[i]-1)*k}
>    if(x[i]==4){
>       err[i] <- b0 + b1*3*(x[i]-1)*k
>     meangrad[i,"b0"] <- 1
>       meangrad[i,"b1"] <- (x[i]-1)*k}
> }
> #  compute analytical derivatives
>   attr(err,"gradient") <- meangrad
>   err
> }
> 
> data.mlfit <- nlme(y ~ meanfunc(x,b0,b1,k),
>   fixed=list(b0 ~ 1, b1 ~ 1, k ~ 1),
>   random=list(b0 ~ 1, b1 ~ 1),
>   groups = ~subj,
>   data=nd,
>   start=list(fixed=c(300,50,1)),
>   method="ML",verbose=T)
> 
> I get the following error message:
> 
> Error in meangrad[i, "b1"] <- (x[i] - 1) * k :
>  number of items to replace is not a multiple of replacement length
> In addition: Warning messages:
> 1: In err[i] <- b0 + b1 * x[i] :
>  number of items to replace is not a multiple of replacement length
> 2: In err[i] <- b0 + b1 * x[i] :
>  number of items to replace is not a multiple of replacement length
> 3: In err[i] <- b0 + b1 * (x[i] - 1) * k :
>  number of items to replace is not a multiple of replacement length
> 
> Any suggestions or insights would be greatly appreciated.
> 
> Thank you,
> Jeff
> 
> -- 
> **********************************************************
> Jeffrey R. Harring, Assistant Professor
> Department of Measurement, Statistics & Evaluation (EDMS)
> 1230 Benjamin Building
> University of Maryland
> College Park, MD 20742-1115
> 
> Phone: 	301.405.3630
> Fax: 	301.314.9245
> Email: 	harring at umd.edu
> Web:  	http://www.education.umd.edu/EDMS/fac/Harring/webpage.html
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

-- 
Andrew Robinson  
Program Manager, ACERA 
Department of Mathematics and Statistics            Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia               (prefer email)
http://www.ms.unimelb.edu.au/~andrewpr              Fax: +61-3-8344-4599
http://www.acera.unimelb.edu.au/