2
votes

I would like to display results of two regression analyses next to each other, let's say a logistic and a COX regression. Variables are presented in rows and corresponding data for p, OR/HR and confidence interval in columns. So, the column names don't match: it's OR on the left, HR on the right.

I tried cbind, but I'm encountering the following problems:

(1) What if one variable (row) is only present in one table due to variable selection? In that case, there would have to a blank row in the other table upon combination of the two.

(2) After both tables are combined, how can you add a header spanning over the four columns of each table in order to label them as logistic or COX?

The output that I'm looking for that is to be displayed as a latex table via xtable:

   | Logistic regression          | COX regression
   | p     | OR    | 2.5% | 97.5% | p     | HR    | 2.5% | 97.5% 
v1 | 0.849 | 0.936 | ...  | ...   |       |       |      |
v2 | 0.249 | 0.595 | ...  | ...   | 0.026 | 1.916 | ...  | ...
v3 |       |       |      |       | 0.023 | 0.140 | ...  | ...

Here is some sample code:

library(xtable);library(survival)
# Creating a sample data frame
set.seed(1234)
event <- as.numeric(round(rnorm(10,5,5))>1);tte<-abs(round(rnorm(10,5,5)))
v1 <- round(rnorm(10,5,5));v2 <- round(rnorm(10,2,5));v3 <- round(rnorm(10,1,5))
df<-data.frame(event, tte, v1, v2, v3)

# Some logistic regression ...
LogReg <- glm(event~v1+v2,family=binomial,data=df)
LogRegTable<-round(cbind(summary(LogReg)$coef[, "Pr(>|z|)"], exp(coef(LogReg)),   exp(confint(LogReg))),3)
colnames(LogRegTable)<-c("p","OR","2.5%","97.5%")
LogRegTable<-LogRegTable[!rownames(LogRegTable)=="(Intercept)",] 

# ... and some COX regression
CoxReg <- coxph(Surv(df$tte, df$event)~v2+v3,data=df)
CoxRegTable<-round(cbind(summary(CoxReg)$coef[, "Pr(>|z|)"], exp(coef(CoxReg)),  exp(confint(CoxReg))),3)
colnames(CoxRegTable)<-c("p","HR","2.5%","97.5%")

# There we go
LogRegTable
CoxRegTable

# Now, how to get them in one table?
xtable(
  cbind(LogRegTable, CoxRegTable)
)
# ... messes up the correct row names   
4
If you want to match side-by-side based on variable names, merge() would be better than cbind. If you want fancy tables, consider using a package like xtables or stargazer to make LaTeX output (or write code to make the LaTeX exactly how you want it). There's simply not enough information here to make more specific recommendations (consider adding sample data and code you've tried)MrFlick
Thanks. In fact, I want to pass the final output through xtable to make it look good in LaTex. When I use merge(), how can I specify to match via variable names, i.e. the first column? Listing column names with "colnames(LogRegTable)" starts with "p", so the second column.moabit21

4 Answers

2
votes

As MrFlick suggests, you can do this with merge() more effectively than with cbind.

You can use merge() as follows:

xtable <- merge(LogRegTable, CoxRegTable, by.x = "p", by.y = "p", all = TRUE)

In this case, the by.x and by.y reference the column you want to merge on. If the column name is the same in both tables (as it is in this case) then you can just use by = "p" but this method (with by.x and by.y) allows you to merge two different column headings.

The all = TRUE will force the merged table to keep all rows from both tables, even if they don't both share a value for p. If you changed this to FALSE, then it would remove any rows with a value for p that doesn't exist in both tables.

If you want to rename the column names after (although your column names should have been retained satisfactorily) you can use colnames() as follows:

colnames(xtable) <- c("col1_name", "col2_name", etc)
1
votes

Here is how I would do it:

#Add column for variable name 
LogRegTable <- cbind(v = rownames(LogRegTable), LogRegTable)
CoxRegTable <- cbind(v = rownames(CoxRegTable), CoxRegTable)

#merge on variable name
bothTable <- merge(LogRegTable, CoxRegTable, by = "v", all = TRUE)

#create an xtable object
forLatex <- xtable(bothTable)

#remove rownames and formatting details which I personally prefer to set up in my actual latex code.
#customize based on ?print.xtable to fit your needs
#and write the table to a text file
cat( print.xtable(forLatex, include.colnames = FALSE, include.rownames = FALSE, only.contents = TRUE), file = "latextable.txt")

