[R] Extracting P-values from the lrm function in the rms library

[David Winsemius]

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