Paul Johnson
2011-Dec-19 21:15 UTC
[R] pls help to print out first row of terms(model) output in example program
Greetings. I've written a convenience function for multicollinearity diagnosis. I'd like to report to the user the formula that is used in a regression. I get output like this:> mcDiagnose(m1)[1] "The following auxiliary models are being estimated and returned in a list:" [1] "`x1` ~ ." formula(fmla)() [1] "`x2` ~ ." I'd like to fill in the period with the variable names that are in there. I know the information is in there, I just need to get it. The terms output for a fitted model has the output at the very top, in the first line, above the attributes. I just want to print out that first line, here:> terms(m2)y ~ log(10 + x1) + poly(x2, 2) attr(,"variables") list(y, log(10 + x1), poly(x2, 2)) attr(,"factors") log(10 + x1) poly(x2, 2) y 0 0 log(10 + x1) 1 0 poly(x2, 2) 0 1 attr(,"term.labels") [1] "log(10 + x1)" "poly(x2, 2)" [snip] In my working example code below , I need the help where I have "##fix me fix me## ##Paul Johnson ## 2011-12-19 ## mcDiagnose.R lmAuxiliary <- function(model){ dat <- as.data.frame(model.matrix(model)) ## ivnames <- attr(delete.response(terms(model)), "term.labels") ## previous does not work with transforms like poly hasIntercept <- attr(terms(model), "intercept") if (hasIntercept) dat <- dat[ , -1] # removes intercept. assumes intercept in column 1 ivnames <- colnames(dat) print("The following auxiliary models are being estimated and returned in a list:") results <- list() for (i in ivnames) { fmla <- paste( "`", i, "`", " ~ ." , sep=""); print(fmla) maux <- lm( formula(fmla), data=dat) results[[ i ]] <- maux print(maux$call[2]) ###fix me fix me ## } results } mcDiagnose <- function(model){ auxRegs <- lmAuxiliary(model) auxRsq <- numeric(length=length(auxRegs)) j <- 0 for ( i in auxRegs ){ j <- j + 1 sumry <- summary(i) auxRsq[j] <- sumry$r.squared } print("Drum roll please!") print("And your R_j Squareds are (auxiliary Rsq)") print(names(auxRegs)) print(auxRsq) } x1 <- rnorm(100) x2 <- rnorm(100) y <- rnorm(100) m1 <- lm(y~x1+x2) mcDiagnose(m1) m2 <- lm(y~log(10+x1) + poly(x2,2)) mcDiagnose(m2) -- Paul E. Johnson Professor, Political Science 1541 Lilac Lane, Room 504 University of Kansas
William Dunlap
2011-Dec-19 22:13 UTC
[R] pls help to print out first row of terms(model) output in example program
Does applying formula() to the output of lm() (or terms()) do what you want? > d <- data.frame(y=1:10, x1=log(1:10), x2=sqrt(1:10), x3=1/(1:10)) > formula(lm(formula("y~."), data=d)) # dot expanded y ~ x1 + x2 + x3 > lm(formula("y~."), data=d)$call[[2]] # your original formula("y~.") > lm(formula("y~."), data=d)$call$formula # a better version of original formula("y~.") Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com> -----Original Message----- > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On Behalf Of Paul Johnson > Sent: Monday, December 19, 2011 1:15 PM > To: R-help > Subject: [R] pls help to print out first row of terms(model) output in example program > > Greetings. > > I've written a convenience function for multicollinearity diagnosis. > I'd like to report to the user the formula that is used in a > regression. I get output like this: > > > mcDiagnose(m1) > [1] "The following auxiliary models are being estimated and returned in a list:" > [1] "`x1` ~ ." > formula(fmla)() > [1] "`x2` ~ ." > > I'd like to fill in the period with the variable names that are in there. > > I know the information is in there, I just need to get it. The terms > output for a fitted model has the output at the very top, in the first > line, above the attributes. I just want to print out that first line, > here: > > > terms(m2) > y ~ log(10 + x1) + poly(x2, 2) > attr(,"variables") > list(y, log(10 + x1), poly(x2, 2)) > attr(,"factors") > log(10 + x1) poly(x2, 2) > y 0 0 > log(10 + x1) 1 0 > poly(x2, 2) 0 1 > attr(,"term.labels") > [1] "log(10 + x1)" "poly(x2, 2)" > [snip] > > In my working example code below , I need the help where I have "##fix > me fix me## > > > ##Paul Johnson > ## 2011-12-19 > ## mcDiagnose.R > > lmAuxiliary <- function(model){ > dat <- as.data.frame(model.matrix(model)) > ## ivnames <- attr(delete.response(terms(model)), "term.labels") > ## previous does not work with transforms like poly > hasIntercept <- attr(terms(model), "intercept") > > if (hasIntercept) dat <- dat[ , -1] # removes intercept. assumes > intercept in column 1 > ivnames <- colnames(dat) > print("The following auxiliary models are being estimated and > returned in a list:") > > results <- list() > > for (i in ivnames) { > fmla <- paste( "`", i, "`", " ~ ." , sep=""); > print(fmla) > maux <- lm( formula(fmla), data=dat) > results[[ i ]] <- maux > print(maux$call[2]) > ###fix me fix me ## > } > results > } > > > mcDiagnose <- function(model){ > auxRegs <- lmAuxiliary(model) > auxRsq <- numeric(length=length(auxRegs)) > j <- 0 > for ( i in auxRegs ){ > j <- j + 1 > sumry <- summary(i) > auxRsq[j] <- sumry$r.squared > } > print("Drum roll please!") > print("And your R_j Squareds are (auxiliary Rsq)") > print(names(auxRegs)) > print(auxRsq) > } > > x1 <- rnorm(100) > x2 <- rnorm(100) > y <- rnorm(100) > > m1 <- lm(y~x1+x2) > > mcDiagnose(m1) > > > m2 <- lm(y~log(10+x1) + poly(x2,2)) > > mcDiagnose(m2) > > -- > Paul E. Johnson > Professor, Political Science > 1541 Lilac Lane, Room 504 > University of Kansas > > ______________________________________________ > R-help at r-project.org 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.