[R] constructing an lm() formula in a function

Michael Friendly friendly at yorku.ca
Wed Sep 12 15:03:01 CEST 2007


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



More information about the R-help mailing list