[R] Capturing warning within user-defined function
Bert Gunter
bgunter.4567 at gmail.com
Tue Mar 6 23:45:30 CET 2018
1. I did not attempt to sort through your voluminous code. But I suspect
you are trying to reinvent wheels.
2. I don't understand this:
"I've failed to find a solution after much searching of various R related
forums."
A web search on "error handling in R" **immediately** brought up ?tryCatch,
which I think is what you want.
If not, you should probably explain why it isn't, so that someone with more
patience than I can muster will sort through your code to help.
Cheers,
Bert
Bert Gunter
"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
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