Katya Mauff
2010-Mar-25 10:11 UTC
[R] 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]]