On Jul 14, 2010, at 2:16 PM, darckeen wrote:
>
> I was hoping that there might be some way to attach significance code like
> the ones from summary.lm to a dataframe. Anyone know how to do something
> like that. Here is the function i'd like to add that functionality to:
>
>
> add1.coef <- function(model,scope,test="F",p.value=1,order.by.p=FALSE)
> {
> num <- length(model$coefficients)
> add <- add1(model,scope,test=test)
> sub <- subset(add,add$'Pr(F)'<=p.value)
> est <- sapply(rownames(sub), function(x) update(model,paste("~ .
> +",x))$coefficients[num+1])
> ret <- data.frame(est,sub$'Pr(F)')
>
> colnames(ret) <- c("Estimate","Pr(F)")
> rownames(ret) <- rownames(sub)
> ret <- format(ret,digits=1,nsmall=1,scientific=F)
>
> if (order.by.p) { ret <- ret[order(ret$'Pr(F)'),]}
> return(ret)
> }
>
> fscope <- as.formula(paste("dep.sign.up ~ ", paste(names(lr)[2:10],
> collapse= "+")))
> rslt <- add1.coef(lm(dep.sign.up ~ 1,
> data=lr),fscope,p.value=1,order.by.p=FALSE)
>
> Estimate Pr(F)
> r1.pop.total 0.02 0.09
> r1.pop.household 0.05 0.09
> r1.pop.female 0.03 0.09
> r1.pop.pct.female 14594.39 0.35
> r1.pop.male 0.04 0.08
> r1.pop.pct.male -13827.51 0.37
> r1.pop.density 0.06 0.09
> r1.age.0.4.pct 14581.65 0.39
> r1.age.5.14.pct 2849.15 0.76
Review the code for print.summary.lm(), in which you will find the use of
printCoefmat(), in which you will find the use of symnum().
Using lm.D9 from example(lm):
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2,10,20, labels=c("Ctl","Trt"))
weight <- c(ctl, trt)
lm.D9 <- lm(weight ~ group)
> printCoefmat(coef(summary(lm.D9)))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.03200 0.22022 22.8501 9.547e-15 ***
groupTrt -0.37100 0.31143 -1.1913 0.249
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
pv <- coef(summary(lm.D9))[, "Pr(>|t|)"]
> pv
(Intercept) groupTrt
9.547128e-15 2.490232e-01
Signif <- symnum(pv, corr = FALSE, na = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", " "))
> Signif
(Intercept) groupTrt
***
attr(,"legend")
[1] 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Note that 'Signif' is:
> str(Signif)
Class 'noquote' atomic [1:2] ***
..- attr(*, "legend")= chr "0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1"
So you will need to coerce it to a vector before cbind()ing to a data frame:
> as.vector(Signif)
[1] "***" " "
HTH,
Marc Schwartz
______________________________________________
R-help@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.