Hope that helps.

1
votes

R doesn't really do the "double header" thing you're asking for (where columns names can be grouped under a header name), but this will get you to something manageable:

You'll have to convert to data.frame for this to work, which means your column names that start with a number will get an X in front of them. You can change this later for xtable. You should add a column that identifies the variable:

> LogRegTable = data.frame(LogRegTable, Variable = rownames(LogRegTable))
> CoxRegTable = data.frame(CoxRegTable, Variable = rownames(CoxRegTable))
> LogRegTable
       p    OR X2.5. X97.5. Variable
v1 0.849 0.936 0.433  2.018       v1
v2 0.249 0.595 0.129  0.995       v2
> CoxRegTable
       p    HR X2.5. X97.5. Variable
v2 0.026 1.916 1.080  3.402       v2
v3 0.023 0.140 0.025  0.766       v3  

Then it's a straightforward merge operation:

> merge(LogRegTable, CoxRegTable, by="Variable", all=T, suffixes = c("Log", "Cox"))
  Variable  pLog    OR X2.5.Log X97.5.Log  pCox    HR X2.5.Cox X97.5.Cox
1       v1 0.849 0.936    0.433     2.018    NA    NA       NA        NA
2       v2 0.249 0.595    0.129     0.995 0.026 1.916    1.080     3.402
3       v3    NA    NA       NA        NA 0.023 0.140    0.025     0.766
1
votes

If you are eventually outputing latex table. You don't have to separate the R output from latex. I recommend using texreg and knitr. As mentionned above, latex output can both be customized using Stargazer or xtable. Yet from my point of view texreg is more efficient that stargazer and easier to customize than xtable. That said I cannot reproduce exactly the same format you require as it seems that the p-values does not seem to be directly available. But the confidence interval is starred if the estimation is outside the confidence interval.

previous version

texreg don't compute odd ratios directly, but it is easy to customize. Once texreg package loaded you can access the extract functions by typing. (In case you prefer log of odds then you can skip directly to the And we are good to go line)

extract.glm and extract.coxph in the command line

and copying it in your script. Then we only have to modify

coefficients <- s$coef[, 1] into coefficients <- exp(s$coef[, 1]) for both function.

It is good to rename the functions as well just in case something goes wrong. I chose

extract.coxphOR and extract.glmOR

After running the function we have to tell texreg to use our modified function instead of the regular ones.

 setMethod("extract", signature = className("glm","stats"),definition = extract.glmOR)
 setMethod("extract", signature = className("coxph","survival"),definition = extract.coxphOR)

And we are good to go:

texreg(list(LogReg,CoxReg),ci.force = T,single.row=T,ci.force.level=0.025,custom.model.names = c("Logistic Regression", "Cox Regression"))

updated version as suggested in the comment, the odds ratios can be computed in an easier way. Though the previous version is not efficient I leave is in case someone needs to modify the extract function.

texreg(list(LogReg,CoxReg),ci.force = TRUE, single.row = TRUE, ci.force.level = 0.025,    custom.model.names = c("Logistic Regression", "Cox Regression"), override.coef = list(exp(coef(LogReg)), exp(coef(CoxReg))))

results:

\begin{table}
\begin{center}
\begin{tabular}{l c c }
\hline
           & Logistic Regression & Cox Regression \\
\hline
(Intercept)    & $1.77 \; [1.72;\ 1.82]^{*}$ &                             \\
v1             & $0.94 \; [0.92;\ 0.95]^{*}$ &                             \\
v2             & $0.59 \; [0.58;\ 0.61]^{*}$ & $1.92 \; [1.91;\ 1.93]^{*}$ \\
v3             &                             & $0.14 \; [0.11;\ 0.17]^{*}$ \\
\hline
AIC            & 13.63                       & 17.36                       \\
BIC            & 14.54                       &                             \\
Log Likelihood & -3.82                       &                             \\
Deviance       & 7.63                        &                             \\
Num. obs.      & 10                          & 10                          \\
R$^2$          &                             & 0.68                        \\
Max. R$^2$     &                             & 0.92                        \\
Num. events    &                             & 7                           \\
Missings       &                             & 0                           \\
PH test        &                             & 0.90                        \\
\hline
\multicolumn{3}{l}{\scriptsize{$^*$ 0 outside the confidence interval}}
\end{tabular}
\caption{Statistical models}
\label{table:coefficients}
\end{center}
\end{table}