G'day Brian,
>>>>> "BDR" == Prof Brian Ripley <[EMAIL PROTECTED]> writes:
BDR> The proposed fix regretably will not work, since one can do
BDR> things like
BDR> library(MASS)
BDR> prcomp(~ dist + dist:climb, hills)
Yes, I had the impression that this would work with the current
version but dismissed it as unintentional. ;-)
I wouldn't expect anyone to specify interaction terms in a principal
components analysis, but perhaps there are valid reasons to do so.
Attached is a slightly modified patch which still allows the above
example but also the changes that I proposed. On my machine, R-devel
with this patch passes "make check FORCE=FORCE" and produces the
following output:
> library(DAAG)
> res <- prcomp(~ . - case - site - Pop - sex, possum)
> res <- princomp(~ . - case - site - Pop - sex, possum)
> res <- prcomp(~ . - case - site - Pop, possum)
Error in prcomp.formula(~. - case - site - Pop, possum) :
PCA applies only to numerical variables
> res <- princomp(~ . - case - site - Pop, possum)
Error in princomp.formula(~. - case - site - Pop, possum) :
PCA applies only to numerical variables
> library(MASS)
Attaching package: 'MASS'
The following object(s) are masked from package:DAAG :
hills
> prcomp(~ dist + dist:climb, hills)
Standard deviations:
[1] 29655.128279 3.101566
Rotation:
PC1 PC2
dist -0.000154139 -0.999999988
dist:climb -0.999999988 0.000154139
Cheers,
Berwin
Index: src/library/stats/R/princomp.R
===================================================================
--- src/library/stats/R/princomp.R (revision 37574)
+++ src/library/stats/R/princomp.R (working copy)
@@ -10,13 +10,15 @@
mf$... <- NULL
mf[[1]] <- as.name("model.frame")
mf <- eval.parent(mf)
+ na.act <- attr(mf, "na.action")
+ mt <- attr(mf, "terms") # allow model.frame to update it
+ attr(mt, "intercept") <- 0
+ mterms <- attr(mt, "factors")
+ mterms <- rownames(mterms)[apply(mterms, 1, any)]
## this is not a `standard' model-fitting function,
## so no need to consider contrasts or levels
- if(any(sapply(mf, function(x) is.factor(x) || !is.numeric(x))))
+ if(any(sapply(mterms, function(x) is.factor(mf[,x]) ||
!is.numeric(mf[,x]))))
stop("PCA applies only to numerical variables")
- na.act <- attr(mf, "na.action")
- mt <- attr(mf, "terms") # allow model.frame to update it
- attr(mt, "intercept") <- 0
x <- model.matrix(mt, mf)
res <- princomp.default(x, ...)
## fix up call to refer to the generic, but leave arg name as `formula'
Index: src/library/stats/R/prcomp.R
===================================================================
--- src/library/stats/R/prcomp.R (revision 37574)
+++ src/library/stats/R/prcomp.R (working copy)
@@ -37,13 +37,15 @@
mf$... <- NULL
mf[[1]] <- as.name("model.frame")
mf <- eval.parent(mf)
+ na.act <- attr(mf, "na.action")
+ mt <- attr(mf, "terms")
+ attr(mt, "intercept") <- 0
+ mterms <- attr(mt, "factors")
+ mterms <- rownames(mterms)[apply(mterms, 1, any)]
## this is not a `standard' model-fitting function,
## so no need to consider contrasts or levels
- if (any(sapply(mf, function(x) is.factor(x) || !is.numeric(x))))
+ if (any(sapply(mterms, function(x) is.factor(mf[,x]) ||
!is.numeric(mf[,x]))))
stop("PCA applies only to numerical variables")
- na.act <- attr(mf, "na.action")
- mt <- attr(mf, "terms")
- attr(mt, "intercept") <- 0
x <- model.matrix(mt, mf)
res <- prcomp.default(x, ...)
## fix up call to refer to the generic, but leave arg name as `formula'
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
