CODEWITHSUNDEEP
rjuliaglm
Guy Sutton
Guy Sutton

Reputation: 45

Julia equivalent of R's 'hatvalue' function

Is there a Julia equivalent of R's hatvalues function? hatvalues calculates leverage values from a fitted linear model. Ultimately, I would like to calculate Standardised Pearson's residuals from a fitted GLM, but it seems like one will have to do this manually in Julia, which requires the calculation of hatvalues. For example, this is exactly the output I want and get from R.

library(tidyverse)
library(glmmTMB)

salamanders <- glmmTMB::Salamanders

mod1 <- glm(count ~ 1 + mined,
            data = salamanders,
            family = poisson(link = "log"))
head(hatvalues(mod1))

>  head(hatvalues(mod1))
1           2           3           4           5           6 
0.003246753 0.003246753 0.003246753 0.002976190 0.002976190 0.002976190

I have tried to use StatsBase.leverage() in Julia, however, this doesn't seem to work for the object type obtained when fitting a GLM (`GLM.jl).

# Load required packages 
using RCall
using Distributions
using Random
using DataFrames
using StatsBase
using GLM
using Statistics
using Compose

# Load in the dataset 
salamander = rcopy(R"glmmTMB::Salamanders")

# Run Poisson GLM
m1 = fit(GeneralizedLinearModel,
            @formula(count ~ 1 + mined),
            salamander,
            Poisson(),
            GLM.LogLink())

# Try to calculate leverage
StatsBase.leverage(m1)

julia> StatsBase.leverage(m1)
ERROR: leverage is not defined for StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},Poisson{Float64},LogLink},GLM.DensePredChol{Float64,Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}.
Stacktrace:
 [1] error(::String) at .\error.jl:33
 [2] leverage(::StatsModels.TableRegressionModel{GeneralizedLinearModel{GLM.GlmResp{Array{Float64,1},Poisson{Float64},LogLink},GLM.DensePredChol{Float64,Cholesky{Float64,Array{Float64,2}}}},Array{Float64,2}}) at C:\Users\s1900332\.julia\packages\StatsBase\XT7PT\src\statmodels.jl:354
 [3] top-level scope at REPL[295]:1

Upvotes: 1

Views: 155

Answers (1)

G. Grothendieck
G. Grothendieck

Reputation: 269694

Run the corresponding linear model and apply hatvalues to that.

library(glmmTMB)
salamanders <- glmmTMB::Salamanders

fo <- count ~ 1 + mined

mod1 <- glm(fo, data = salamanders, family = poisson(link = "log"))
hv1 <- hatvalues(mod1)

mod2 <- lm(fo, data = salamanders)
hv2 <- hatvalues(mod2)

all.equal(hv1, hv2)
## [1] TRUE

Upvotes: 2

Related Questions