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

[Frank E Harrell Jr]

On 06/19/2010 10:30 AM, David Winsemius wrote:
>
> 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"
>

For once all the solutions offered are incorrect.  They only work in the 
special case where each variable has only one degree of freedom.

Just do a <- anova(mod.poly3) and treat the result as matrix.  You'll 
get the needed multiple degree of freedom test and P-value.

Frank
-- 
Frank E Harrell Jr   Professor and Chairman        School of Medicine
                      Department of Biostatistics   Vanderbilt University