Viitama Viitama - 2 months ago 7
R Question

How to export all coefficients of a penlized regression model from package `penalized`? Need them for reporting rolling regression estimate

I have been able to run regression with some coefficients constrained to positive territory, but I'm doing alot of rolling regressions where I face the problem. Here is my sample code:

library(penalized)
set.seed(1)
x1=rnorm(100)*10
x2=rnorm(100)*10
x3=rnorm(100)*10
y=sin(x1)+cos(x2)-x3+rnorm(100)

data <- data.frame(y, x1, x2, x3)

win <- 10

coefs <- matrix(NA, ncol=5, nrow=length(y))

for(i in 1:(length(y)-win)) {
d <- data[(1+i):(win+i),]
p <- win+i
# Linear Regression
coefs[p,] <- as.vector(coef(penalized(y, ~ x1 + x2 + x3, ~1,
lambda1=0, lambda2=0, positive = c(F, F, T), data=data)))}


This is how I usually populate matrix with coefs from rolling regression and now I receive error:

Error in coefs[p, ] <- as.vector(coef(penalized(y, ~x1 + x2 + x3, ~1, :
number of items to replace is not a multiple of replacement length


I assume that this error is produced because there is not always Intercept + 3 coefficients coming out of that penalized regression function. Is there away to get
penalized
function to show 0 coefs as well? or other way to populated matrix / data.frame?

Answer

You at most have 4 coefficients. (Seems like you know this by saying "intercept + 3 coefficients" but not sure why you in your code set up a matrix of 5 columns.)

Also, perhaps you are unaware of the which argument for coef for "penfit" object. Have a look at:

getMethod(coef, "penfit")

#function (object, ...) 
#{
#    .local <- function (object, which = c("nonzero", "all", "penalized", 
#        "unpenalized"), standardize = FALSE) 
#    {
#        coefficients(object, which, standardize)
#    }
#    .local(object, ...)
#}
#<environment: namespace:penalized>

We can set which = "all" to report all coefficients. The default is which = "nonzero", which is causing the "replacement length differs" issue.

The following works:

library(penalized)
set.seed(1)
x1 = rnorm(100)*10
x2 = rnorm(100)*10
x3 = rnorm(100)*10
y = sin(x1) + cos(x2) - x3 + rnorm(100)

data <- data.frame(y, x1, x2, x3)

win <- 10

coefs <- matrix(NA, ncol=4, nrow=length(y))

for(i in 1:(length(y)-win)) {
  d <- data[(1+i):(win+i),]
  p <- win + i
  pen <- penalized(y, ~ x1 + x2 + x3, ~1, lambda1 = 0, lambda2 = 0,
                   positive = c(F, F, T), data = data)
  beta <- coef(pen, which = "all")
  coefs[p,] <- unname(beta)
  }