Multinomial and Cox redux; also fixing showarray error

Project.toml

@@ -3,6 +3,7 @@ uuid = "8d5ece8b-de18-5317-b113-243142960cc6"
 version = "0.5.2"
 
 [deps]
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

README.md

@@ -7,7 +7,7 @@
 
 ## Quick start
 
-To fit a basic model:
+To fit a basic regression model:
 
 ```julia
 julia> using GLMNet
@@ -17,63 +17,137 @@ julia> y = collect(1:100) + randn(100)*10;
 julia> X = [1:100 (1:100)+randn(100)*5 (1:100)+randn(100)*10 (1:100)+randn(100)*20];
 
 julia> path = glmnet(X, y)
-Least Squares GLMNet Solution Path (55 solutions for 4 predictors in 163 passes):
-55x3 DataFrame:
-         df  pct_dev        λ
-[1,]      0      0.0  27.1988
-[2,]      1 0.154843  24.7825
-[3,]      1 0.283396  22.5809
-  :
-[53,]     2 0.911956 0.215546
-[54,]     2 0.911966 0.196397
-[55,]     2 0.911974  0.17895
+Least Squares GLMNet Solution Path (74 solutions for 4 predictors in 832 passes):
+74x3 DataFrame
+| Row | df | pct_dev  | λ         |
+|-----|----|----------|-----------|
+| 1   | 0  | 0.0      | 29.6202   |
+| 2   | 1  | 0.148535 | 26.9888   |
+| 3   | 1  | 0.271851 | 24.5912   |
+| 4   | 1  | 0.37423  | 22.4066   |
+⋮
+| 70  | 4  | 0.882033 | 0.0482735 |
+| 71  | 4  | 0.882046 | 0.043985  |
+| 72  | 4  | 0.882058 | 0.0400775 |
+| 73  | 4  | 0.882067 | 0.0365171 |
+| 74  | 4  | 0.882075 | 0.033273  |
 ```
 
 `path` represents the Lasso or ElasticNet fits for varying values of λ. The value of the intercept for each λ value are in `path.a0`. The coefficients for each fit are stored in compressed form in `path.betas`.
 
 ```julia
 julia> path.betas
-4x55 CompressedPredictorMatrix:
- 0.0  0.083706  0.159976  0.22947  …  0.929157    0.929315
- 0.0  0.0       0.0       0.0         0.00655753  0.00700862
- 0.0  0.0       0.0       0.0         0.0         0.0
- 0.0  0.0       0.0       0.0         0.0         0.0
+4x74 CompressedPredictorMatrix:
+ 0.0  0.091158  0.174218  …   0.913497   0.915593   0.917647
+ 0.0  0.0       0.0           0.128054   0.127805   0.127568
+ 0.0  0.0       0.0          -0.126211  -0.128015  -0.129776
+ 0.0  0.0       0.0           0.108217   0.108254   0.108272
 ```
 
 This CompressedPredictorMatrix can be indexed as any other AbstractMatrix, or converted to a Matrix using `convert(Matrix, path.betas)`.
 
