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
|
# modguide.model2.R:
#
# linmod code from Stephen Milborrow "Guidelines for S3 Regression Models"
# This is called Model 2 in that document.
#
## A simple linear model (extended from Friedrich Leisch's tutorial).
## Functions like print.linmod in the tutorial don't appear in the code below.
linmod <- function(...) UseMethod("linmod")
linmod.fit <- function(x, y) # internal function, not for the casual user
{ # first column of x is the intercept (all 1s)
y <- as.vector(as.matrix(y)) # necessary when y is a data.frame
qx <- qr(x) # QR-decomposition of x
coef <- solve.qr(qx, y) # compute (x'x)^(-1) x'y
df.residual <- nrow(x) - ncol(x) # degrees of freedom
sigma2 <- sum((y - x %*% coef)^2) / df.residual # variance of residuals
vcov <- sigma2 * chol2inv(qx$qr) # covar mat is sigma^2 * (x'x)^(-1)
colnames(vcov) <- rownames(vcov) <- colnames(x)
fitted.values <- qr.fitted(qx, y)
names(fitted.values) <- rownames(x)
fit <- list(coefficients = coef,
residuals = y - fitted.values,
fitted.values = fitted.values,
vcov = vcov,
sigma = sqrt(sigma2),
df.residual = df.residual)
class(fit) <- "linmod"
fit
}
linmod.default <- function(x, y, ...)
{
fit <- linmod.fit(cbind("(Intercept)"=1, as.matrix(x)), y)
fit$call <- match.call()
fit
}
linmod.formula <- function(formula, data=parent.frame(), ...)
{
mf <- model.frame(formula=formula, data=data)
terms <- attr(mf, "terms")
fit <- linmod.fit(model.matrix(terms, mf), model.response(mf))
fit$call <- match.call()
fit$terms <- terms
fit
}
predict.linmod <- function(object, newdata=NULL, ...)
{
if(is.null(newdata))
y <- fitted(object)
else {
if(is.null(object$terms)) # x,y interface
x <- cbind(1, as.matrix(newdata))
else { # formula interface
terms <- delete.response(object$terms)
x <- model.matrix(terms, model.frame(terms, as.data.frame(newdata)))
}
y <- as.vector(x %*% coef(object))
}
y
}
|
