Michael Friendly
2013-Jul-03  21:55 UTC
[R] converting a list of loglin terms to a model formula
I'm developing some functions to create symbolic specifications for 
loglinear models of different types.
I don't really know how to 'compute' with model formulas, so
I've done
this in the notation
for stats::loglin(), which is a list of high-order terms in the model.
What I'd like is a function to turn the results of these into a model 
formula, suitable for
MASS::loglm.  That's the reverse of what loglm does.
For example, the simplest versions of models for 3-way tables for joint,
  conditional, and marginal independence can be computed as follows. 
After each, I indicated
the WANTED model formula I'd like from the result
 > joint(3)
$term1
[1] 1 2
$term2
[1] 3
WANTED:  ~ 1:2 + 3
 > condit(3)
$term1
[1] 1 3
$term2
[1] 2 3
WANTED: ~ 1:2 + 2:3
 > mutual(3)
$term1
[1] 1
$term2
[1] 2
$term3
[1] 3
WANTED: ~ 1 + 2 + 3
In case anyone want to play with the code, here are the current, not too 
elegant definitions
of the functions, and some further test cases,
# models of joint independence
   joint <- function(nf, factors=1:nf, with=nf) {
     if (nf == 1) return (list(term1=factors[1]))
     if (nf == 2) return (list(term1=factors[1], term2=factors[2]))
     others <- setdiff(1:nf, with)
     result <- list(term1=factors[others], term2=factors[with])
     result
   }
# conditional independence
   condit <- function(nf, factors=1:nf, with=nf) {
     if (nf == 1) return (list(term1=factors[1]))
     if (nf == 2) return (list(term1=factors[1], term2=factors[2]))
     main <- setdiff(1:nf, with)
     others <- matrix(factors[with], length(with), length(main))
     result <- rbind(factors[main], others)
     result <- as.list(as.data.frame(result, stringsAsFactors=FALSE))
     names(result) <- paste('term', 1:length(result), sep='')
     result
   }
# mutual independence
   mutual <- function(nf, factors=1:nf) {
     result <- sapply(factors[1:nf], list)
     names(result) <- paste('term', 1:length(result), sep='')
     result
   }
### some comparisons
loglin(HairEyeColor, list(c(1, 2), c(1, 3), c(2, 3)))$lrt
loglm(~1:2 + 1:3 +2:3, HairEyeColor)
# use factor names
joint(3, factors=names(dimnames(HairEyeColor)))
condit(3, factors=names(dimnames(HairEyeColor)))
loglin(HairEyeColor, joint(3))$lrt
loglm(~1:2 + 3, HairEyeColor)
loglin(HairEyeColor, condit(3))$lrt
loglm(~1:3 + 2:3, HairEyeColor)
-- 
Michael Friendly     Email: friendly AT yorku DOT ca
Professor, Psychology Dept. & Chair, Quantitative Methods
York University      Voice: 416 736-2100 x66249 Fax: 416 736-5814
4700 Keele Street    Web:   http://www.datavis.ca
Toronto, ONT  M3J 1P3 CANADA
Henrique Dallazuanna
2013-Jul-03  22:38 UTC
[R] converting a list of loglin terms to a model formula
Try this:
 as.formula(sprintf(" ~ %s", do.call(paste, c(lapply(mutual(3), paste,
collapse = ":"), sep = "+"))))
On Wed, Jul 3, 2013 at 6:55 PM, Michael Friendly <friendly@yorku.ca>
wrote:
> I'm developing some functions to create symbolic specifications for
> loglinear models of different types.
> I don't really know how to 'compute' with model formulas, so
I've done
> this in the notation
> for stats::loglin(), which is a list of high-order terms in the model.
>
> What I'd like is a function to turn the results of these into a model
> formula, suitable for
> MASS::loglm.  That's the reverse of what loglm does.
>
> For example, the simplest versions of models for 3-way tables for joint,
>  conditional, and marginal independence can be computed as follows. After
> each, I indicated
> the WANTED model formula I'd like from the result
>
> > joint(3)
> $term1
> [1] 1 2
>
> $term2
> [1] 3
>
> WANTED:  ~ 1:2 + 3
>
> > condit(3)
> $term1
> [1] 1 3
>
> $term2
> [1] 2 3
>
> WANTED: ~ 1:2 + 2:3
>
> > mutual(3)
> $term1
> [1] 1
>
> $term2
> [1] 2
>
> $term3
> [1] 3
>
> WANTED: ~ 1 + 2 + 3
>
> In case anyone want to play with the code, here are the current, not too
> elegant definitions
> of the functions, and some further test cases,
>
> # models of joint independence
>   joint <- function(nf, factors=1:nf, with=nf) {
>     if (nf == 1) return (list(term1=factors[1]))
>     if (nf == 2) return (list(term1=factors[1], term2=factors[2]))
>     others <- setdiff(1:nf, with)
>     result <- list(term1=factors[others], term2=factors[with])
>     result
>   }
> # conditional independence
>   condit <- function(nf, factors=1:nf, with=nf) {
>     if (nf == 1) return (list(term1=factors[1]))
>     if (nf == 2) return (list(term1=factors[1], term2=factors[2]))
>     main <- setdiff(1:nf, with)
>     others <- matrix(factors[with], length(with), length(main))
>     result <- rbind(factors[main], others)
>     result <- as.list(as.data.frame(result, stringsAsFactors=FALSE))
>     names(result) <- paste('term', 1:length(result),
sep='')
>     result
>   }
> # mutual independence
>   mutual <- function(nf, factors=1:nf) {
>     result <- sapply(factors[1:nf], list)
>     names(result) <- paste('term', 1:length(result),
sep='')
>     result
>   }
>
> ### some comparisons
>
> loglin(HairEyeColor, list(c(1, 2), c(1, 3), c(2, 3)))$lrt
> loglm(~1:2 + 1:3 +2:3, HairEyeColor)
>
> # use factor names
> joint(3, factors=names(dimnames(**HairEyeColor)))
> condit(3, factors=names(dimnames(**HairEyeColor)))
>
> loglin(HairEyeColor, joint(3))$lrt
> loglm(~1:2 + 3, HairEyeColor)
>
> loglin(HairEyeColor, condit(3))$lrt
> loglm(~1:3 + 2:3, HairEyeColor)
>
>
>
> --
> Michael Friendly     Email: friendly AT yorku DOT ca
> Professor, Psychology Dept. & Chair, Quantitative Methods
> York University      Voice: 416 736-2100 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<https://stat.ethz.ch/mailman/listinfo/r-help>
> PLEASE do read the posting guide http://www.R-project.org/**
> posting-guide.html <http://www.R-project.org/posting-guide.html>
> and provide commented, minimal, self-contained, reproducible code.
>
-- 
Henrique Dallazuanna
Curitiba-Paraná-Brasil
25° 25' 40" S 49° 16' 22" O
	[[alternative HTML version deleted]]