Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

backport of QR to 1.x #537

Open
wants to merge 9 commits into
base: v1
Choose a base branch
from
Open

backport of QR to 1.x #537

Changes from 1 commit
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
cd7b2b1
`QR` decomposition method in GLM.
mousum-github May 4, 2023
2d91725
Link(::TableRegressionModel)
palday Jun 7, 2023
97cab57
oh contravariance
palday Jun 7, 2023
ae88551
oh TableRegressionModel
palday Jun 7, 2023
173d3b7
predict got more complicated
palday Jun 7, 2023
5269e7a
test updates
palday Jun 7, 2023
23321a8
version bump
palday Jun 7, 2023
92ec3fa
doctest fixes
palday Jun 7, 2023
cd187fb
Apply suggestions from code review
palday Jun 7, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 19 additions & 22 deletions src/lm.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -256,8 +256,23 @@ Valid values of `interval` are `:confidence` delimiting the uncertainty of the
predicted relationship, and `:prediction` delimiting estimated bounds for new data points.
"""
function predict(mm::LinearModel, newx::AbstractMatrix;
interval::Union{Symbol,Nothing}=nothing, level::Real = 0.95)
retmean = newx * coef(mm)
interval::Union{Symbol,Nothing}=nothing, level::Real=0.95)
retmean = similar(view(newx, :, 1))
if interval === nothing
res = retmean
predict!(res, mm, newx)
else
res = (prediction=retmean, lower=similar(retmean), upper=similar(retmean))
predict!(res, mm, newx, interval=interval, level=level)
end
return res
end
function StatsModels.predict!(res::Union{AbstractVector,
NamedTuple{(:prediction, :lower, :upper),
<:NTuple{3, AbstractVector}}},
mm::LinearModel, newx::AbstractMatrix;
interval::Union{Symbol, Nothing}=nothing,
level::Real=0.95)
if interval === :confint
Base.depwarn("interval=:confint is deprecated in favor of interval=:confidence", :predict)
interval = :confidence
Expand Down Expand Up @@ -286,7 +301,6 @@ function predict(mm::LinearModel, newx::AbstractMatrix;
length(prediction) == length(lower) == length(upper) == size(newx, 1) ||
throw(DimensionMismatch("length of vectors in `res` must equal the number of rows in `newx`"))
length(mm.rr.wts) == 0 || error("prediction with confidence intervals not yet implemented for weighted regression")

dev = deviance(mm)
dofr = dof_residual(mm)
ret = diag(newx*vcov(mm)*newx')
Expand All @@ -300,27 +314,10 @@ function predict(mm::LinearModel, newx::AbstractMatrix;
lower .= prediction .+ ret
upper .= prediction -+ ret
end
length(mm.rr.wts) == 0 || error("prediction with confidence intervals not yet implemented for weighted regression")
chol = cholesky!(mm.pp)
# get the R matrix from the QR factorization
if chol isa CholeskyPivoted
ip = invperm(chol.p)
R = chol.U[ip, ip]
else
R = chol.U
end
residvar = ones(size(newx,2)) * deviance(mm)/dof_residual(mm)
if interval == :confidence
retvariance = (newx/R).^2 * residvar
elseif interval == :prediction
retvariance = (newx/R).^2 * residvar .+ deviance(mm)/dof_residual(mm)
else
error("only :confidence and :prediction intervals are defined")
end
retinterval = quantile(TDist(dof_residual(mm)), (1. - level)/2) * sqrt.(retvariance)
(prediction = retmean, lower = retmean .+ retinterval, upper = retmean .- retinterval)
return res
end


Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Intended?

function confint(obj::LinearModel; level::Real=0.95)
hcat(coef(obj),coef(obj)) + stderror(obj) *
quantile(TDist(dof_residual(obj)), (1. - level)/2.) * [1. -1.]
Expand Down