I'm working on some functions for generalized canonical discriminant
analysis in conjunction with the heplots package. I've written a
candisc.mlm function that takes an mlm object and computes a
candisc object containing canonical scores, coeficients, etc.
But I'm stumped on how to construct a mlm for the canonical scores,
in a function using the *same* right-hand-side of the model formula that
was used in the original mlm for the data.
I can illustrate step-by-step and where I'm stumped.
> # fit randomized block model, .~Block+Gp, where Gp = Contour:Depth
> # Soils data is in car package
> library(car)
> soils.mod1 <- lm(cbind(pH,N,Dens,P,Ca,Mg,K,Na,Conduc) ~ Block + Gp,
+ data=Soils)
> # Do the canonical discriminant analysis for Gp:
> can1 <-candisc(soils.mod1, term="Gp")
>
> str(can1$scores)
'data.frame': 48 obs. of 11 variables:
$ Block: Factor w/ 4 levels
"1","2","3","4": 1 2 3 4 1 2 3 4 1 2 ...
$ Gp : Factor w/ 12 levels "D0","D1","D3",..:
9 9 9 9 10 10 10 10 11
11 ...
$ Can1 : num 8.87 5.57 6.31 7.01 6.64 ...
$ Can2 : num -3.76 -4.78 -4.63 -4.06 -3.13 ...
$ Can3 : num 1.241 -0.561 -1.299 0.642 2.217 ...
$ Can4 : num 1.313 -0.402 -1.631 2.481 0.384 ...
$ Can5 : num -1.913 -1.103 -0.428 1.134 -0.937 ...
$ Can6 : num -0.4219 -0.3593 1.2070 0.0652 -0.6424 ...
$ Can7 : num 0.701 0.243 -0.728 1.147 -0.149 ...
$ Can8 : num 0.0728 1.4491 -0.1546 -1.5191 -0.2374 ...
$ Can9 : num -2.4318 1.2773 -0.0905 1.0646 -1.8069 ...
> str(can1$terms)
chr [1:2] "Block" "Gp"
OK, so now in a function I want to do the equivalent of fitting
lm( cbind(Can1, Can2, ... Can9) ~ Block+Gp, data=can$scores)
The following function shows two ways I've tried to do this, neither of
which works:
fitlm.candisc <- function( can ) {
factors <- can$factors # factor variable from candisc
terms <- can$terms # terms from original lm
scores <- can$scores # scores data fram
# can.mod <-
lm(as.matrix(scores[,paste("Can",1:can$rank,sep="")])
# ~ paste(can1$terms,collapse="+"), data=scores)
can.mod <-
lm(as.matrix(scores[,paste("Can",1:can$rank,sep="")])
~ scores[,terms], data=scores)
}
> fitlm.candisc(can1)
Error in model.frame(formula, rownames, variables, varnames, extras,
extranames, :
invalid type (list) for variable 'scores[, terms]'
>
I would like the terms in the model can.mod to be Block and Gp.
What am I missing here?
thanks,
--
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 http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT M3J 1P3 CANADA
