Solutions des TPs
645
C- Cr´ eation d’une classe lm2 pour la r´ egression lin´ eaire avec deux
r´ egresseurs
lm2 <- function(...) {
obj <- lm(...)
if(ncol(model.frame(obj))!=3)
stop("two independent variables are required!")
class(obj) <- c("lm2",class(obj)) # ou c("lm2","lm")
obj
}
n <- 20
x1 <- runif(n,-5,5)
x2 <- runif(n,-50,50)
y <- .3+2 * x1+2 * x2+rnorm(n,0,20)
reg2 <- lm2(y∼x1+x2)
summary(reg2)
plot3d.lm2 <- function(obj,radius=1,lines=TRUE,windowRect,...) {
matreg <- model.frame(obj)
colnames(matreg) <- c("y","x1","x2")
predlim <- cbind(c(range(matreg[,2]),rev(range(matreg[,2]))),
rep(range(matreg[,3]),c(2,2)))
predlim <- cbind(predlim,apply(predlim,1,
function(l) sum(c(1,l) * coef(obj))
))
if(missing(windowRect)) windowRect=c(2,2,500,500)
open3d(windowRect=windowRect,...)
bg3d(color = "gray")
plot3d(formula(obj),type="n")
spheres3d(formula(obj),radius=radius,specular="green")
quads3d(predlim,color="blue",alpha=0.7,shininess=128)
quads3d(predlim,color="cyan",size=5,front="lines",back="lines",
lit=FALSE)
if(lines) {
matpred <- cbind(matreg[2:3],model.matrix(obj)% * %coef(obj))
points3d(matpred)
colnames(matpred) <- c("x1","x2","y")
matlines <- rbind(matreg[,c(2:3,1)],matpred)
nr <- nrow(matreg)
matlines <- matlines[rep(1:nr,rep(2,nr))+c(0,nr),]
segments3d(matlines)
}
}