+One can visualize the path by
+
+```julia
+julia> using Gadfly
+
+julia> plot(path, Guide.xlabel("||β||₁"), Guide.ylabel("βᵢ"), x=:norm1)
+```
+![regression-lasso-path](https://rawgit.com/linxihui/Misc/master/Images/GLMNet.jl/regression_lasso_path.svg)
+
+One can see that the LASSO path is piecewise linear.
+
 To predict the output for each model along the path for a given set of predictors, use `predict`:
 
 ```julia
 julia> predict(path, [22 22+randn()*5 22+randn()*10 22+randn()*20])
-1x55 Array{Float64,2}:
- 51.7098  49.3242  47.1505  45.1699  …  25.1036  25.0878  25.0736
+1x74 Array{Float64,2}:
+ 50.8669  48.2689  45.9017  …  21.9344  21.9377  21.9407
 ```
 
 To find the best value of λ by cross-validation, use `glmnetcv`:
 
 ```julia
 julia> cv = glmnetcv(X, y)
 Least Squares GLMNet Cross Validation
-55 models for 4 predictors in 10 folds
-Best λ 0.343 (mean loss 76.946, std 12.546)
+74 models for 4 predictors in 10 folds
+Best λ 0.450 (mean loss 129.720, std 14.871)
 
 julia> argmin(cv.meanloss)
 48
 
-julia> cv.path.betas[:, 48]
+julia> cv.path.betas[:, 46]
 4-element Array{Float64,1}:
- 0.926911
- 0.00366805
+ 0.781119
+ 0.128094
  0.0
+ 0.103008
+
+julia> coef(cv)
+4-element Array{Float64,1}:
+ 0.781119
+ 0.128094
  0.0
+ 0.103008
+```
+
+### A classification Example
+
+```julia
+julia> using RDatasets
+
+julia> iris = dataset("datasets", "iris");
+
+julia> X = convert(Matrix, iris[:, 1:4]);
+
+julia> y = convert(Vector, iris[:Species]);
+
+julia> iTrain = sample(1:size(X,1), 100, replace = false);
+
+julia> iTest = setdiff(1:size(X,1), iTrain);
+
+julia> iris_cv = glmnetcv(X[iTrain, :], y[iTrain])
+Multinomial GLMNet Cross Validation
+100 models for 4 predictors in 10 folds
+Best λ 0.001 (mean loss 0.124, std 0.046)
+
+julia> yht = round(predict(iris_cv, X[iTest, :], outtype = :prob), 3);
+
+julia> DataFrame(target=y[iTest], set=yht[:,1], ver=yht[:,2], vir=yht[:,3])[5:5:50,:]
+10x4 DataFrame
+| Row | target       | set   | ver   | vir   |
+|-----|--------------|-------|-------|-------|
+| 1   | "setosa"     | 0.999 | 0.001 | 0.0   |
+| 2   | "setosa"     | 1.0   | 0.0   | 0.0   |
+| 3   | "setosa"     | 1.0   | 0.0   | 0.0   |
+| 4   | "versicolor" | 0.0   | 0.983 | 0.017 |
+| 5   | "versicolor" | 0.002 | 0.961 | 0.037 |
+| 6   | "versicolor" | 0.0   | 0.067 | 0.933 |
+| 7   | "versicolor" | 0.0   | 0.993 | 0.007 |
+| 8   | "virginica"  | 0.0   | 0.0   | 1.0   |
+| 9   | "virginica"  | 0.0   | 0.397 | 0.603 |
+| 10  | "virginica"  | 0.0   | 0.025 | 0.975 |
+
+julia> plot(iris_cv.path, Scale.x_log10, Guide.xlabel("λ"), Guide.ylabel("βᵢ"))
+```
+![iris-lasso-path](https://rawgit.com/linxihui/Misc/master/Images/GLMNet.jl/iris_lasso_path.svg) 
+
+```julia
+julia> plot(iris_cv)
 ```
+![iris-cv](https://rawgit.com/linxihui/Misc/master/Images/GLMNet.jl/iris_cv.svg)
 
 ## Fitting models
 
-`glmnet` has two required parameters: the m x n predictor matrix `X` and the dependent variable `y`. It additionally accepts an optional third argument, `family`, which can be used to specify a generalized linear model. Currently, only `Normal()` (least squares, default), `Binomial()` (logistic), and `Poisson()` are supported, although the glmnet Fortran code also implements a Cox model. For logistic models, `y` is a m x 2 matrix, where the first column is the count of negative responses for each row in `X` and the second column is the count of positive responses. For all other models, `y` is a vector.
+`glmnet` has two required parameters: the m x n predictor matrix `X` and the dependent variable `y`. It additionally accepts an optional third argument, `family`, which can be used to specify a generalized linear model. Currently, `Normal()` (least squares, default), `Binomial()` (logistic), `Poisson()` , `Multinomial()`, `CoxPH()` (Cox model) are supported. 
+
+- For linear and Poisson models, `y` is a numerical vector.
+- For logistic models, `y` is either a string vector or a m x 2 matrix, where the first column is the count of negative responses for each row in `X` and the second column is the count of positive responses. 
+- For multinomial models, `y` is etiher a string vector (with at least 3 unique values) or a m x k matrix, where k is number of unique values (classes).
+- For Cox models, `y` is a 2-column matrix, where the first column is survival time and second column is (right) censoring status. Indeed, For survival data, `glmnet` has another method `glmnet(X::Matrix, time::Vector, status::Vector)`. Same for `glmnetcv`.
+
 
-`glmnet` also accepts many optional parameters, described below:
+`glmnet` also accepts many optional keyword parameters, described below:
 
  - `weights`: A vector of weights for each sample of the same size as `y`.
  - `alpha`: The tradeoff between lasso and ridge regression. This defaults to `1.0`, which specifies a lasso model.