CODEWITHSUNDEEP
rsummarysignificance
rr19
rr19

Reputation: 111

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

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()

Upvotes: 0

Views: 430

Answers (1)

DaveArmstrong
DaveArmstrong

Reputation: 22034

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)

Upvotes: 2

Related Questions