How to add stars for significance level with odds ratio for polr?-CodePudding

My question is at the very bottom of this post.

Here is an example of codes that is very similar to the method I use:

url <- "http://peopleanalytics-regression-book.org/data/soccer.csv"
soccer <- read.csv(url)

head(soccer)
##   discipline n_yellow_25 n_red_25 position result country level
## 1       None           4        1        S      D England     1
## 2       None           2        2        D      W England     2
## 3       None           2        1        M      D England     1
## 4       None           2        1        M      L Germany     1
## 5       None           2        0        S      W Germany     1
## 6       None           3        2        M      W England     1

str(soccer)
## 'data.frame':    2291 obs. of  7 variables:
##  $ discipline : chr  "None" "None" "None" "None" ...
##  $ n_yellow_25: int  4 2 2 2 2 3 4 3 4 3 ...
##  $ n_red_25   : int  1 2 1 1 0 2 2 0 3 3 ...
##  $ position   : chr  "S" "D" "M" "M" ...
##  $ result     : chr  "D" "W" "D" "L" ...
##  $ country    : chr  "England" "England" "England" "Germany" ...
##  $ level      : int  1 2 1 1 1 1 2 1 1 1 ...

# convert discipline to ordered factor
soccer$discipline <- ordered(soccer$discipline, 
                             levels = c("None", "Yellow", "Red"))

# check conversion
str(soccer)

# apply as.factor to four columns
cats <- c("position", "country", "result", "level")
soccer[ ,cats] <- lapply(soccer[ ,cats], as.factor)

# check again
str(soccer)

# run proportional odds model
library(MASS)
model <- polr(
  formula = discipline ~ n_yellow_25   n_red_25   position   
    country   level   result, 
  data = soccer
)

# get summary
summary(model)

## Call:
## polr(formula = discipline ~ n_yellow_25   n_red_25   position   
##     country   level   result, data = soccer)
## 
## Coefficients:
##                   Value Std. Error t value
## n_yellow_25     0.32236    0.03308  9.7456
## n_red_25        0.38324    0.04051  9.4616
## positionM       0.19685    0.11649  1.6899
## positionS      -0.68534    0.15011 -4.5655
## countryGermany  0.13297    0.09360  1.4206
## level2          0.09097    0.09355  0.9724
## resultL         0.48303    0.11195  4.3147
## resultW        -0.73947    0.12129 -6.0966
## 
## Intercepts:
##             Value   Std. Error t value
## None|Yellow  2.5085  0.1918    13.0770
## Yellow|Red   3.9257  0.2057    19.0834
## 
## Residual Deviance: 3444.534 
## AIC: 3464.534


# get coefficients (it's in matrix form)
coefficients <- summary(model)$coefficients

# calculate p-values
p_value <- (1 - pnorm(abs(coefficients[ ,"t value"]), 0, 1))*2

# bind back to coefficients
coefficients <- cbind(coefficients, p_value)

# calculate odds ratios
odds_ratio <- exp(coefficients[ ,"Value"])

# combine with coefficient and p_value
(coefficients <- cbind(
  coefficients[ ,c("Value", "p_value")],
  odds_ratio
))

Doing this i get the following output:

##                      Value      p_value odds_ratio
## n_yellow_25     0.32236030 0.000000e 00  1.3803820
## n_red_25        0.38324333 0.000000e 00  1.4670350
## positionM       0.19684666 9.105456e-02  1.2175573
## positionS      -0.68533697 4.982908e-06  0.5039204
## countryGermany  0.13297173 1.554196e-01  1.1422177
## level2          0.09096627 3.308462e-01  1.0952321
## resultL         0.48303227 1.598459e-05  1.6209822
## resultW        -0.73947295 1.083595e-09  0.4773654
## None|Yellow     2.50850778 0.000000e 00 12.2865822
## Yellow|Red      3.92572124 0.000000e 00 50.6896241

My question
However I want stars with odds ratios. How is this possible with THIS method? If possible I would like to also add the standard error.
I do not want to use modelsummary() or gtsummary()

CodePudding user response：

How about this:

url <- "http://peopleanalytics-regression-book.org/data/soccer.csv"
soccer <- read.csv(url)

soccer$discipline <- ordered(soccer$discipline, 
                             levels = c("None", "Yellow", "Red"))
cats <- c("position", "country", "result", "level")
soccer[ ,cats] <- lapply(soccer[ ,cats], as.factor)

library(MASS)
#> 
#> Attaching package: 'MASS'
#> The following object is masked _by_ '.GlobalEnv':
#> 
#>     cats
model <- polr(
  formula = discipline ~ n_yellow_25   n_red_25   position   
    country   level   result, 
  data = soccer
)

coefficients <- summary(model)$coefficients
#> 
#> Re-fitting to get Hessian

# calculate p-values
p_value <- (1 - pnorm(abs(coefficients[ ,"t value"]), 0, 1))*2

# bind back to coefficients
coefficients <- cbind(coefficients, p_value)

# calculate odds ratios
coefficients <- cbind(coefficients, odds_ratio = exp(coefficients[ ,"Value"]))

# combine with coefficient and p_value
printCoefmat(coefficients[ ,c("Value", "Std. Error", "odds_ratio", "p_value")], 
             P.values=TRUE, 
             has.Pvalue=TRUE)
#>                    Value Std. Error odds_ratio   p_value    
#> n_yellow_25     0.322360   0.033078     1.3804 < 2.2e-16 ***
#> n_red_25        0.383243   0.040505     1.4670 < 2.2e-16 ***
#> positionM       0.196847   0.116487     1.2176   0.09105 .  
#> positionS      -0.685337   0.150112     0.5039 4.983e-06 ***
#> countryGermany  0.132972   0.093599     1.1422   0.15542    
#> level2          0.090966   0.093547     1.0952   0.33085    
#> resultL         0.483032   0.111951     1.6210 1.598e-05 ***
#> resultW        -0.739473   0.121293     0.4774 1.084e-09 ***
#> None|Yellow     2.508508   0.191826    12.2866 < 2.2e-16 ***
#> Yellow|Red      3.925721   0.205714    50.6896 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

^{Created on 2022-12-06 by the reprex package (v2.0.1)}