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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
