3

I've a series of lme4 models I want to run in R on different outcomes in different subsets (each model is run on an Intention To Treat (ITT) and Per-Protocol (PP) subset, and I've different outcomes) and I use texreg() to print a LaTeX table to compare the results and print nicely within knitr documents.

Since I've lots of subtly different models to run I figured the sensible thing to do would be to warp my models and texreg() calls into a function, to which end I have written...

pleasant.regression <- function(data       = proportion,
                                time.frame = "September 2013",
                                outcome    = "unscheduled",
                                family     = binomial,
                                caption    = "ITT and PP Linear Mixed Effects Model Coefficients and Standard Errors for \\emph{any} Unscheduled Visits in September 2013.",
                                label    = "logistic-unscheduled",
                                ...){
    ## Require packages 
    require(lme4)
    require(ResourceSelection)
    require(texreg)
    ## Set defaults for texreg tables, can be modified with additional arguments

    texreg.digits        <- 2
    texreg.table         <- TRUE
    texreg.caption.above <- TRUE
    texreg.booktabs      <- TRUE
    texreg.dcolumn       <- TRUE
    texreg.use.packages  <- FALSE
    texreg.float.pos     <- "!htpb"
    texreg.ci.force      <- TRUE
    texreg.ci.test       <- 0
    ## Parse the outcome into a formula with the desired mixed
    ## effects model
    .formula <- reformulate(response = outcome, termlabels = c("allocation", "gender", "age", "n.unscheduled.2012", "(1 | pracid)"))
    ## Logistic Regresion
    if(family == "binomial"){
        ## ITT
        itt.mixed <- lmer(.formula,
                          data   = subset(data,
                                          period == time.frame & age.include == TRUE),
                          family = family)
        ## PP
        pp.mixed <- lmer(.formula,
                         data   = subset(data,
                                         period == time.frame & age.include == TRUE & pp == 1),
                         family = family)
    }
    ## Negative Binomial
    else if(family == "negbin"){
        ## ITT
        itt.mixed <- glmer.nb(.formula,
                              data   = subset(data,
                                              period == time.frame & age.include == TRUE))
        ## PP
        pp.mixed <- glmer.nb(.formula,
                             data   = subset(data,
                                             period == time.frame & age.include == TRUE & pp == 1))
    }
    ## Save table comparing ITT to PP using texreg()
    results <- invisible(texreg(list(itt.mixed, pp.mixed), 
                                custom.model.names = c("ITT", "PP"),
                                custom.coef.names  = c("Intercept", "Allocation (Letter)", "Gender (Female)", "Age", "$N_{Unscheduled}$ September 2012"),
                                digits             = texreg.digits,
                                caption            = caption,
                                table              = texreg.table,
                                caption.above      = texreg.caption.above,
                                label              = label,
                                booktabs           = texreg.booktabs,
                                dcolumn            = texreg.dcolumn,
                                use.packages       = texreg.use.packages,
                                float.pos          = texreg.float.pos,
                                ci.force           = texreg.ci.force,
                                ci.test            = texreg.ci.test))
    return(results)
}

When I make a call to this though texreg() prints out the results from within the function, but doesn't return the table as an object for printing...

> my.results <- pleasant.regression(data       = proportion,
                               time.frame = "September 2013",
                               outcome    = "unscheduled",
                               family     = "binomial",
                               caption    = "ITT and PP Linear Mixed Effects Model Coefficients and Standard Errors for \\emph{any} Unscheduled Visits in September 2013.",
                               label    = "logistic-unscheduled") ## Not expecting any output
Computing profile confidence intervals ...
Confidence intervals not available for this model. Using naive p values instead.
Computing profile confidence intervals ...
Confidence intervals not available for this model. Using naive p values instead.

\begin{table}[!htpb]
\caption{ITT and PP Linear Mixed Effects Model Coefficients and Standard Errors for \emph{any} Unscheduled Visits in September 2013. \textbf{NB} P-values are not explicitly calculated in favour of 95\% confidence intervals}
\begin{center}
\begin{tabular}{l D{.}{.}{5.11}@{} D{.}{.}{5.11}@{} }
\toprule
                                    & \multicolumn{1}{c}{ITT} & \multicolumn{1}{c}{PP} \\
