phildeutsch - 1 year ago 112
R Question

# R package `penalized`: error when using predict() with matrix input for "penalized" and "unpenalized"

I am trying to use the penalized package to calibrate a penalized linear regression such that the coefficients of a subset of the variables are positive. I managed to calibrate the model but I am failing to use it for making new predictions.

Here's a toy example:

``````require(dplyr)
require(penalized)
require(ggplot2)

# A tibble: 6 x 10
#  carat       cut color clarity depth table price     x     y     z
#  <dbl>     <ord> <ord>   <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
#1  0.23     Ideal     E     SI2  61.5    55   326  3.95  3.98  2.43
#2  0.21   Premium     E     SI1  59.8    61   326  3.89  3.84  2.31
#3  0.23      Good     E     VS1  56.9    65   327  4.05  4.07  2.31
#4  0.29   Premium     I     VS2  62.4    58   334  4.20  4.23  2.63
#5  0.31      Good     J     SI2  63.3    58   335  4.34  4.35  2.75
#6  0.24 Very Good     J    VVS2  62.8    57   336  3.94  3.96  2.48

response = diamonds\$price
penalized_vars = "x"
unpenalized_vars = "depth"
fit_penalized = penalized(response=as.matrix(response),
penalized = model.matrix(~., select_(diamonds, .dots = penalized_vars)),
unpenalized = model.matrix(~., select_(diamonds, .dots = unpenalized_vars)),
model="linear",
positive=TRUE,
maxiter=25)
# nonzero coefficients: 3

show(fit_penalized)
#Penalized linear regression object
#4 regression coefficients of which 3 are non-zero

#Loglikelihood =     -482648.2

#         1          2          3          4          5          6
#-1679.6983 -1924.0012 -1515.2681  -863.6912  -393.7956 -1668.7105
``````

So far so good. How can I actually use this to predict values based on new information? I have tried

``````predict(fit_penalized,
penalized = model.matrix(~., select_(vars, .dots = penalized_vars)),
unpenalized = model.matrix(~., select_(vars, .dots = unpenalized_vars)) )
# Error in terms.default(object@formula\$unpenalized) :
#   no terms component nor attribute
``````

Regarding the latest `penalized_0.9-47` updated on 2016-05-27

As will be seen from my test below, `unpenalized` can only be specified via "formula" not "matrix", in order to get `predict` work. This seems to be a bug in the package. The bug is not with `penalized` but with `predict`. On one hand, `?penalized::predict` says:

``````In particular, if penalized and/or unpenalized
was specified in matrix form, a matrix must be given with the new
subjects' data. The columns of these matrices must be exactly the
same as in the matrices supplied in the original call that
produced the ‘penfit’ object. If either penalized or unpenalized
was given as a ‘formula’ in the original call, the user of
‘predict’ must supply a new ‘data’ argument.
``````

It seems legitimate to pass matrix to `unpenalized`, but my practical test denies this. So try contact package author.

As a side note, your current specification will have a penalized intercept, as well as a free intercept. You should drop one of them for identifiability. In the following, I drop the free intercept.

``````library(penalized)
library(ggplot2)
X1 <- model.matrix (~x, diamonds)    ## model matrix for penalized terms
X2 <- model.matrix(~ depth - 1, diamonds)    ## model matrix for free terms
vars <- diamonds[1:5, ]    ## prediction dataset
Xp1 <- model.matrix(~x, vars)    ## prediction matrix for penalized terms
Xp2 <- model.matrix(~ depth - 1, vars)    ## prediction matrix for free terms

## use "formula" for both
fit <- penalized (price, ~ x, ~ depth - 1, data = diamonds, model = "linear", positive = TRUE)
predict(fit, ~ x, ~ depth - 1, vars)
#          mu  sigma2
#1 -1523.5643 3600667
#2 -1328.9947 3600667
#3  -183.9204 3600667
#4  -950.1526 3600667
#5  -717.6854 3600667

## "matrix" for `penalized` and "formula" for `unpenalized`
fit <- penalized (price, X1, ~ depth - 1, data = diamonds, model = "linear", positive = TRUE)
predict(fit, Xp1, ~ depth - 1, vars)
#          mu  sigma2
#1 -1523.5643 3600667
#2 -1328.9947 3600667
#3  -183.9204 3600667
#4  -950.1526 3600667
#5  -717.6854 3600667

## "formula" for `penalized` and "matrix" for `unpenalized`
fit <- penalized (price, ~ x, X2, data = diamonds, model = "linear", positive = TRUE)
predict(fit, ~ x, Xp2, vars)
# Error in terms.default(object@formula\$unpenalized) :
#  no terms component nor attribute

## "matrix" for both
fit <- penalized (price, X1, X2, data = diamonds, model = "linear", positive = TRUE)
predict(fit, Xp1, Xp2, vars)
# Error in terms.default(object@formula\$unpenalized) :
#  no terms component nor attribute
``````
Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download