R help - Mar 2010 - help with breaking loops used to fit covariates in nlme model building procedure

Dear All
 
I'm attempting to speed up my model building procedure, but need some help
with the loops I've created...please bear with me through the explanation!
 
My basic model call is something like: 
m0sulf.nlme<-nlme(conc~beta0*exp(-beta1*day)+beta2*exp(-beta3*day),
                        data=m0sulf,
                        fixed=(beta0+beta1+beta2+beta3~1),
                        random=pdDiag(beta0+beta1+beta2+beta3~1),
                        start=c(beta0=-300,beta1=15, beta2=300,beta3=0.5),
                        )
 
In order to update this model with covariates I would use something like: 
 
modelname=update(m0sulf.nlme, fixed=list(beta0~X,beta1+beta2+beta3~1),
start=some.vector.of.starting.estimates)
 
where X is the covariate of choice. 
 
I have over 16 covariates, and running each possible combination is tedious and
horribly time-consuming, as I'd be running:
 
model with beta0~X, beta1~X.... (beta0+beta1)~X, (beta0+beta2)~X...
(beta0+beta1+beta2+beta3)~X, for each X!
 
My starting estimates must be of the form (e.g. if beta0~X, beta1..beta3~1): 
 
start=c(m0sulf.fix[1],0,m0sulf.fix[2],m0sulf.fix[3],m0sulf.fix[4]) 
 
where m0sulf.fix[j] are the fixed estimates for the jth parameter in my previous
model.
 
The number of zero's is dependent on the number of levels of the covariate
(1 if covariate is continuous or binary, m-1 for cov with m levels), and their
placement in the vector is dependent on which parameter I'm looking at, so
e.g. if I were looking at the model (beta0~1, beta1~X, beta2+beta3~1) instead, I
would need:
 
start=c(m0sulf.fix[1],m0sulf.fix[2],0,m0sulf.fix[3],m0sulf.nlme) .
 
I've thus created the loop: (for the first tier of modelling-covariate only
on one parameter at a time)
 
vec=vector("list", 1) #ist tier, only need 1 entry
vec <- lapply(1:length(vec), function(vec) vector("list",4)) # 4
parameters
vec <-  lapply(1:length(vec[[1]]), function(vec) vector("list",6))
# maximum 6 levels to variable
vec

        for (p in 1) { # p is 1 only because for now I'm looking at the
first tier only
        for (b in c(1,2,3,4)) { #parameters
 
            for (j in
1:length(m0sulf[,c("preg","trimester","site")]) ) 
{
 
#2 levels, 3 levels and 6 levels respectively, so have all types of covariates
covered-for my dataset

               
k=(m0sulf[,c("preg","trimester","site")][[j]])
 
                y=nlevels(k)
m=abs(y)
 
                
        my.vect=function(b,p,m){
            if (b==1) {
           
c(m0sulf.fix[1],vector(mode="numeric",p*(m-1)),m0sulf.fix[2],m0sulf.fix[3],m0sulf.fix[4])
                } else {
                  if (b==2) {
                 
c(m0sulf.fix[1],m0sulf.fix[2],vector(mode="numeric",p*(m-1)),m0sulf.fix[3],m0sulf.fix[4])
 
                      } else{
                        if (b==3) {
                       
c(m0sulf.fix[1],m0sulf.fix[2],m0sulf.fix[3],vector(mode="numeric",p*(m-1)),m0sulf.fix[4])
 
                            } else {
                              if (b==4) {
                             
c(m0sulf.fix[1],m0sulf.fix[2],m0sulf.fix[3],m0sulf.fix[4],vector(mode="numeric",p*(m-1)))
                              }}}}}
        vec[[p]][[b]][[m]]=try(my.vect(b,p,m),TRUE)
                            }}}
 
 
#which creates: see attach1.txt
 
followed by: again only for the first tier-cov only on one parameter at a time
 
mod<- vector("list", 1)
mod<- lapply(1:length(vec), function(vec) vector("list",4))
mod<-  lapply(1:length(vec[[1]]), function(vec) vector("list",4))
mod
#modelling loop
#1st tier
        for (p in 1) {
            for (b in c(1,2,3,4)) {
 
                for (j in m0sulf[,c("preg", "trimester")]) 
{ #just trying it for 2 of 16 covariates for now
                    X=j
                    y=nlevels(X)
m=abs(y)
                    my.funct=function(conc,day,beta0,beta1,beta2,beta3,X) {
                        if (b==1) {
                        update(m0sulf.nlme,
fixed=list(beta0~X,beta1+beta2+beta3~1), start=vec[[p]][[b]][[m]])
                        }else {
                            if (b==2) {
                            update(m0sulf.nlme,
fixed=list(beta0~1,beta1~X,beta2+beta3~1), start=vec[[p]][[b]][[m]])
                            }else {
                                if (b==3) {
                                update(m0sulf.nlme,
fixed=list(beta0~1,beta1~1,beta2~X,beta3~1), start=vec[[p]][[b]][[m]])
                                }else {
                                    if (b==4) {
                                    update(m0sulf.nlme,
fixed=list(beta0~1,beta1~1,beta2~1,beta3~X), start=vec[[p]][[b]][[m]])
                                    }}
                            }}}
 
        mod=try(my.funct(conc,day,beta0,beta1,beta2,beta3,X),TRUE)
        }}}
 
This runs perfectly-the problem is that some of these models don't work-they
return an error or do not converge. The ones that return an error are dealt
with, because I end up with a list of either model fits or error messages, but
the ones which don't converge take hours to run, and I don't have that
kind of time.
 
what I'd like to do, is break out the current model if it is not converging
within a certain time and move on to the next... any suggestions??
 
Suggestions for simplifying the horrific code above would also be very welcome!
It was the only way I could think to do it, but I know that it's clumsy...
 
Thanks in advance
Katya Mauff
 


 
______________________________________________________________________________________________

UNIVERSITY OF CAPE TOWN 

This e-mail is subject to the UCT ICT policies and e-mail disclaimer published
on our website at http://www.uct.ac.za/about/policies/emaildisclaimer/ or
obtainable from +27 21 650 4500. This e-mail is intended only for the person(s)
to whom it is addressed. If the e-mail has reached you in error, please notify
the author. If you are not the intended recipient of the e-mail you may not use,
disclose, copy, redirect or print the content. If this e-mail is not related to
the business of UCT it is sent by the sender in the sender's individual
capacity.

_____________________________________________________________________________________________________


	[[alternative HTML version deleted]]