using MultinomialRegression
using DataFrames
using RDatasets
# Data
iris = dataset("datasets", "iris")
# Unregularized fit
model = fit(@formula(Species ~ 1 + SepalWidth), iris) # levels(iris.Species)[1] is the reference category
opts = Dict(:iterations => 250, :f_abstol => 1e-9) # Same terminology as Optim.Options
model = fit(@formula(Species ~ 1 + SepalWidth), iris; opts=opts)
# Predict
xnew = [1.0, iris.SepalWidth[1]]
pred = predict(model, xnew)
pred2 = zeros(3)
predict!(pred2, model, xnew)
# Model-level diagnostics
isregularized(model)
nobs(model)
loglikelihood(model)
aic(model)
aicc(model)
bic(model)
# Coefficient-level diagnostics
nparams(model)
coef(model)
coeftable(model)
stderror(model)
vcov(model)
coefcor(model)
# Get values using row names and column names
B = coef(model)
B["(Intercept)", "virginica"]
B["(Intercept)", :]
B[:, "virginica"]
B[1, 2] # Integer indices also work
B2 = [18.85825175905592 12.997166560384215; -6.118905216314172 -4.079050714520578]
isapprox(B, B2; atol=1e-10) # Reproducible result
#=
Regularized fit
Note: Standard errors calculated from the hessian evaluated at the parameters.
Since the parameters of regularized models are not MLEs, the standard errors are approximate.
=#
model_L1 = fit(@formula(Species ~ 1 + SepalWidth), iris; reg=L1(0.5))
model_L2 = fit(@formula(Species ~ 1 + SepalWidth), iris; reg=L2(0.5))
model_ElNet = fit(@formula(Species ~ 1 + SepalWidth), iris; reg=ElasticNet(0.5, 0.5))
# Weighted fit
w = collect(0.25:0.01:1.75)
splice!(w, findfirst(==(1.0), w))
weighted_fit = fit(@formula(Species ~ 1 + SepalWidth), iris; wts=w)
# Under the hood, weights are scaled to sum to the number of observations
w2 = 2*w
weighted_fit2 = fit(@formula(Species ~ 1 + SepalWidth), iris; wts=w2)
isapprox(coef(weighted_fit2), coef(weighted_fit); atol=1e-10)