-
Notifications
You must be signed in to change notification settings - Fork 0
/
coef_standardized.R
154 lines (143 loc) · 5.1 KB
/
coef_standardized.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
#' **[!!]** Compute Standardized Regression Coefficients
#'
#' Compute the standardized regression coefficients (beta) from an object of class `lm`).
#'
#' @param obj (`lm` object) A result of function `lm()`.
#'
#' @return
#' Object of class `lm_beta` which is a list with 2 named fields:
#' \itemize{
#' \item `b` named numeric vector with regression coefficients (not standardized);
#' \item `beta` named numeric vector with standardized coefficients from `lm()` model.
#' }
#'
#'
#' @details
#' This function is inspired by function \pkg{QuantPsyc}`::lm.beta()` written by Thomas D. Fletcher. \cr
#' `coef_standardized()` provides standardized coefficients even when interaction members are present. This is achieved by computing whole model matrix (with all right-hand side members of formula used in call of `lm()`) and calculating standard deviations of each regressor (including interaction members) based on these columns. \cr
#' `coef_standardized()` does not fail if intercept is not present.\cr
#'
#' The remaining calculations are the same as in \pkg{QuantPsyc}`::lm.beta()`.
#'
#'
#'
#' @author
#' Vilmantas Gegzna (modified function written by Thomas D. Fletcher).
#'
#' @seealso
#' \itemize{
#' \item [stats::lm()] in package \pkg{stats}.
#' \item [QuantPsyc::lm.beta()] in package \pkg{QuantPsyc}.
#' }
#' @keywords models
#' @export
#' @examples
#'
#' data(USJudgeRatings)
#' us <- USJudgeRatings
#' names(us)
#'
#' lm1 <- lm(CONT ~ INTG + DMNR + log(DILG), data = us)
#' coef_standardized(lm1)
#'
#' lm2 <- lm(CONT ~ INTG + DMNR * DILG, data = us)
#' coef_standardized(lm2)
#'
#' summary(coef_standardized(lm2))
#'
#'
#' # Do not include intercept
#' lm3 <- lm(CONT ~ 0 + INTG, data = us)
#' coef_standardized(lm3)
#' summary(coef_standardized(lm3))
coef_standardized <- function(obj) {
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
checkmate::assert_class(obj, "lm")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
COEFS <- summary(obj)$coef
ind_a <- which(rownames(COEFS) == "(Intercept)")
if (length(ind_a) == 0) {
# If intercept is not present
a <- 0
b <- COEFS[, 1]
b_names <- rownames(COEFS)
} else {
# If intercept is present
a <- COEFS["(Intercept)", 1]
b <- COEFS[-1, 1]
b_names <- rownames(summary(obj)$coef)[-1]
b_names <- rownames(COEFS)[-1]
}
# Extracts all members of right-hand side of formulae:
data_ <- model.matrix(as.formula(obj$call$formula), data = obj$model)
# Make correct order of columns:
data_ <- as.data.frame(data_)[, b_names, drop = FALSE]
# Do the standardization:
sx_ <- sapply(data_, sd)
sy_ <- sapply(obj$model[1], sd)
beta <- b * sx_ / sy_
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
structure(list(a = a, b = b, beta = beta),
class = c("lm_beta", "list")
)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
}
#' @name deprecated_functions
#' @description Depreceted functions. They exist for compatibility with
#' previous versions.
#' @title Deprecated functions in `biostat` package
#' @param obj `lm` object.
#' @export
standardized_coef <- function(obj) {
.Deprecated("coef_standardized")
coef_standardized(obj)
}
#' @rdname coef_standardized
#'
#' @param x `lm_beta` object.
#' @param object `lm_beta` object.
#' @param digits (integer) number of decimal places to round the answer to.
#' Default is 3.
#' @param ... further parameters to `print` method.
#'
#' @export
print.lm_beta <- function(x, ..., digits = 3) {
cat("Standardized Regression Coefficients:\n")
print(unclass(round(x$beta, digits = digits)), ...)
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#' @rdname coef_standardized
#' @export
summary.lm_beta <- function(object, ..., digits = 3) {
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
b <- rm_names(object$b)
beta <- rm_names(object$beta)
rez <- data.frame(
regressor = c("(Intercept)", names(object$beta)),
coeff = c(object$a, b),
standardized_coeff = c(NA, beta),
influence_rank = c(NA, rank(-abs(beta)))
)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
structure(rez,
class = c("lm_beta_summary", "data.frame")
)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#' @rdname coef_standardized
#' @export
print.lm_beta_summary <- function(x, ..., digits = 3) {
cat("Summary of Standardized Regression Coefficients:\n")
x$standardized_coeff <- round(x$standardized_coeff, digits = digits)
print(data.frame(x))
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# # # This part is the same as QuantPsyc::lm.beta() ~~~~~~~~~~~~~~~~~~~~~~
# b <- summary(obj)$coef[-1, 1]
#
# sx <- sapply(obj$model[-1], sd)
# sy <- sapply(obj$model[1], sd)
# beta_0 <- b * sx / sy
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~