0

R users, I've still been hashing out bits and pieces related to my initial question as seen here and now I'm quite stuck.

http://stackoverflow.com/questions/12270578/skipping-over-an-error-warning-in-an-lme-loop-in-r

Here is the code using mtcars as an example dataset. I want to save the lower and upper confidence intervals for every response variable as listed below (but not for the intercept, though I still need its other coefficients for the lme but I've got that already), all in one go (my real dataset is very large and I'm trying to automate it as much as possible)

library(log10)
library(nlme) 
library(lattice)

responseVariables = c("mpg",                                   
 "disp",             
 "hp")

carModels <- list()
carModelNames <- list()
coint<- list()  
coint2<- list()
coint$fixed <- list()
lower<- list() 
upper<- list() 
carCIlower <- list()
carCIupper <- list()


for (i in responseVariables){
    print("Doing: ")
    print(i)
    mtcars$tmp <- as.numeric(mtcars[,i])
    tmpLme <- lme( log10(tmp) ~ I(log10(wt)), random = ~1 | carb / gear / am, data=mtcars,na.action=na.omit ) 
    carModels <- append(carModels, list(tmpLme))
    carModelNames <- append(carModelNames,i) 


    coint <- try(intervals(tmpLme))   
        if (inherits(coint, "try-error")) {
        tmpLme <- lme(log10(tmp) ~ log10(wt), random = ~1 | carb / gear, data=mtcars, na.action=na.omit);
        coint <- try(intervals(tmpLme));

    } else if (inherits(coint, "try-error")) {
        tmpLme <- lme(log10(tmp) ~ log10(wt), random = ~1 | carb / gear, data=mtcars, na.action=na.omit, method="ML");
        coint <- try(intervals(tmpLme));

    } else if (inherits(coint, "try-error")) {
        tmpLme <- lme(log10(tmp) ~ log10(wt), random = ~1 | carb, data=mtcars, na.action=na.omit);
        coint <- try(intervals(tmpLme));
           #} 

    coint2<- append(coint, list(tmpLme))   
    lower <- dim(coint2$fixed)[1]   
    upper <- dim(coint2$fixed)[1]
    carCIlower <- append(carCIlower, coint2$fixed[2,1],lower) 
    carCIupper <- append(carCIupper, coint2$fixed[2,3],upper)  
    vs_wt <- cbind(carModelNames , carCIlower , carCIlower )
    }
        }

Currently I can get the CI values if I run the commands for each response variable, but not as part of the loop. The loop doesn't proceed past the coint2 statement. Well it does but it doesn't give me answers for coint2 and beyond. Or if I run those lines again it'll only give me values for the last item in the loop (i=hp). Also I see there's an unbalanced curly bracket (I've commented it out to show you) but if I use it gives me the lower CI for "disp" for the lower and upper CI for all resp variables.

Can someone point out what's missing?

AWC
  • 13
  • 4
  • 1
    I think you will benefit if you learn how to use `debug()`, `browser()` and `traceback()` among others. http://stackoverflow.com/questions/4442518/general-suggestions-for-debugging-r This will allow you to stop and step through the code, line by line, examining each object at a time, giving you surgical precision for writing code. – Roman Luštrik Oct 10 '12 at 11:03
  • 1
    I'm skeptical that your if .. elseif... elseif construct is doing what you want. But as Roman said, use debugging tools, or at the very least try each `lme` version by hand to see what the result is. – Carl Witthoft Oct 10 '12 at 11:29
  • Carl- I have done and it generates results fine. There's definitely a problem with the loop, or at least how I'm referring to the data afterwards. – AWC Oct 10 '12 at 16:39
  • Ok so the if elseif construct definitely isn't doing what I want. I did fix it up better and am now getting confidence intervals of most things, but the loop does break somewhere in the middle and starts to reloop. I think that section needs to be made into a function. FYI unfortunately my code doesn't work with mtcars in the way I wanted, so I've returned to using my original data. – AWC Oct 29 '12 at 18:45
  • Right, so I'm giving up :) I'll have to do the ones that don't work by hand. – AWC Oct 30 '12 at 18:05

0 Answers0