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