Merge pull request #8 from dmbates/zygote

Use Zygote for gradient and lowrankupdate! for least squares solution

Package checks are failing because this package is not yet registered and therefore can't be installed for testing.

.appveyor.yml

@@ -0,0 +1,31 @@
+# Documentation: https://github.com/JuliaCI/Appveyor.jl
+environment:
+  matrix:
+  - julia_version: 1.2
+  - julia_version: nightly
+platform:
+  - x86
+  - x64
+matrix:
+  allow_failures:
+    - julia_version: nightly
+branches:
+  only:
+    - master
+    - /release-.*/
+notifications:
+  - provider: Email
+    on_build_success: false
+    on_build_failure: false
+    on_build_status_changed: false
+install:
+  - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1"))
+build_script:
+  - echo "%JL_BUILD_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%"
+test_script:
+  - echo "%JL_TEST_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%"
+on_success:
+  - echo "%JL_CODECOV_SCRIPT%"
+  - C:\julia\bin\julia -e "%JL_CODECOV_SCRIPT%"

.gitignore

@@ -0,0 +1,8 @@
+*.jl.*.cov
+*.jl.cov
+*.jl.mem
+.DS_Store
+/Manifest.toml
+/dev/
+/docs/build/
+/docs/site/

.travis.yml

@@ -1,13 +1,26 @@
-language: cpp
-compiler:
-  - clang
+# Documentation: http://docs.travis-ci.com/user/languages/julia/
+language: julia
+os:
+  - linux
+  - osx
+julia:
+  - 1.2
+  - nightly
+matrix:
+  allow_failures:
+    - julia: nightly
+  fast_finish: true
 notifications:
   email: false
-before_install:
-  - sudo add-apt-repository ppa:staticfloat/julia-deps -y
-  - sudo add-apt-repository ppa:staticfloat/julianightlies -y
-  - sudo apt-get update -qq -y
-  - sudo apt-get install libpcre3-dev julia -y
-script:
-  - julia -e 'Pkg.init(); run(`ln -s $(pwd()) $(Pkg.dir())/NLreg`); Pkg.resolve()'
-  - julia -e 'using NLreg; @assert isdefined(:NLreg); @assert typeof(NLreg) === Module'
+after_success:
+  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'
+jobs:
+  include:
+    - stage: Documentation
+      julia: 1.2
+      script: julia --project=docs -e '
+          using Pkg;
+          Pkg.develop(PackageSpec(path=pwd()));
+          Pkg.instantiate();
+          include("docs/make.jl");'
+      after_success: skip

Project.toml

@@ -1,25 +1,17 @@
 name = "NLreg"
 uuid = "6aa54777-d00a-57a2-a775-234c624c12d3"
 authors = ["Douglas Bates <[email protected]>"]
+version = "0.2.0"
 
 [deps]
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
-DataFrames = "0.18"
-Distributions = "0.19, 0.20"
-ForwardDiff = "0.10"
-StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30"
+Distributions = "0.19, 0.20, 0.21"
+StatsBase = "0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31"
+Zygote = "0.3"
 julia = "1"
-
-[extras]
-CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[targets]
-test = ["CSV", "Test"]

README.md

@@ -1,50 +1,14 @@
-# Nonlinear regression models in [Julia](http://julialang.org)
+# Nonlinear regression
 
