2015-02-13 3 views
0

Я реализую нулевой завышены отрицательное биномиальное в R. Код здесь:Обработка ошибок в R при выполнении теста ассоциации

> ICHP<-read.table("ichip_data_recodeA.raw",header=TRUE) 
ICHPdt<-data.table(ICHP) 
covfile<-read.table("sorted.covfile.to.glm.out",header=TRUE) 
covfiledt<-data.table(covfile) 
library(pscl) 
fhandle<-file("ichip_zi_nb_model_scoretest.csv","a") 
for (i in seq(7, ncol(ICHPdt), 1)) { 
    notna<-which(!is.na(ICHPdt[[i]])) 
    string<-eval(parse(text = paste("ICHPdt$", colnames(ICHPdt)[i], sep=""))) 
    nullglmmod<-zeroinfl(formula=OverllTot0[notna] ~ EurAdmix[notna] + Sex[notna] + DisDurMonths[notna] + BMI[notna] + Group[notna] + SmokingStatus[notna], data=covfiledt, dist="negbin") 
    nullsum<-coef(summary(nullglmmod)) 
    glmmod<-zeroinfl(formula=OverllTot0[notna] ~ EurAdmix[notna] + Sex[notna] + DisDurMonths[notna] + BMI[notna] + Group[notna] + SmokingStatus[notna] + ICHPdt[[i]][notna], data=covfiledt, dist="negbin") 
    summ <- coef(summary(glmmod)) 
    rownames(summ$zero)[8] <- paste0("ICHPdt$", colnames(ICHPdt)[i]) 
    rownames(summ$count)[8] <- paste0("ICHPdt$", colnames(ICHPdt)[i]) 
    writeLines("zero", con=fhandle) 
    writeLines(colnames(ICHPdt)[i], fhandle) 
    write.table(round(summ$zero, 4), file=fhandle) 
    writeLines("count", con=fhandle) 
    writeLines(colnames(ICHPdt)[i], fhandle) 
    write.table(round(summ$count, 4), file=fhandle) 
} 

Ошибки сценария, когда я = 9246, а также вопросы следующие:

Error in solve.default(as.matrix(fit$hessian)) : 
    system is computationally singular: reciprocal condition number = 1.12288e-19 

В целом, мне нужно пройти через ~ 100 000 маркеров, поэтому я должен ожидать ~ 11 таких ошибок.

Я хотел бы помочь реализовать варианты, например, tryCatch() для обнаружения такой ошибки, пропущения этого маркера и перемещения.

+0

Я полагаю, самый простой способ выражения вопрос «как я могу изменить скрипт таким образом, что любое значение, которое приводит к вычислительно-сингулярной системе, будет пропущено, и скрипт перейдет к следующей строке. Я пытаюсь tryCatch, но без особого успеха ... –

ответ

2

Я рекомендую прочитать this page для быстрой грунтовки и this page для более полного объяснения обработки ошибок, и вы должны в конечном итоге прочитать ?conditions, но в общем, есть два способа обработки ошибок. Во-первых, с примерки уловом, как:

AS.NUMERIC <- function(x){ 

    # for use in the warning handler 
    expectedWarning <- FALSE 

    result = tryCatch({ 
     # a calculation that might raise an error or warning 
     as.numeric(x) 
    }, warning = function(w) { 

     # the typical way to identify the type of 
     # warning is via it's message attribure 
     if(grepl('^NAs introduced by coercion',w$message)){ 
      cat('an expected warning was issued\n') 
      # assign the expected value using the scoping assignment 
      expectedWarning <<- TRUE 
     }else 
      cat('an unexpected warning was issued\n') 
     # reissue the warning 
     warning(w) 
    }, error = function(e) { 
     cat('an error occured\n') 
     # similar things go here but for handling errors 

    }, finally = { 
     # stuff goes here that should happen no matter what, 
     # such as closing connections or resetting global 
     # options such as par(ask), etc. 
    }) 

    # you can handle errors similarly 
    if(expectedWarning) 
     result <- 5 

    return(result) 
} 


AS.NUMERIC('5') 
#> [1] 5 
AS.NUMERIC('five') # raises a warning 
#> an expected warning was issued 
#> [1] 5 
#> Warning message: 
#> In doTryCatch(return(expr), name, parentenv, handler) : 
#> NAs introduced by coercion 

Второй способ заключается в использовании try(), который менее нюансами:

x = try(stop('arbitrary error'),# raise an error 
     silent=TRUE) 

# if there is an error, x will be an object with class 'try-error' 
if(inherits(x,'try-error')) 
    # set the default value for x here 
    x = 5 
Смежные вопросы