Gotmadstacks Gotmadstacks - 7 days ago 6
R Question

Coefficient table does not have NA rows in rank-deficient fit; how to insert them?

library(lmPerm)
x <- lmp(formula = a ~ b * c + d + e, data = df, perm = "Prob")

summary(x) # truncated output, I can see `NA` rows here!

#Coefficients: (1 not defined because of singularities)
# Estimate Iter Pr(Prob)
#b 5.874 51 1.000
#c -30.060 281 0.263
#b:c NA NA NA
#d1 -31.333 60 0.633
#d2 33.297 165 0.382
#d3 -19.096 51 1.000
#e 1.976 NA NA


I want to pull out the
Pr(Prob)
results for everything, but

y <- summary(x)$coef[, "Pr(Prob)"]

#(Intercept) b c d1 d2
# 0.09459459 1.00000000 0.26334520 0.63333333 0.38181818
# d3 e
# 1.00000000 NA


This is not what I want. I need
b:c
row, too, in the right position.


An example of the output I would like from the above would be:

# (Intercept) b c b:c d1 d2
# 0.09459459 1.00000000 0.26334520 NA 0.63333333 0.38181818
# d3 e
# 1.00000000 NA


I also would like to pull out the
Iter
column that corresponds to each variable. Thanks.

Answer

lmp is based on lm and summary.lmp also behaves like summary.lm, so I will first use lm for illustration, then show that we can do the same for lmp.


lm and summary.lm

Have a read on ?summary.lm and watch out for the following returned values:

coefficients: a p x 4 matrix with columns for the estimated
              coefficient, its standard error, t-statistic and
              corresponding (two-sided) p-value.  Aliased coefficients are
              omitted.

     aliased: named logical vector showing if the original coefficients are
              aliased.

When you have rank-deficient models, NA coefficients are omitted in the coefficient table, and they are called aliased variables. Consider the following small, reproducible example:

set.seed(0)
zz <- xx <- rnorm(10)
yy <- rnorm(10)
fit <- lm(yy ~ xx + zz)

coef(fit)  ## we can see `NA` here
#(Intercept)          xx          zz 
#  0.1295147   0.2706560          NA 

a <- summary(fit)  ## it is also printed to screen
#Coefficients: (1 not defined because of singularities)
#            Estimate Std. Error t value Pr(>|t|)
#(Intercept)   0.1295     0.3143   0.412    0.691
#xx            0.2707     0.2669   1.014    0.340
#zz                NA         NA      NA       NA

b <- coef(a)  ## but no `NA` returned in the matrix / table
#             Estimate Std. Error   t value  Pr(>|t|)
#(Intercept) 0.1295147  0.3142758 0.4121051 0.6910837
#xx          0.2706560  0.2669118 1.0140279 0.3402525

d <- a$aliased
#(Intercept)          xx          zz 
#      FALSE       FALSE        TRUE 

If you want to pad NA rows to coefficient table / matrix, we can do

## an augmented matrix of `NA`
e <- matrix(nrow = length(d), ncol = ncol(b),
            dimnames = list(names(d), dimnames(b)[[2]]))
## fill rows for non-aliased variables
e[!d] <- b

#             Estimate Std. Error   t value  Pr(>|t|)
#(Intercept) 0.1295147  0.3142758 0.4121051 0.6910837
#xx          0.2706560  0.2669118 1.0140279 0.3402525
#zz                 NA         NA        NA        NA

lmp and summary.lmp

Nothing needs be changed.

library(lmPerm)
fit <- lmp(yy ~ xx + zz, perm = "Prob")
a <- summary(fit)  ## `summary.lmp`
b <- coef(a)

#              Estimate Iter  Pr(Prob)
#(Intercept) -0.0264354  241 0.2946058
#xx           0.2706560  241 0.2946058

d <- a$aliased
#(Intercept)          xx          zz 
#      FALSE       FALSE        TRUE 

e <- matrix(nrow = length(d), ncol = ncol(b),
            dimnames = list(names(d), dimnames(b)[[2]]))
e[!d] <- b

#              Estimate Iter  Pr(Prob)
#(Intercept) -0.0264354  241 0.2946058
#xx           0.2706560  241 0.2946058
#zz                  NA   NA        NA

If you, want to extract Iter and Pr(Prob), just do

e[, 2]  ## e[, "Iter"]
#(Intercept)          xx          zz 
#        241         241          NA 

e[, 3]  ## e[, "Pr(Prob)"]
#(Intercept)          xx          zz 
#  0.2946058   0.2946058          NA 
Comments