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
|
## if(getRversion() < "2.13") {
## nobs <- function (object, ...) UseMethod("nobs")
## ## also used for mlm fits *and* lmrob :
## nobs.lm <- function(object, ...)
## if(!is.null(w <- object$weights)) sum(w != 0) else NROW(object$residuals)
## ## for glmrob :
## nobs.glm <- function(object, ...) sum(!is.na(object$residuals))
## }
## Here and in NAMESPACE:
if(getRversion() < "3.1.0") {
## cut'n'paste from R's source src/library/stats/R/confint.R
format.perc <- function(probs, digits)
## Not yet exported, maybe useful in other contexts:
## quantile.default() sometimes uses a version of it
paste(format(100 * probs, trim = TRUE, scientific = FALSE, digits = digits),
"%")
confint.lm <- function(object, parm, level = 0.95, ...)
{
cf <- coef(object)
pnames <- names(cf)
if(missing(parm)) parm <- pnames
else if(is.numeric(parm)) parm <- pnames[parm]
a <- (1 - level)/2
a <- c(a, 1 - a)
fac <- qt(a, object$df.residual) # difference from default method
pct <- format.perc(a, 3)
ci <- array(NA, dim = c(length(parm), 2L),
dimnames = list(parm, pct))
ses <- sqrt(diag(vcov(object)))[parm] # gives NA for aliased parms
ci[] <- cf[parm] + ses %o% fac
ci
}
## cut'n'paste from R's source src/library/stats/R/dummy.coef.R
dummy.coef.lm <- function(object, use.na=FALSE, ...)
{
Terms <- terms(object)
tl <- attr(Terms, "term.labels")
int <- attr(Terms, "intercept")
facs <- attr(Terms, "factors")[-1, , drop=FALSE]
Terms <- delete.response(Terms)
vars <- all.vars(Terms)
xl <- object$xlevels
if(!length(xl)) { # no factors in model
return(as.list(coef(object)))
}
nxl <- setNames(rep.int(1, length(vars)), vars)
tmp <- unlist(lapply(xl, length)) ## ?? vapply(xl, length, 1L)
nxl[names(tmp)] <- tmp
lterms <- apply(facs, 2L, function(x) prod(nxl[x > 0]))
nl <- sum(lterms)
args <- setNames(vector("list", length(vars)), vars)
for(i in vars)
args[[i]] <- if(nxl[[i]] == 1) rep.int(1, nl)
else factor(rep.int(xl[[i]][1L], nl), levels = xl[[i]])
dummy <- do.call("data.frame", args)
pos <- 0
rn <- rep.int(tl, lterms)
rnn <- rep.int("", nl)
for(j in tl) {
i <- vars[facs[, j] > 0]
ifac <- i[nxl[i] > 1]
if(length(ifac) == 0L) { # quantitative factor
rnn[pos+1] <- j
} else if(length(ifac) == 1L) { # main effect
dummy[ pos+1L:lterms[j], ifac ] <- xl[[ifac]]
rnn[ pos+1L:lterms[j] ] <- as.character(xl[[ifac]])
} else { # interaction
tmp <- expand.grid(xl[ifac])
dummy[ pos+1L:lterms[j], ifac ] <- tmp
rnn[ pos+1L:lterms[j] ] <-
apply(as.matrix(tmp), 1L, function(x) paste(x, collapse=":"))
}
pos <- pos + lterms[j]
}
## some terms like poly(x,1) will give problems here, so allow
## NaNs and set to NA afterwards.
mf <- model.frame(Terms, dummy, na.action=function(x)x, xlev=xl)
mm <- model.matrix(Terms, mf, object$contrasts, xl)
if(any(is.na(mm))) {
warning("some terms will have NAs due to the limits of the method")
mm[is.na(mm)] <- NA
}
coef <- object$coefficients
if(!use.na) coef[is.na(coef)] <- 0
asgn <- attr(mm,"assign")
res <- setNames(vector("list", length(tl)), tl)
for(j in seq_along(tl)) {
keep <- asgn == j
ij <- rn == tl[j]
res[[j]] <-
setNames(drop(mm[ij, keep, drop=FALSE] %*% coef[keep]), rnn[ij])
}
if(int > 0) {
res <- c(list("(Intercept)" = coef[int]), res)
}
class(res) <- "dummy_coef"
res
}
}# if R <= 3.1.0
## Not exported; used for faster checking, e.g., on CRAN
doExtras <- function() {
interactive() || nzchar(Sys.getenv("R_robustbase_check_extra")) ||
identical("true", unname(Sys.getenv("R_PKG_CHECKING_doExtras")))
}
if(getRversion() < "3.5") {
isTRUE <- function(x) is.logical(x) && length(x) == 1L && !is.na(x) && x
isFALSE <- function(x) is.logical(x) && length(x) == 1L && !is.na(x) && !x
if(getRversion() < "3.3") {
sigma <- function(object, ...) UseMethod("sigma")
## For completeness, and when comparing with nlrob() results:
sigma.nls <- function(object, ...)
## sqrt ( sum( R_i ^ 2) / (n - p) ) :
sqrt( deviance(object) / (nobs(object) - length(coef(object))) )
}
}
## shortcut -- used often in print() etc:
pasteK <- function(...) paste(..., collapse = ", ")
## stopifnot(..) helper :
is.1num <- function(x) is.numeric(x) && length(x) == 1L
##' return 'x' unless it is NULL where you'd use 'orElse'
`%||%` <- function(x, orElse) if(!is.null(x)) x else orElse
##' Augment a vcov-matrix by NA rows & cols when needed; from */R/src/library/stats/R/vcov.R
.vcov.aliased <-
asNamespace("stats")$.vcov.aliased %||%
function(aliased, vc, complete = TRUE) {
## Checking for "NA coef": "same" code as in print.summary.lm() in ./lm.R :
if(complete && nrow(vc) < (P <- length(aliased)) && any(aliased)) {
## add NA rows and columns in vcov
cn <- names(aliased)
VC <- matrix(NA_real_, P, P, dimnames = list(cn,cn))
j <- which(!aliased)
VC[j,j] <- vc
VC
} else # default
vc
}
## e.g. for once-per-session warnings:
.optEnv <- new.env(parent = emptyenv(), hash = FALSE)
|