\midrule
Intercept                           & -0.73^{*}       & -0.71^{*}       \\
                                    & [-0.95;\ -0.51] & [-0.95;\ -0.47] \\
Allocation (Letter)                 & -0.11           & -0.12           \\
                                    & [-0.33;\ 0.12]  & [-0.36;\ 0.12]  \\
Gender (Female)                     & 0.06            & 0.06            \\
                                    & [-0.03;\ 0.15]  & [-0.03;\ 0.16]  \\
Age                                 & -0.01           & -0.01           \\
                                    & [-0.03;\ 0.00]  & [-0.03;\ 0.00]  \\
$N_{Unscheduled}$ September 2012    & 1.18^{*}        & 1.15^{*}        \\
                                    & [1.12;\ 1.25]   & [1.08;\ 1.22]   \\
\midrule
AIC                                 & 12828.97        & 11597.78        \\
BIC                                 & 12873.10        & 11641.29        \\
Log Likelihood                      & -6408.49        & -5792.89        \\
Deviance                            & 12816.97        & 11585.78        \\
Num. obs.                           & 11547           & 10415           \\
Number of Practices                 & 142             & 136             \\
Variance : Practice                 & 0.37            & 0.40            \\
Variance : Residual                 & 1.00            & 1.00            \\
\bottomrule
\multicolumn{3}{l}{\scriptsize{$^*$ 0 outside the confidence interval}}
\end{tabular}
\label{logistic-unscheduled}
\end{center}
\end{table}

> print(my.results) ## Would expect this to hold the above table and print it again
NULL

I tried wrapping the texreg() call in invisible() based on a StackOverflow thread here, so that the results weren't printed and were returned but that doesn't seem to have worked as I expected.

I expect I'm missing something obvious but can't work it out, any pointers/suggestions would be gratefully received.

Community
  • 1
  • 1
slackline
  • 2,295
  • 4
  • 28
  • 43

1 Answers1

3

For whatever reason, the texreg() function chooses to print the results directly to the screen via cat() if you do not provide a file name for output. This is not really how most functions in R work, but with user submitted packages, a package author can do whatever they like.

So technically by default texreg() returns nothing. You can have it return a string by setting return.string=TRUE but it will still print to the screen as well. The easiest way to prevent the automatic printing is by wrapping the call with capture.output(). That will suppress the screen output and turn the result in to a character vector with an entry for each line of output. Thus you can change the end of your function to

results <- paste(capture.output(texreg(... rest of code ...)), collapse="\n")
return(results)

That should return the character value you expect and not print to the screen.

MrFlick
  • 195,160
  • 17
  • 277
  • 295
  • Thanks MrFlick, thats sorted out the problem I was having. – slackline Jul 21 '14 at 05:35
  • This works in an interactive terminal, but when I include a call to my function in my Knitr `.Rnw` file it doesn't produce any output. Assigning the returned object (e.g. `my.results <- pleasant.regression()`) doesn't help either since if I then try to display them (`print(my.results)`) I get an error message saying `’print’: Error: object ’my.results’ not found`. I'll have a play around and do more reading and see if I can work out why Knitr doesn't like this. – slackline Jul 21 '14 at 12:56
  • perhaps `knitr` is already intercepting the output. I'm not sure. I don't use `knitr` myself. – MrFlick Jul 21 '14 at 13:00
  • Ok, I'll have a look at this (and appreciate your help despite not using `knitr`), although if I simply have a call to my function and don't assign the results the returned character value isn't displayed either. I expect its something `knitr` specific and probably worthy of its own thread as you've provided a working solution to my original question. Thanks again. – slackline Jul 21 '14 at 13:30
  • It would be easiest to help if you posted a complete knitr document that reproduces the problem. Perhaps you could simplify the function so it's not that large. Keep the the code as minimal as possible while still reproducing the exact error is helpful. – MrFlick Jul 21 '14 at 15:45
  • 1
    Sorry, I do realise the utility of minimal reproducible examples but don't have enough work hours at the moment. I found the problem it was knitr related and caused by including a reference within the function to a file that couldn't be read, now sorted. As an aside, I asked Philip Liefield the author of `texreg` about using `cat()` to return objects and he has already updated his package, see https://r-forge.r-project.org/forum/forum.php?thread_id=28874&forum_id=4325&group_id=1420 – slackline Jul 22 '14 at 12:52