/
ols-f-test.R
145 lines (120 loc) · 3.63 KB
/
ols-f-test.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
#' F test
#'
#' @description
#' Test for heteroskedasticity under the assumption that the errors are
#' independent and identically distributed (i.i.d.).
#'
#' @param model An object of class \code{lm}.
#' @param fitted_values Logical; if TRUE, use fitted values of regression model.
#' @param rhs Logical; if TRUE, specifies that tests for heteroskedasticity be
#' performed for the right-hand-side (explanatory) variables of the fitted
#' regression model.
#' @param vars Variables to be used for for heteroskedasticity test.
#' @param ... Other arguments.
#'
#' @return \code{ols_test_f} returns an object of class \code{"ols_test_f"}.
#' An object of class \code{"ols_test_f"} is a list containing the
#' following components:
#'
#' \item{f}{f statistic}
#' \item{p}{p-value of \code{f}}
#' \item{fv}{fitted values of the regression model}
#' \item{rhs}{names of explanatory variables of fitted regression model}
#' \item{numdf}{numerator degrees of freedom}
#' \item{dendf}{denominator degrees of freedom}
#' \item{vars}{variables to be used for heteroskedasticity test}
#' \item{resp}{response variable}
#' \item{preds}{predictors}
#'
#' @references
#' Wooldridge, J. M. 2013. Introductory Econometrics: A Modern Approach. 5th ed. Mason, OH: South-Western.
#'
#' @examples
#' # model
#' model <- lm(mpg ~ disp + hp + wt + qsec, data = mtcars)
#'
#' # using fitted values
#' ols_test_f(model)
#'
#' # using all predictors of the model
#' ols_test_f(model, rhs = TRUE)
#'
#' # using fitted values
#' ols_test_f(model, vars = c('disp', 'hp'))
#'
#' @family heteroskedasticity tests
#'
#' @importFrom stats pf
#'
#' @export
#'
ols_test_f <- function(model, fitted_values = TRUE, rhs = FALSE, vars = NULL, ...) UseMethod("ols_test_f")
#' @export
#'
ols_test_f.default <- function(model, fitted_values = TRUE, rhs = FALSE, vars = NULL, ...) {
check_model(model)
check_logic(fitted_values)
check_logic(rhs)
if (length(vars) > 0) {
check_modelvars(model, vars)
fitted_values <- FALSE
}
l <- ols_prep_avplot_data(model)
nam <- names(l)[-1]
resp <- names(l)[1]
n <- nrow(l)
if (rhs) {
fitted_values <- FALSE
k <- frhs(nam, model, n, l)
result <- ftest_result(k)
} else {
if (fitted_values) {
k <- ffit(model)
result <- ftest_result(k)
} else {
k <- fvar(n, l, model, vars)
result <- ftest_result(k)
}
}
out <- list(dendf = result$dendf, f = result$f, fv = fitted_values,
numdf = result$numdf, p = result$p, preds = nam,
resp = resp, rhs = rhs, vars = vars)
class(out) <- "ols_test_f"
return(out)
}
#' @export
#'
print.ols_test_f <- function(x, ...) {
print_ftest(x)
}
frhs <- function(nam, model, n, l) {
np <- length(nam)
var_resid <- (model_rss(model) / n) - 1
ind <- (residuals(model) ^ 2) / var_resid
l <- cbind(l, ind)
mdata <- l[-1]
summary(lm(ind ~ ., data = mdata))$fstatistic
}
fvar <- function(n, l, model, vars) {
var_resid <- (model_rss(model) / n) - 1
ind <- (residuals(model) ^ 2) / var_resid
mdata <- l[-1]
dl <- mdata[, vars]
dk <- as.data.frame(cbind(ind, dl))
summary(lm(ind ~ ., data = dk))$fstatistic
}
ffit <- function(model) {
pred <- fitted(model)
pred_len <- length(pred)
resid <- model$residuals ^ 2
avg_resid <- sum(resid) / pred_len
scaled_resid <- resid / avg_resid
summary(lm(scaled_resid ~ pred))$fstatistic
}
ftest_result <- function(k) {
f <- k[[1]]
numdf <- k[[2]]
dendf <- k[[3]]
p <- pf(f, numdf, dendf, lower.tail = F)
list(f = f, numdf = numdf, dendf = dendf, p = p)
}