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.