-[![Build Status](https://travis-ci.org/dmbates/NLreg.jl.png)](https://travis-ci.org/dmbates/NLreg.jl)
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://dmbates.github.io/NLreg.jl/stable)
+[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://dmbates.github.io/NLreg.jl/dev)
+[![Build Status](https://travis-ci.com/dmbates/NLreg.jl.svg?branch=master)](https://travis-ci.com/dmbates/NLreg.jl)
+[![Build Status](https://ci.appveyor.com/api/projects/status/github/dmbates/NLreg.jl?svg=true)](https://ci.appveyor.com/project/dmbates/NLreg-jl)
+[![Codecov](https://codecov.io/gh/dmbates/NLreg.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/dmbates/NLreg.jl)
 
-In this [Julia](http://julialang.org) package nonlinear regression models are formulated as Julia types that inherit from `NLregMod`.
-A simple example is the predicted concentration in a 1 compartment model with a single bolus dose at time 0.
-```jl
-conc = exp(logV) * exp(-exp(logK)*t)
-```
-where `logV` and `logK` are the logarithms of the volume of distribution and `K` is the elimination rate constant and `t` is the time of the measurement.
-
-The `logsd1` type represents this model and the data to which it is to be fit.
-The fields of this type are `t`, the vector of times at which samples are drawn, `y`, the vector of measured concentrations, `mu` the mean responses at the current parameter values, `resid` the residuals at the current parameter values, and `tgrad`, the transpose of the gradient matrix.
-The external constructors for this model allow it to be specified from `t` and `y` or in a Formula/Data specification.
-
-A nonlinear regression model must provide methods for `pnames`, the parameter names, `updtmu`, update the mean response, residuals and `tgrad` from new parameter values, and `initpars`, determine initial parameter estimates from the data.
-
-```jl
-julia> using DataFrames, NLreg
-
-julia> const sd1 = readtable(Pkg.dir("NLreg","data","sd1.csv.gz"));
-
-julia> nl = fit(BolusSD1(CONC ~ TIME, sd1))
-Nonlinear least squares fit to 580 observations
-
-     Estimate Std.Error t value Pr(>|t|)
-V     1.14296 0.0495656 23.0595  < eps()
-K    0.245688 0.0202414 12.1379  < eps()
-
-Residual sum of squares at estimates: 110.597
-Residual standard error = 0.43743 on 578 degrees of freedom
-```
-
-## Plans for the near future
+This package is an experiment in using the [`Zygote`](https://github.com/FluxML/Zygote.jl) automatic differentiation package and the `lowrankupdate!` function in the `LinearAlgebra` package to solve the linear least squares problem for a Gauss-Newton update.
 
-- Nonlinear mixed-effects models fit using the Laplace approximation to the log-likelihood
-
-- Specification of partially linear models
-
-- Composite models consisting of a parameter transformation and a nonlinear model.
-
-## Partially linear models (this used to work but is now broken)
-
-Partially linear models (those models with some parameters that occur
-linearly in the model expression) are expressed as types that inherit
-from the `PLregMod` abstract type.  A instance of a model type is
-created from the values of any covariates in the model.
+The data are represented as a `Tables.RowTable`, which is a vector of `NamedTuple`s.  The model parameters are also a `NamedTuple`.  The model function is given as a function of two arguments - the parameters and a data row.
 
 ## Example - a Michaelis-Menten fit
 
@@ -64,124 +28,44 @@ To fit such a model we create a `MicMen` object from the vector of
 observed concentrations and a `PLregFit` object from this model and
 the responses.
 ```julia
-julia> pur = dataset("datasets","Puromycin");
-
-julia> purtrt = sub(pur, pur[:State] .== "treated")
-12x3 SubDataFrame{Array{Int64,1}}
-|-------|------|------|-----------|
-| Row # | Conc | Rate | State     |
-| 1     | 0.02 | 76   | "treated" |
-| 2     | 0.02 | 47   | "treated" |
-| 3     | 0.06 | 97   | "treated" |
-| 4     | 0.06 | 107  | "treated" |
-| 5     | 0.11 | 123  | "treated" |
-| 6     | 0.11 | 139  | "treated" |
-| 7     | 0.22 | 159  | "treated" |
-| 8     | 0.22 | 152  | "treated" |
-| 9     | 0.56 | 191  | "treated" |
-| 10    | 0.56 | 201  | "treated" |
-| 11    | 1.1  | 207  | "treated" |
-| 12    | 1.1  | 200  | "treated" |
-
-julia> pm1 = fit(MicMen(Rate ~ Conc, purtrt), true)
-Iteration:  1, rss = 1679.58, cvg = 0.257787 at [201.837,0.0484065]
-   Incr: [0.012547]  f = 1.0, rss = 1211.92
-Iteration:  2, rss = 1211.92, cvg = 0.0122645 at [210.623,0.0609536]
-   Incr: [0.00280048]  f = 1.0, rss = 1195.66
-Iteration:  3, rss = 1195.66, cvg = 0.000161307 at [212.448,0.063754]
-   Incr: [0.000331041]  f = 1.0, rss = 1195.45
-Iteration:  4, rss = 1195.45, cvg = 1.56158e-6 at [212.661,0.0640851]
-   Incr: [3.27084e-5]  f = 1.0, rss = 1195.45
-Iteration:  5, rss = 1195.45, cvg = 1.4557e-8 at [212.681,0.0641178]
-   Incr: [3.15934e-6]  f = 1.0, rss = 1195.45
-Iteration:  6, rss = 1195.45, cvg = 1.35192e-10 at [212.684,0.0641209]
-   Incr: [3.04477e-7]  f = 1.0, rss = 1195.45
-
-Nonlinear least squares fit to 12 observations
-
-      Estimate  Std.Error t value Pr(>|t|)
-Vm     212.684    6.94715 30.6145  3.2e-11
-K    0.0641212 0.00828095 7.74323   1.6e-5
-
-Residual sum of squares at estimates: 1195.45
-Residual standard error = 10.9337 on 10 degrees of freedom
-julia> using RDatasets, NLreg
-
-julia> purtrt = sub(dataset("datasets","Puromycin"),:(State .== "treated"));
-
-julia> pm1 = fit(MicMen(Conc ~ Rate, purtrt),true)
-Iteration:  1, rss = 0.188744, cvg = 0.0888416 at [-0.0786133,-220.728]
-   Incr: [-4.66305]  f = 1.0, rss = 0.173277
-Iteration:  2, rss = 0.173277, cvg = 0.00102418 at [-0.0995101,-225.391]
-   Incr: [-0.688546]  f = 1.0, rss = 0.173117
-Iteration:  3, rss = 0.173117, cvg = 6.54049e-6 at [-0.10249,-226.08]
-   Incr: [0.0574836]  f = 1.0, rss = 0.173116
-Iteration:  4, rss = 0.173116, cvg = 7.53229e-8 at [-0.102242,-226.022]
-   Incr: [-0.00614653]  f = 1.0, rss = 0.173116
-Iteration:  5, rss = 0.173116, cvg = 8.30647e-10 at [-0.102269,-226.028]
-   Incr: [0.000645718]  f = 1.0, rss = 0.173116
-
-Nonlinear least squares fit to 12 observations
-
-      Estimate Std.Error  t value Pr(>|t|)
-Vm   -0.102266 0.0315309 -3.24335   0.0088
-K     -226.028   7.08463 -31.9039  2.2e-11
-
-Residual sum of squares at estimates: 0.173116
-Residual standard error = 0.131574 on 10 degrees of freedom
-```
-
-We can also use parameter transformations
-
-```julia
-julia> pm2 = fit([LogTr] * MicMen(Rate ~ Conc, purtrt), true)
-Iteration:  1, rss = 1679.58, cvg = 0.257787 at [201.837,-3.02812]
-   Incr: [0.259201]  f = 1.0, rss = 1198.55
-Iteration:  2, rss = 1198.55, cvg = 0.00234006 at [211.785,-2.76892]
-   Incr: [0.0198568]  f = 1.0, rss = 1195.48
-Iteration:  3, rss = 1195.48, cvg = 2.12492e-5 at [212.598,-2.74906]
-   Incr: [0.00188326]  f = 1.0, rss = 1195.45
-Iteration:  4, rss = 1195.45, cvg = 1.96812e-7 at [212.675,-2.74718]
-   Incr: [0.000181183]  f = 1.0, rss = 1195.45
-Iteration:  5, rss = 1195.45, cvg = 1.82666e-9 at [212.683,-2.747]
-   Incr: [1.74545e-5]  f = 1.0, rss = 1195.45
-
-Nonlinear least squares fit to 12 observations
-
-        Estimate Std.Error  t value Pr(>|t|)
-Vm       212.684   6.94715  30.6145  3.2e-11
-log(K)  -2.74698  0.129145 -21.2705   1.2e-9
-
-Residual sum of squares at estimates: 1195.45
-Residual standard error = 10.9337 on 10 degrees of freedom
+julia> using CSV, DataFrames, NLreg
+
+julia> datadir = normpath(joinpath(dirname(pathof(NLreg)), "..", "data"));
+
+julia> PurTrt = first(groupby(CSV.read(joinpath(datadir, "Puromycin.csv")), :state))
+12×3 SubDataFrame
+│ Row │ conc    │ rate    │ state   │
+│     │ Float64 │ Float64 │ String  │
+├─────┼─────────┼─────────┼─────────┤
+│ 1   │ 0.02    │ 76.0    │ treated │
+│ 2   │ 0.02    │ 47.0    │ treated │
+│ 3   │ 0.06    │ 97.0    │ treated │
+⋮
+│ 9   │ 0.56    │ 191.0   │ treated │
+│ 10  │ 0.56    │ 201.0   │ treated │
+│ 11  │ 1.1     │ 207.0   │ treated │
+│ 12  │ 1.1     │ 200.0   │ treated │
+
+julia> pm1 = fit(NLregModel, PurTrt, :rate, (p,d) -> p.Vm * d.conc/(p.K + d.conc),
+                  (Vm = 200., K = 0.05))
+Nonlinear regression model fit by maximum likelihood
+
+Data schema (response variable is rate)
+Tables.Schema:
+ :conc   Float64
+ :rate   Float64
+ :state  String
+
+ Sum of squared residuals at convergence: 1195.4488145417758
+ Achieved convergence criterion:          8.798637504793927e-6
+
+ Number of observations:                  12
+
+ Parameter estimates
+───────────────────────────────────────
+      Estimate   Std.Error  t-statistic
+───────────────────────────────────────
+Vm  212.684     6.94715        30.6145
+K     0.064121  0.00828092      7.74322
+───────────────────────────────────────
 ```
-## Creating a PLregMod type
-
-A `PLregMod` type contains the transposed gradient, usually called
-`tgrad` with the conditionally linear parameters first, the
-three-dimensional Jacobian array, usually called `MMD`, with each face
-corresponding to the partial derivative of the conditionally linear
-rows of `tgrad` with respect to the nonlinear parameters, and the
-values of any covariates needed to evaluate the model.  The
-model-matrix function, `mmf`, is a function of two read-only
-arguments; `nlp`, the nonlinear parameters in the model function as a
-vector and x the covariates for a single observation, also as a
-vector, and two arguments, `tg` and `MMD` that are updated in the
-function.  The `tg` vector is updated with the model matrix for the
-conditionally linear parameters and the `MMD` slice, considered as a
-matrix, is updated with the gradient.
-
-For the Michaelis-Menten model the model-matrix function is
-```julia
-function MicMenmmf(nlp,x,tg,MMD)
-    x1 = x[1]
-    denom = nlp[1] + x1
-    MMD[1,1] = -(tg[1] =  x1/denom)/denom
-end
-```
-The arguments are untyped, to allow for submatrices or views of
-matrices and arrays, but they should be treated as three vectors and a
-matrix, for the purposes of indexing.
-
-
-

doc/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

doc/make.jl

@@ -0,0 +1,17 @@
+using Documenter, NLreg
+
+makedocs(;
+    modules=[NLreg],
+    format=Documenter.HTML(),
+    pages=[
+        "Home" => "index.md",
+    ],
+    repo="https://github.com/dmbates/NLreg.jl/blob/{commit}{path}#L{line}",
+    sitename="NLreg.jl",
+    authors="Douglas Bates <[email protected]>",
+    assets=String[],
+)
+
+deploydocs(;
+    repo="github.com/dmbates/NLreg.jl",
+)

doc/src/index.md

@@ -0,0 +1,8 @@
+# NLreg.jl
+
+```@index
+```
+
+```@autodocs
+Modules = [NLreg]
+```