Source Code for Module mvpa.clfs.glmnet

  1  # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- 
  2  # vi: set ft=python sts=4 ts=4 sw=4 et: 
  3  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  4  # 
  5  #   See COPYING file distributed along with the PyMVPA package for the 
  6  #   copyright and license terms. 
  7  # 
  8  ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## 
  9  """GLM-Net (GLMNET) regression classifier.""" 
 10   
 11  __docformat__ = 'restructuredtext' 
 12   
 13  # system imports 
 14  import numpy as N 
 15   
 16  import mvpa.base.externals as externals 
 17   
 18  # do conditional to be able to build module reference 
 19  if externals.exists('rpy', raiseException=True) and \ 
 20     externals.exists('glmnet', raiseException=True): 
 21      import rpy 
 22      rpy.r.library('glmnet') 
 23   
 24  # local imports 
 25  from mvpa.clfs.base import Classifier 
 26  from mvpa.measures.base import Sensitivity 
 27  from mvpa.misc.param import Parameter 
 28   
 29  if __debug__: 
 30      from mvpa.base import debug 
 31   
 32 -def _label2indlist(labels, ulabels): 
 33      """Convert labels to list of unique label indicies starting at 1. 
 34      """ 
 35   
 36      # allocate for the new one-of-M labels 
 37      new_labels = N.zeros(len(labels), dtype=N.int) 
 38   
 39      # loop and convert to one-of-M 
 40      for i, c in enumerate(ulabels): 
 41          new_labels[labels == c] = i+1 
 42   
 43      return [str(l) for l in new_labels.tolist()] 
 44   
 45   
 46 -class _GLMNET(Classifier): 
 47      """GLM-Net regression (GLMNET) `Classifier`. 
 48   
 49      GLM-Net is the model selection algorithm from: 
 50   
 51      Friedman, J., Hastie, T. and Tibshirani, R. (2008) Regularization 
 52      Paths for Generalized Linear Models via Coordinate 
 53      Descent. http://www-stat.stanford.edu/~hastie/Papers/glmnet.pdf 
 54   
 55      To make use of GLMNET, you must have R and RPy installed as well 
 56      as both the glmnet contributed package. You can install the R and 
 57      RPy with the following command on Debian-based machines: 
 58   
 59      sudo aptitude install python-rpy python-rpy-doc r-base-dev 
 60   
 61      You can then install the glmnet package by running R 
 62      as root and calling: 
 63   
 64      install.packages() 
 65   
 66      """ 
 67   
 68      _clf_internals = [ 'glmnet', 'linear', 'has_sensitivity', 
 69                         'does_feature_selection' 
 70                         ] 
 71   
 72      family = Parameter('gaussian', 
 73                         allowedtype='basestring', 
 74                         choices=["gaussian", "multinomial"], 
 75                         doc="""Response type of your labels 
 76                         (either regression or classification.""") 
 77   
 78      alpha = Parameter(1.0, min=0.01, max=1.0, allowedtype='float', 
 79                        doc="""The elastic net mixing parameter. 
 80                        Larger values will give rise to 
 81                        less L2 regularization, with alpha=1.0 
 82                        as a true lasso penalty.""") 
 83   
 84      nlambda = Parameter(100, allowedtype='int', min=1, 
 85                          doc="""Maximum number of lambdas to calculate 
 86                          before stopping if not converged.""") 
 87   
 88      standardize = Parameter(True, allowedtype='bool', 
 89                              doc="""Whether to standardize the variables 
 90                              prior to fitting.""") 
 91   
 92      thresh = Parameter(1e-4, min=1e-10, max=1.0, allowedtype='float', 
 93               doc="""Convergence threshold for coordinate descent.""") 
 94   
 95      pmax = Parameter(None, min=1, allowedtype='None or int', 
 96               doc="""Limit the maximum number of variables ever to be 
 97               nonzero.""") 
 98   
 99      maxit = Parameter(100, min=10, allowedtype='int', 
100               doc="""Maximum number of outer-loop iterations for 
101               'multinomial' families.""") 
102   
103      model_type = Parameter('covariance', 
104                             allowedtype='basestring', 
105                             choices=["covariance", "naive"], 
106               doc="""'covariance' saves all inner-products ever 
107               computed and can be much faster than 'naive'. The 
108               latter can be more efficient for 
109               nfeatures>>nsamples situations.""") 
110   
111 -    def __init__(self, **kwargs): 
112          """ 
113          Initialize GLM-Net. 
114   
115          See the help in R for further details on the parameters 
116          """ 
117          # init base class first 
118          Classifier.__init__(self, **kwargs) 
119   
120          # pylint friendly initializations 
121          self.__weights = None 
122          """The beta weights for each feature.""" 
123          self.__trained_model = None 
124          """The model object after training that will be used for 
125          predictions.""" 
126          self.__trained_model_dict = None 
127          """The model object in dict form after training that will be 
128          used for predictions.""" 
129   
130          # It does not make sense to calculate a confusion matrix for a 
131          # regression 
132          if self.params.family == 'gaussian': 
133              self.states.enable('training_confusion', False) 
134   
135  #     def __repr__(self): 
136  #         """String summary of the object 
137  #         """ 
138  #         return """ENET(lm=%s, normalize=%s, intercept=%s, trace=%s, max_steps=%s, enable_states=%s)""" % \ 
139  #                (self.__lm, 
140  #                 self.__normalize, 
141  #                 self.__intercept, 
142  #                 self.__trace, 
143  #                 self.__max_steps, 
144  #                 str(self.states.enabled)) 
145   
146   
147 -    def _train(self, dataset): 
148          """Train the classifier using `data` (`Dataset`). 
149          """ 
150          # process the labels based on the model family 
151          if self.params.family == 'gaussian': 
152              # do nothing, just save the labels as a list 
153              labels = dataset.labels.tolist() 
154              pass 
155          elif self.params.family == 'multinomial': 
156              # turn lables into list of range values starting at 1 
157              labels = _label2indlist(dataset.labels, 
158                                      dataset.uniquelabels) 
159          self.__ulabels = dataset.uniquelabels.copy() 
160   
161          # process the pmax 
162          if self.params.pmax is None: 
163              # set it to the num features 
164              pmax = dataset.nfeatures 
165          else: 
166              # use the value 
167              pmax = self.params.pmax 
168   
169          # train with specifying max_steps 
170          # must not convert trained model to dict or we'll get segfault 
171          rpy.set_default_mode(rpy.NO_CONVERSION) 
172          self.__trained_model = rpy.r.glmnet(dataset.samples, 
173                                              labels, 
174                                              family=self.params.family, 
175                                              alpha=self.params.alpha, 
176                                              nlambda=self.params.nlambda, 
177                                              standardize=self.params.standardize, 
178                                              thresh=self.params.thresh, 
179                                              pmax=pmax, 
180                                              maxit=self.params.maxit, 
181                                              type=self.params.model_type) 
182          rpy.set_default_mode(rpy.NO_DEFAULT) 
183   
184          # get a dict version of the model 
185          self.__trained_model_dict = rpy.r.as_list(self.__trained_model) 
186   
187          # save the lambda of last step 
188          self.__last_lambda = self.__trained_model_dict['lambda'][-1] 
189   
190          # set the weights to the last step 
191          weights = rpy.r.coef(self.__trained_model, s=self.__last_lambda) 
192          if self.params.family == 'multinomial': 
193              self.__weights = N.hstack([rpy.r.as_matrix(weights[str(i)])[1:] 
194                                         for i in range(1,len(self.__ulabels)+1)]) 
195          elif self.params.family == 'gaussian': 
196              self.__weights = rpy.r.as_matrix(weights)[1:] 
197   
198   
199 -    def _predict(self, data): 
200          """ 
201          Predict the output for the provided data. 
202          """ 
203          # predict with standard method 
204          values = rpy.r.predict(self.__trained_model, 
205                                 newx=data, 
206                                 type='link', 
207                                 s=self.__last_lambda) 
208   
209          # predict with the final state (i.e., the last step) 
210          classes = None 
211          if self.params.family == 'multinomial': 
212              # remove last dimension of values 
213              values = values[:,:,0] 
214   
215              # get the classes too (they are 1-indexed) 
216              rpy.set_default_mode(rpy.NO_CONVERSION) 
217              class_ind = rpy.r.predict(self.__trained_model, 
218                                        newx=data, 
219                                        type='class', 
220                                        s=self.__last_lambda) 
221              rpy.set_default_mode(rpy.NO_DEFAULT) 
222              class_ind = rpy.r.as_vector(class_ind) 
223   
224              # convert the strings to ints and subtract 1 
225              class_ind = N.array([int(float(c))-1 for c in class_ind]) 
226   
227              # convert to actual labels 
228              classes = self.__ulabels[class_ind] 
229          else: 
230              # is gaussian, so just remove last dim of values 
231              values = values[:,0] 
232   
233          if not classes is None: 
234              # set the values and return none 
235              self.values = values 
236              return classes 
237          else: 
238              # return the values as predictions 
239              return values 
240   
241   
242 -    def _getFeatureIds(self): 
243          """Return ids of the used features 
244          """ 
245          return N.where(N.abs(self.__weights)>0)[0] 
246   
247   
248   
249 -    def getSensitivityAnalyzer(self, **kwargs): 
250          """Returns a sensitivity analyzer for GLMNET.""" 
251          return GLMNETWeights(self, **kwargs) 
252   
253      weights = property(lambda self: self.__weights) 
254   
255   
256   
257 -class GLMNETWeights(Sensitivity): 
258      """`SensitivityAnalyzer` that reports the weights GLMNET trained 
259      on a given `Dataset`. 
260      """ 
261   
262      _LEGAL_CLFS = [ _GLMNET ] 
263   
264 -    def _call(self, dataset=None): 
265          """Extract weights from GLMNET classifier. 
266   
267          GLMNET always has weights available, so nothing has to be computed here. 
268          """ 
269          clf = self.clf 
270          weights = clf.weights 
271   
272          if __debug__: 
273              debug('GLMNET', 
274                    "Extracting weights for GLMNET - "+ 
275                    "Result: min=%f max=%f" %\ 
276                    (N.min(weights), N.max(weights))) 
277   
278          return weights 
279   
280 -class GLMNET_R(_GLMNET): 
281      """ 
282      GLM-NET Gaussian Regression Classifier. 
283   
284      This is the GLM-NET algorithm from 
285   
286      Friedman, J., Hastie, T. and Tibshirani, R. (2008) Regularization 
287      Paths for Generalized Linear Models via Coordinate 
288      Descent. http://www-stat.stanford.edu/~hastie/Papers/glmnet.pdf 
289   
290      parameterized to be a regression. 
291   
292      See GLMNET_C for the multinomial classifier version. 
293   
294      """ 
295   
296      _clf_internals = _GLMNET._clf_internals + ['regression'] 
297   
298 -    def __init__(self,  **kwargs): 
299          """ 
300          Initialize GLM-Net. 
301   
302          See the help in R for further details on the parameters 
303          """ 
304          # make sure they didn't specify regression 
305          if not kwargs.pop('family', None) is None: 
306              warning('You specified the "family" parameter, but we ' 
307                      'force this to be "gaussian".') 
308   
309          # init base class first, forcing regression 
310          _GLMNET.__init__(self, family='gaussian', **kwargs) 
311   
312   
313 -class GLMNET_C(_GLMNET): 
314      """ 
315      GLM-NET Multinomial Classifier. 
316   
317      This is the GLM-NET algorithm from 
318   
319      Friedman, J., Hastie, T. and Tibshirani, R. (2008) Regularization 
320      Paths for Generalized Linear Models via Coordinate 
321      Descent. http://www-stat.stanford.edu/~hastie/Papers/glmnet.pdf 
322   
323      parameterized to be a multinomial classifier. 
324   
325      See GLMNET_Class for the gaussian regression version. 
326   
327      """ 
328   
329      _clf_internals = _GLMNET._clf_internals + ['multiclass', 'binary'] 
330   
331 -    def __init__(self,  **kwargs): 
332          """ 
333          Initialize GLM-Net multinomial classifier. 
334   
335          See the help in R for further details on the parameters 
336          """ 
337          # make sure they didn't specify regression 
338          if not kwargs.pop('family', None) is None: 
339              warning('You specified the "family" parameter, but we ' 
340                      'force this to be "multinomial".') 
341   
342          # init base class first, forcing regression 
343          _GLMNET.__init__(self, family='multinomial', **kwargs) 
344