Julia GLM - используя devresid для построения графиков - PullRequest
Новая
       11

Julia GLM - используя devresid для построения графиков

1 голос
/ 21 апреля 2020

Я хотел бы сделать некоторый остаточный анализ для GLM.

Моя модель имеет вид 

using GLM
model = glm(@formula(y ~ x), data, Binomial(), LogitLink())

Мой учебник предлагает, чтобы остаточный анализ в GLM проводился с использованием остаточных отклонений. Я был рад видеть, что GLM Джулии имеет функцию devresid(), и что она предлагает, как использовать ее для построения графиков (sign(y - μ) * sqrt(devresid(D, y, μ))). Тем не менее, я в полной растерянности относительно того, какими должны быть входные аргументы. Глядя на строку do c: 

?devresid
devresid(D, y, μ::Real)

Return the squared deviance residual of μ from y for distribution D

The deviance of a GLM can be evaluated as the sum of the squared deviance residuals. This is the principal use for these values. The actual deviance residual, say for plotting, is the signed square root of this value

sign(y - μ) * sqrt(devresid(D, y, μ))

Examples

julia> devresid(Normal(), 0, 0.25) ≈ abs2(0.25)
true

julia> devresid(Bernoulli(), 1, 0.75) ≈ -2*log(0.75)
true

julia> devresid(Bernoulli(), 0, 0.25) ≈ -2*log1p(-0.25)
true
  • D: Я предполагаю, что в моем случае это Binomial()
  • y: Я предполагаю это переменная индикатора для отдельного случая, т.е. 1 или 0
  • μ: что это такое?

Как я могу использовать эту функцию для создания таких вещей, как графики Остаточное отклонение по нормальной шкале вероятности и от подгоночных значений?

Вот данные, которые я использую в форме CSV 

x,y
400,0
220,1
490,0
210,1
500,0
270,0
200,1
470,0
480,0
310,1
240,1
490,0
420,0
330,1
280,1
210,1
300,1
470,1
230,0
430,0
460,0
220,1
250,1
200,1
390,0

1 Ответ

1 голос
/ 21 апреля 2020

Я понимаю, что это то, что вы хотите: 

julia> data = DataFrame(X=[1,1,1,2,2], Y=[1,1,0,0,1])
5×2 DataFrame
│ Row │ X     │ Y     │
│     │ Int64 │ Int64 │
├─────┼───────┼───────┤
│ 1   │ 1     │ 1     │
│ 2   │ 1     │ 1     │
│ 3   │ 1     │ 0     │
│ 4   │ 2     │ 0     │
│ 5   │ 2     │ 1     │

julia> model = glm(@formula(Y ~ X), data, Binomial(), LogitLink())
StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},Binomial{Float64},LogitLink},GLM.DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}

Y ~ 1 + X

Coefficients:
─────────────────────────────────────────────────────────────────────────────
              Estimate  Std. Error    z value  Pr(>|z|)  Lower 95%  Upper 95%
─────────────────────────────────────────────────────────────────────────────
(Intercept)   1.38629      2.82752   0.490286    0.6239   -4.15554    6.92813
X            -0.693146     1.87049  -0.37057     0.7110   -4.35923    2.97294
─────────────────────────────────────────────────────────────────────────────

julia> p = predict(model)
5-element Array{Float64,1}:
 0.6666664218508201
 0.6666664218508201
 0.6666664218508201
 0.5
 0.5

julia> y = data.Y
5-element Array{Int64,1}:
 1
 1
 0
 0
 1

julia> @. sign(y - p) * sqrt(devresid(Bernoulli(), y, p))
5-element Array{Float64,1}:
  0.9005170462928523
  0.9005170462928523
 -1.4823033118905455
 -1.1774100225154747
  1.1774100225154747

(это то, что вы получите от вызова residuals(model, type="deviance") в R)

Обратите внимание, что в последней строке я использую @. для векторизации всей строки. В качестве альтернативы вы могли бы написать это как: 

julia> sign.(y .- p) .* sqrt.(devresid.(Bernoulli(), y, p))
5-element Array{Float64,1}:
  0.9005170462928523
  0.9005170462928523
 -1.4823033118905455
 -1.1774100225154747
  1.1774100225154747
...