Publishable Tables for T-Test, ANOVA and Regression Results

knitr::opts_chunk$set(echo = TRUE, warning=FALSE,comment = NA, message=FALSE,
                      fig.height=4, fig.width=6)

Perform the following

Use the data on price pass through in meat to estimate a regression of price change at retail on price change at the farm gate. Report a regression table from R. Interpret the coefficient on price change at the farm gate. Interpret the t-statistics and the p-value reported by R. Is this effect of price change at the farm gate statistically significant at the 5% level?

Install the following packages

library(stargazer)
library(tidyr)
library(ggplot2)
library(plyr)        

Import the dataset

retail <- read.csv("C:\\Users\\user\\Downloads\\retail.csv")

View the first few observations

head(retail,5)

Attacha the dataset

attach(retail)

Estimate the Regression Equation

my_model <- lm(retail$change_retail_value~retail$change_farm_value, data=retail)

View the Results Using Stargazer

stargazer(my_model, report = "vc*stp",type = "text",out = "./q7results.txt")

===============================================
                        Dependent variable:    
                    ---------------------------
                        change_retail_value    
-----------------------------------------------
change_farm_value            0.314***          
                              (0.044)          
                             t = 7.120         
                             p = 0.000         
                                               
Constant                     0.952***          
                              (0.336)          
                             t = 2.834         
                             p = 0.005         
                                               
-----------------------------------------------
Observations                    623            
R2                             0.075           
Adjusted R2                    0.074           
Residual Std. Error      8.369 (df = 621)      
F Statistic           50.689*** (df = 1; 621)  
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

library(texreg)

screenreg(my_model, digits = 2, title = "Regression Results", caption = "Source: retail dataset")

====================================
                          Model 1   
------------------------------------
(Intercept)                 0.95 ** 
                           (0.34)   
retail$change_farm_value    0.31 ***
                           (0.04)   
------------------------------------
R^2                         0.08    
Adj. R^2                    0.07    
Num. obs.                 623       
====================================
*** p < 0.001; ** p < 0.01; * p < 0.05

Presenting Nice Looking T-stat Tables

Install the following library

library(rempsyc)
library(flextable)
library(report)
library(broom)
library(gapminder)
library(xtable)

Load the mtcars data

data("mtcars")

View the first few observations

head(mtcars,6)

Example One

nice_t_test(data = mtcars,
            response = "mpg",
            group = "vs",
            warning = FALSE)

Exammple Two

nice_t_test(
  data = mtcars,
  response = names(mtcars)[1:6],
  group = "am",
  warning = FALSE)

Example Three

t.test.results <- nice_t_test(
  data = mtcars,
  response = names(mtcars)[1:6],
  group = "am",
  warning = FALSE)
t.test.results

Use Flextable to Display Results

my_table <- nice_table(t.test.results)
my_table

One sided T-test

nice_t_test(data = mtcars,
            response = "mpg",
            group = "am",
            alternative = "greater",
            warning = FALSE) |> 
  nice_table()

One Sample T-test

nice_t_test(data = mtcars,
            response = "mpg",
            mu = 17,
            alternative = "less",
            warning = FALSE) |> 
  nice_table()

Paired T-test

Load the dataset

data("ToothGrowth")

View the first few observations

head(ToothGrowth,5)

nice_t_test(data = ToothGrowth,
            response = "len",
            group = "supp",
            paired = TRUE) |> 
  nice_table()

Integrations

There are other ways to do t-tests and format the results properly, should you wish—for example through the broom and report packages. Examples below.

Run the test

model <- t.test(mpg ~ am, alternative="two.sided", data = mtcars)

Broom Table

stats.table <- tidy(model, conf.int = TRUE)

Present the results

nice_table(stats.table, broom = "t.test")

(stats.table <- as.data.frame(report(model)))

nice_table(stats.table, report = "t.test")

Publication Ready ANOVA Table

Load the following libraries

library(readxl)
library(dplyr)
library(tibble)
library(flextable)

Import the dataset

data_anova <- read.csv("C:\\Users\\user\\Downloads\\data_anova.csv")

View the first few observations

head(data_anova,5)

Convert categorical variables to factor variables

data_anova$Rep = as.factor(data_anova$Rep)
data_anova$Water = as.factor(data_anova$Water)
data_anova$Priming = as.factor(data_anova$Priming)

Attach the dataset

attach(data_anova)

DESIGN: RCBD TWO FACTOR FACTORIAL (TWO-FACTOR ANOVA)

for(i in 1:ncol(data_anova[-c(1:3)])) {
          cols <- names(data_anova)[4:ncol(data_anova)]
          aov.model <- lapply(X = cols, FUN = function(x) 
                    aov(reformulate(termlabels = "Rep + Water*Priming", 
                                    response = x), 
                        data = data_anova))
          
          # print df, MS and Pvalue
          final = anova(aov.model[[i]])[,c(1,3,5)]
          
          # Getting rownames
          rnames = rownames(final)
          
          # Setting column names
          colnames(final) = c("DF", "MS", "P-value")
          colnames(final)[2] = cols[i]
          
          # Rounding values to 3 decimal place
          final = as.data.frame(round(final, digits = 2))
          
          # Assign astericks according to p values
          final$sign[final$`P-value` < 0.1] <- "*" 
          final$sign[final$`P-value` < 0.05] <- "**" 
          final$sign[final$`P-value` < 0.01] <- "***"
          final$sign[final$`P-value` > 0.01] <- "ns"
          
          # Merge MS and significance column together
          final[[2]] = paste(final[[2]], 
                             ifelse(is.na(final[[4]]), "", final[[4]]))
          final = final[-c(3,4)]
          anova = writexl::write_xlsx(final, 
                                      path = paste(cols[i], '-ANOVA.xlsx'))
          
          # Print final ANOVA table ----
          file.list <- list.files(pattern='*-ANOVA.xlsx')
          
          df.list <- lapply(X = file.list, FUN = read_excel)
          
          # Combined ANOVA table for all variables
          aov.table = rlist::list.cbind(df.list)
          
          # Remove duplicate columns for DF
          dup.cols = which(duplicated(names(aov.table)))
          aov.table = aov.table[,-dup.cols]

          # Names for sources of variation in ANOVA
          rownames(aov.table) = rnames
}          

Printing ANOVA table

table = flextable(data = aov.table %>%
                    rownames_to_column("SOV")) 

bold(table, bold = TRUE, part = "header")