On 2/9/2012 6:24 PM, ilai wrote:
You do not provide mlm.influence() so your code can't be reproduced.
Or did you mean to put lm.influence() in the formals to your hatvalues.mlm ?
If yes, then 1) you have a typo 2) lm.influence doesn't allow you to
pass on arguments, maybe try influence.lm instead.
No, I've written an S3 method, influence.mlm and a computational method,
mlm.influence, both of which take an m= argument. I didn't post all the
code because I thought the question might have an easy answer based on
what I provided.
I include all the code and a test case in the attached .txt file that
can be sourced.
-Michael
Elai
On Thu, Feb 9, 2012 at 1:42 PM, Michael Friendly<frien...@yorku.ca> wrote:
I'm trying to write some functions extending influence measures to
multivariate linear models and also
allow subsets of size m>=1 to be considered for deletion diagnostics. I'd
like these to work roughly parallel
to those functions for the univariate lm where only single case deletion
(m=1) diagnostics are considered.
Corresponding to stats::hatvalues.lm, the S3 method for class "lm" objects,
hatvalues<-function (model, ...)
UseMethod("hatvalues")
hatvalues.lm<-
function (model, infl = lm.influence(model, do.coef = FALSE), ...)
{
hat<- infl$hat
names(hat)<- names(infl$wt.res)
hat
}
I have, for class "mlm" objects
hatvalues.mlm<- function(model, m=1, infl=mlm.influence(model, m=m, do.coef
= FALSE), ...)
{
hat<- infl$H
m<- infl$m
names(hat)<- if(m==1) infl$subsets else apply(infl$subsets,1, paste,
collapse=',')
hat
}
where mlm.influence() does the calculations, but also allows the m= argument
to specify subset size.
Yet when I test this I can't seem to pass the m= argument directly, so that
it gets stuffed in to the infl=
call to mlm.influence.
# fit an mlm
library(heplots)
Rohwer2<- subset(Rohwer, subset=group==2)
rownames(Rohwer2)<- 1:nrow(Rohwer2)
Rohwer.mod<- lm(cbind(SAT, PPVT, Raven) ~ n+s+ns+na+ss, data=Rohwer2)
class(Rohwer.mod)
[1] "mlm" "lm"
## this doesn't work, as I would like it to, calling the hatvalues.mlm
method, but passing m=2:
hatvalues(Rohwer.mod, m=2)
Error in UseMethod("hatvalues") :
no applicable method for 'hatvalues' applied to an object of class
"c('double', 'numeric')"
I don't understand why this doesn't just call hatvalues.mlm, since
Rohwer.mod is of class "mlm".
# These work -- calling hatvalues.mlm explicitly, or passing the infl=
argument with the call to
# mlm.influence
hatvalues.mlm(Rohwer.mod, m=2)
hatvalues(Rohwer.mod, infl=mlm.influence(Rohwer.mod,m=2))
Can someone help me understand what is wrong and how to make the .mlm method
allow m= to be passed
directly to the infl= computation?
thx,
-Michael
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street Web: http://www.datavis.ca
Toronto, ONT M3J 1P3 CANADA
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
--
Michael Friendly Email: friendly AT yorku DOT ca
Professor, Psychology Dept.
York University Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street Web: http://www.datavis.ca
Toronto, ONT M3J 1P3 CANADA
# S3 method for class "mlm"
influence.mlm <- function(model, m=1, do.coef=TRUE, ...)
mlm.influence(model, m=m, do.coef = do.coef, ...)
# main computation a la lm.influence
mlm.influence <-
function (model, m=1, do.coef = TRUE, ...)
{
# helper functions
vec <- function(M) {
R <- matrix(M, ncol=1)
if (is.vector(M)) return(R)
nn<-expand.grid(dimnames(M))[,2:1]
rownames(R) <- apply(as.matrix(nn), 1, paste, collapse=":")
R
}
X <- model.matrix(model)
data <- model.frame(model)
Y <- as.matrix(model.response(data))
r <- ncol(Y)
n <- nrow(X)
p <- ncol(X)
labels <- rownames(X)
call <- model$call
B <- coef(model)
E <- residuals(model)
XPXI <- solve(crossprod(X))
EPEI <- solve(crossprod(E))
vB <- vec(t(B))
S <- crossprod(E)/(n-p);
V <- solve(p*S) %x% crossprod(X)
subsets <- t(combn(n, m))
nsub <- nrow(subsets)
Beta <- as.list(rep(0, nsub))
R <- L <- H <- Q <- as.list(rep(0, nsub))
CookD <- as.vector(rep(0, nsub))
# FIXME: naive, direct computations can use qr() or update()
for (i in seq(nsub)) {
I <- c(subsets[i,])
rows <- which(!(1:n) %in% I)
XI <- X[rows,]
YI <- Y[rows,]
BI <- solve(crossprod(XI)) %*% t(XI) %*% YI
EI <- (Y - X %*% BI)[I, , drop=FALSE]
CookD[i] <- t(vec(B - BI)) %*% V %*% vec(B - BI)
H[[i]] <- X[I, , drop=FALSE] %*% XPXI %*% t(X[I, , drop=FALSE])
Q[[i]] <- EI %*% EPEI %*% t(EI)
if (do.coef) Beta[[i]] <- BI
L[[i]] <- if(m==1) H[[i]] / (1-H[[i]]) else H[[i]] %*%
solve(diag(m) - H[[i]])
R[[i]] <- if(m==1) Q[[i]] / (1-H[[i]])
else {
IH <- mpower(diag(m)-H[[i]], -1/2)
IH %*% Q[[i]] %*% IH
}
}
if(m==1) {
H <- unlist(H)
Q <- unlist(Q)
L <- unlist(L)
R <- unlist(R)
subsets <- c(subsets)
}
result <- list(m=m, H=H, Q=Q, CookD=CookD, L=L, R=R, subsets=subsets,
labels=labels, call=call)
if (do.coef) result <- c(result, list(Beta=Beta))
class(result) <-"inflmlm"
result
}
#####################
# extractor functions -- trying to find a way to be consistent with the .lm
versions
hatvalues.mlm <- function(model, m=1, infl=mlm.influence(model, m=m, do.coef =
FALSE), ...)
{
hat <- infl$H
m <- infl$m
names(hat) <- if(m==1) infl$subsets else apply(infl$subsets,1, paste,
collapse=',')
hat
}
cooks.distance.mlm <-
function (model, infl = mlm.influence(model, do.coef = FALSE), ...)
{
cookd <- infl$CookD
m <- infl$m
names(cookd) <- if(m==1) infl$subsets else apply(infl$subsets,1, paste,
collapse=',')
cookd
}
TESTME <- TRUE
if(TESTME) {
library(heplots)
Rohwer2 <- subset(Rohwer, subset=group==2)
rownames(Rohwer2)<- 1:nrow(Rohwer2)
Rohwer.mod <- lm(cbind(SAT, PPVT, Raven) ~ n+s+ns+na+ss, data=Rohwer2)
## this doesn't work, as I would like it to, calling the hatvalues.mlm method,
but passing m=2:
hatvalues(Rohwer.mod, m=2)
## Error in UseMethod("hatvalues") :
## no applicable method for 'hatvalues' applied to an object of class
"c('double', 'numeric')"
#I don't understand why this doesn't just call hatvalues.mlm, since Rohwer.mod
is of class "mlm".
# These work -- calling hatvalues.mlm explicitly, or passing the infl= argument
with the call to
# mlm.influence
hatvalues.mlm(Rohwer.mod, m=2)
hatvalues(Rohwer.mod, infl=mlm.influence(Rohwer.mod,m=2))
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
