[R] Capturing warning within user-defined function
William Dunlap
wdunlap at tibco.com
Tue Mar 6 23:41:59 CET 2018
You can capture warnings by using withCallingHandlers. Here is an example,
its help file has more information.
dataList <- list(
A = data.frame(y=c(TRUE,TRUE,TRUE,FALSE,FALSE), x=1:5),
B = data.frame(y=c(TRUE,TRUE,FALSE,TRUE,FALSE), x=1:5),
C = data.frame(y=c(FALSE,FALSE,TRUE,TRUE,TRUE), x=1:5))
withWarnings <- function(expr) {
.warnings <- NULL # warning handler will append to this using '<<-'
value <- withCallingHandlers(expr,
warning=function(e) {
.warnings <<- c(.warnings,
conditionMessage(e))
invokeRestart("muffleWarning")
})
structure(value, warnings=.warnings)
}
z <- lapply(dataList, function(data) withWarnings(coef(glm(data=data, y ~
x, family=binomial))))
z
The last line produces
> z
$A
(Intercept) x
160.80782 -45.97184
attr(,"warnings")
[1] "glm.fit: fitted probabilities numerically 0 or 1 occurred"
$B
(Intercept) x
3.893967 -1.090426
$C
(Intercept) x
-115.02321 45.97184
attr(,"warnings")
[1] "glm.fit: fitted probabilities numerically 0 or 1 occurred"
and lapply(z, attr, "warnings") will give you the warnings themselves.
Bill Dunlap
TIBCO Software
wdunlap tibco.com
On Tue, Mar 6, 2018 at 2:26 PM, Jen <plessthanpointohfive at gmail.com> wrote:
> Hi, I am trying to automate the creation of tables for some simply
> analyses. There are lots and lots of tables, thus the creation of a
> user-defined function to make and output them to excel.
>
> My problem is that some of the analyses have convergence issues, which I
> want captured and included in the output so the folks looking at them know
> how to view those estimates.
>
> I am successfully able to do this in a straightforward set of steps.
> However, once I place those steps inside a function it fails.
>
> Here's the code (sorry this is a long post):
>
> # create data
> wt <- rgamma(6065, 0.7057511981, 0.0005502062)
> grp <- sample(c(replicate(315, "Group1"), replicate(3672, "Group2"),
> replicate(1080, "Group3"), replicate(998, "Group4")))
> dta <- data.frame(grp, wt)
> head(dta)
> str(dta)
>
> # declare design
> my.svy <- svydesign(ids=~1, weights=~wt, data=dta)
>
> # subset
> grp1 <- subset(my.svy, grp == "Group1")
>
> # set options and clear old warnings
> options(warn=0)
> assign("last.warning", NULL, envir = baseenv())
>
> ## proportions and CIs
> p <- ((svyciprop(~grp, grp1, family=quasibinomial))[1])
>
> # save warnings
> wrn1 <- warnings(p)
>
> ci_l <- (confint(svyciprop(~grp, grp1, family=quasibinomial), 'ci')[1])
> ci_u <- (confint(svyciprop(~grp, grp1, family=quasibinomial), 'ci')[2])
>
> ## sample counts
> n <- unwtd.count(~grp, grp1)[1]
>
> ## combine into table
> overall <- data.frame(n, p, ci_l, ci_u)
> colnames(overall) <- c("counts", "Group1", "LL", "UL")
>
> ## add any warnings
> ind <- length(wrn1)
> ind
>
> if (ind == 0) { msg <- "No warnings" }
> if (ind > 0) {msg <- names(warnings()) }
> overall[1,5] <- msg
>
> print(overall)
>
> Here's the output from the above:
>
> > # set options and clear old warnings
> > options(warn=0)
> > assign("last.warning", NULL, envir = baseenv())
> >
> > ## proportions and CIs
> > p <- ((svyciprop(~grp, grp1, family=quasibinomial))[1])
> Warning message:
> glm.fit: algorithm did not converge
> >
> > # save warnings
> > wrn1 <- warnings(p)
> >
> > ci_l <- (confint(svyciprop(~grp, grp1, family=quasibinomial), 'ci')[1])
> Warning message:
> glm.fit: algorithm did not converge
> > ci_u <- (confint(svyciprop(~grp, grp1, family=quasibinomial), 'ci')[2])
> Warning message:
> glm.fit: algorithm did not converge
> >
> > ## sample counts
> > n <- unwtd.count(~grp, grp1)[1]
> >
> > ## combine into table
> > overall <- data.frame(n, p, ci_l, ci_u)
> > colnames(overall) <- c("counts", "Group1", "LL", "UL")
> >
> > ## add any warnings
> > ind <- length(wrn1)
> > ind
> [1] 1
> >
> > if (ind == 0) { msg <- "No warnings" }
> > if (ind > 0) {msg <- names(warnings()) }
> > overall[1,5] <- msg
> >
> > print(overall)
> counts Group1 LL UL
> V5
> counts 315 2.364636e-12 2.002372e-12 2.792441e-12 glm.fit: algorithm did
> not converge
>
> Here's the function:
>
> est <- function(var) {
>
> ## set up formula
> formula <- paste ("~", var)
>
> ## set options and clear old warning
> options(warn=0)
> assign("last.warning", NULL, envir = baseenv())
>
> ## proportions and CIs
> p <- ((svyciprop(as.formula(formula), grp1, family=quasibinomial))[1])
>
> ## save warnings
> wrn1 <- warnings(p)
>
> ci_l <- (confint(svyciprop(as.formula(formula) , grp1,
> family=quasibinomial), 'ci')[1])
> ci_u <- (confint(svyciprop(as.formula(formula) , grp1,
> family=quasibinomial), 'ci')[2])
>
> ## sample counts
> n <- unwtd.count(as.formula(formula), grp1)[1]
>
> ## combine into table
> overall <- data.frame(n, p, ci_l, ci_u)
> colnames(overall) <- c("counts", "Group1", "LL", "UL")
>
>
> ## add any warnings
> ind <- length(warnings(p))
> print(ind)
>
> if (ind == 0) { msg <- "No warnings" }
> if (ind > 0) {msg <- names(warnings()) }
> overall[1,5] <- msg
>
> print(overall)
>
> }
>
> Here's the output from running the function:
>
> > est("grp")
> [1] 0
> counts Group1 LL UL V5
> counts 315 2.364636e-12 2.002372e-12 2.792441e-12 No warnings
> Warning messages:
> 1: glm.fit: algorithm did not converge
> 2: glm.fit: algorithm did not converge
> 3: glm.fit: algorithm did not converge
>
> So, the warnings are showing up in the output at the end of the function
> but they're not being captured like they are when run outside of the
> function. Note the 0 output from print(ind) and V7 has "No warnings".
> I know a lot of things "behave" differently inside functions. Case in
> point, the use of "as.formula(var)" rather than just "~grp" being passed to
> the function.
>
> I've failed to find a solution after much searching of various R related
> forums. I even posted this to stackoverflow but with no response. So, if
> anyone can help, I'd be appreciative.
>
> (sidenote: I used rgamma to create my sampling weights because that's what
> most resembles the distribution of my weights and it's close enough to
> reproduce the convergence issue. If I used rnorm or even rlnorm or rweibull
> I couldn't reproduce it. Just FYI.)
>
> Best,
>
> Jen
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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.
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list