On Jun 19, 2010, at 7:45 AM, Christos Argyropoulos wrote:
Hi,
mod.poly3$coef/sqrt(diag(mod.poly3$var))
will give you the Wald stat values, so
pnorm(abs(mod.poly3$coef/sqrt(diag(mod.poly3$var))),lower.tail=F)*2
will yield the corresponding p-vals
It will, although it may appear as magic to those untutored in
examining R model objects. Josh B should also consider looking at the
output of str(mod.poly3) and then tracing through the logic used in
the rms/Design function, print.lrm(). It's not a hidden function, so
simply tyoing its name at the console will let him see what steps
Harrell uses. They are a bit different, but are mathematically
equivalent. Stripped of quite a bit of code that is not essential in
this case:
print.lrm.simple <- function (x, digits = 4)
{
# first a couple of utility formatting functions:
sg <- function(x, d) {
oldopt <- options(digits = d)
on.exit(options(oldopt))
format(x)
rn <- function(x, d) format(round(as.single(x), d))
# Then the extraction of compoents from the model, "x"
cof <- x$coef # using name completion, since the full name is x
$coefficients
vv <- diag(x$var) # the diagonal of the variance-covariance matrix
z <- cof/sqrt(vv) # the Wald Z value
stats <- cbind(sg(cof, digits),
sg(sqrt(vv), digits),
rn(cof/sqrt(vv),
2))
stats <- cbind(stats,
# This is the step that calculates the p-values
rn(1 - pchisq(z^2, 1), 4))
#
dimnames(stats) <- list(names(cof), c("Coef", "S.E.", "Wald Z",
"P"))
print(stats, quote = FALSE)
cat("\n")
# the regular print.lrm does not return anything, ... it just prints,
# but if you add this line you will be able to access the components
of:
invisible(stats)
}
> print.lrm.simple(mod.poly3)[ , 4] # still prints first
Coef S.E. Wald Z P
Intercept -5.68583 5.23295 -1.09 0.2772
x1 1.87020 2.14635 0.87 0.3836
x1^2 -0.42494 0.48286 -0.88 0.3788
x1^3 0.02845 0.03120 0.91 0.3618
x2 3.49560 3.54796 0.99 0.3245
x2^2 -0.94888 0.82067 -1.16 0.2476
x2^3 0.06362 0.05098 1.25 0.2121
# the 4th column are the p-values:
Intercept x1 x1^2 x1^3 x2 x2^2 x2^3
"0.2772" "0.3836" "0.3788" "0.3618" "0.3245" "0.2476" "0.2121"
--
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
