selecting penalty fractions

[yonicd]

from discuss forum: http://disq.us/p/2a36lz7

Hello, I have a question. To use 'lmmen' function, I should decide penalty level for fixed and random effects first. How can I find the most proper fraction? I guess golden_section and golden_section_2d function in your packages are useful to do it, I can't find how to connect between two. I really want to use your packages, please help me. Thanks.

[yonicd]

thanks for trying out the package!

You are correct in using the golden* functions.

They are wrappers on the main package function lmmen

The user is intended to supply the data, initial fixed effects and which penalty to search on.

There is an example in each of the functions' documentation

dat <- initialize_example(n.i = 5,n = 30,q=4,seed=1)
init <- init.beta(dat,method='glmnet')
golden_section(dat,init,pen.effect = 'FE.L1')

e.g. you can see the calls to lmmen here: https://github.com/yonicd/lmmen

• golden_section

• golden_section_2d

[GUMIRYEONG]

I appreciate your help :) I think I can get lower BIC value when I use golden_section_2d function, because it can consider not only L1 penalty and L2 penalty simultaneously. Am I right? In golden_section_2d, argument : l2 (L2 fixed and random effect) default values are c(1,1). If I set L2 penalty as default, can I get proper result? Any big differences?

[yonicd]

Yes.

The 2d should give lower bic, but is much more computationally intensive.

[yonicd]

regarding L2.

Much like a regular EN CV over ridge levels you would still need to CV across a range of L2 [0,1]x[0,1].

The farther away from 1 the more grouping you are allowing across the variables.

[GUMIRYEONG]

Thanks a lot :) Tuning L2 penalty with L1 penalty sounds good but I think it is computationally intensive. Is it okay to tune only L1 letting L2 penalty on default value? I think it is similar situation in elastic net, we can just make one of tuning parameter, alpha, which can decide how much we will consider L1 and L2 penalty 0.5.

[GUMIRYEONG]

And how can you set "opt.tau"? Default is (sqrt(5)-1)/2). I was trying to analyze using your package using default value, I got error message:

Error in [<-(*tmp*, start.point:end.point, (q (i - 1) + 1):(q : subscript out of bounds

I think it is because initial value of gs.x is smaller than 0. gs.x[1] = -0.73, gs.x[2]=1.73 actually. How can I change "opt.tau"?

[yonicd]

Changing this value would render the function to not use the golden section ratio, which is the intention of the functions. Based on this paper:

J. Kiefer. Sequential minimax search for a maximum. Proceedings of the American mathematical society, 4(3):502–506, 1953

If you want to change it opt.tau is an argument of the golden_section* functions.

see here.

[yonicd]

Regarding alpha penalty (ridge (0) to lasso (1)) selection:

In lmmen it is not set directly by the user.

  # alpha = l1.f / (l1.f + l2.f)
  # l1.f = (alpha * l1.f) + (alpha * l2.f)
  # l2.f = (l1.f - (alpha * l1.f)) / alpha

[GUMIRYEONG]

Thanks! One more question. I want assess model fit. To do that, I should predict by using test data set, I can't find that kinds of function in your package. So I tried to use general predict function, but I failed. Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "lmmen" Do you have any ideas about it? I really want to find solutions.

[yonicd]

unfortunately i was a young developer when i wrote this package and never wrote a predict method.

I wrote that manually for the accompanying paper.

You are welcome to create a Pull Request if you write one before I have time again to dive into this package again.

[GUMIRYEONG]

Oh, I see. thank you. According to your paper, you compared the probabilistic forecast accuracy of the two methods using the Brier Score. I'm not good at programming, so I'd like to calculate prediction error manually. To do that, I need estimated random effects. Through lmmen, I can get estimated random effects covariance matrix standard deviations. Could you please explain how to get estimated random effects ? Thank you so much !

[yonicd]

I dont have the bandwidth at the moment to write up the full code, but the general idea is that you would want to solve for equation (2) in the paper.

Given the assumed distribution of b_i ~ N(0,hat{D}), ie you could use the estimate of D from the output to draw b_i.

The output of lmmen gives all the building blocks you would need for this calculation.

see: https://github.com/yonicd/lmmen/blob/master/R/lmmen.R#L332_L334

Hope this helps

[GUMIRYEONG]

I see. After getting hat{D}, I can generate bi for the multivariate normal distribution with mean 0 and covariance matrix hat{D} randomly. Am I right? Can I one more question? I think 'lmmen' function is optimized in case of random effects model. I want use that function in case of random intercept model. In this case, random effect is only one. How can I change lmmen function code? Ask for your help, thanks!