# converting a list of loglin terms to a model formula

4 messages
Open this post in threaded view
|

## 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.caToronto, ONT  M3J 1P3 CANADA ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: 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 <[hidden email]> 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 > > ______________________________**________________ > [hidden email] mailing list > https://stat.ethz.ch/mailman/**listinfo/r-help > PLEASE do read the posting guide 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]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.