9
9
# ' 'rotate', 'perm', 'pboot' and 'boot' are defined by \code{\link{resid_rotate}},
10
10
# ' \code{\link{resid_perm}}, \code{\link{resid_pboot}} and \code{\link{resid_boot}}
11
11
# ' respectively
12
- # ' @param ... other arguments passedd onto \code{method}.
12
+ # ' @param additional whether to compute additional measures: standardized
13
+ # ' residuals and leverage
14
+ # ' @param ... other arguments passed onto \code{method}.
13
15
# ' @return a function that given \code{data} generates a null data set.
14
16
# ' For use with \code{\link{lineup}} or \code{\link{rorschach}}
15
17
# ' @export
16
- # ' @importFrom stats lm predict
18
+ # ' @importFrom stats lm predict deviance df.residual lm.influence
17
19
# ' @seealso null_permute, null_dist
18
20
# ' @examples
19
21
# ' data(tips)
23
25
# ' ggplot(lineup(null_lm(tip ~ total_bill, method = 'rotate'), tips.reg)) +
24
26
# ' geom_point(aes(x = total_bill, y = .resid)) +
25
27
# ' facet_wrap(~ .sample)
26
- null_lm <- function (f , method = " rotate" , ... ) {
28
+ null_lm <- function (f , method = " rotate" , additional = FALSE , ... ) {
27
29
n <- NULL
28
30
if (is.character(method )) {
29
31
method <- match.fun(paste(" resid" , method , sep = " _" ))
@@ -33,9 +35,15 @@ null_lm <- function(f, method = "rotate", ...) {
33
35
resp_var <- all.vars(f [[2 ]])
34
36
35
37
resid <- method(model , df , ... )
36
- fitted <- stats :: predict(model , df )
38
+ fitted <- predict(model , df )
37
39
df [" .resid" ] <- resid
38
40
df [" .fitted" ] <- fitted
41
+ if (additional ){
42
+ s <- sqrt(deviance(model )/ df.residual(model ))
43
+ hii <- lm.influence(model , do.coef = FALSE )$ hat
44
+ df [" .leverage" ] <- dropInf(hii , hii )
45
+ df [" .stdresid" ] <- dropInf(resid / (s * sqrt(1 - hii )), hii )
46
+ }
39
47
df [[resp_var ]] <- fitted + resid
40
48
df
41
49
}
@@ -103,7 +111,20 @@ resid_boot <- function(model, data) {
103
111
# '
104
112
# ' @param model to extract residuals from
105
113
# ' @importFrom stats resid
114
+ # ' @param data used to fit model
106
115
# ' @export
107
116
resid_perm <- function (model , data ) {
108
117
sample(stats :: resid(model ))
109
118
}
119
+
120
+
121
+ # Helper function for leverages, adapted from plot.lm
122
+ dropInf <- function (x , h ) {
123
+ if (any(isInf <- h > = 1 )) {
124
+ warning(gettextf(" not plotting observations with leverage greater than one:\n %s" ,
125
+ paste(which(isInf ), collapse = " , " )), call. = FALSE ,
126
+ domain = NA )
127
+ x [isInf ] <- NaN
128
+ }
129
+ x
130
+ }
0 commit comments