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]]