# Domain error when running user-defined function, but not when code is run outside function
**URL:** https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771
**Category:** New to Julia
**Tags:** function, glm
**Created:** [March 23, 2021, 7:39am UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771 "2021-03-23T07:39:44Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@guysutton](https://discourse.julialang.org/u/guysutton)
#### Post date: [March 23, 2021, 7:39am UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771/1 "2021-03-23T07:39:45Z")
I am trying to write a function that will calculate deviance residuals for a general linear model (specified using `GLM.jl`). The function extracts predictions from the fitted model (`pred`), a vector containing the response variable (`response_y`) and the error distribution specified for the fitted model (`distType`). I then use a conditional if-elseif-else statement to calculate the deviance residuals depending on the output of `distType`.
```julia
# Write function
using Test
using Distributions
using Random
using DataFrames
using StatsBase
using GLM
using Gadfly
using Statistics
using Compose
# Define function to calculate deviance residuals
function calcDevResids(;model, df)
###################################################
# Section #1: Extract response vector and model residuals
###################################################
# Extract fitted/predicted values from model object
pred = GLM.predict(model)
# Extract vector containing response variable
response_y = df.y
###################################################
# Section #2: Calculate deviance residuals
# - requires conditional statements to
# loop through different data distributions
###################################################
# Extract statistical distribution specified in model
distType = typeof(model).parameters[1].parameters[1].parameters[2]
# Use if-elseif-else statements to calculate deviance residuals depending
# on which errors were specified during model specification
# If model was specified with Poisson errors
if distType == "Poisson"
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Poisson(), response_y, pred))
# If model was specified with Gamma errors
elseif distType == "Gamma"
return sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Gamma(), response_y, pred))
# If model was specified with Bernoulli errors
elseif distType == "Bernoulli"
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Bernoulli(), response_y, pred))
# If model was specified with Binomial errors
else distType == "Binomial"
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Binomial(), response_y, pred))
end
end
```
And here is a basic model to test the function:
```julia
# Set reproducible seed
Random.seed!(123)
# Simulate n = 100 data points from a Poisson distribution (lambda = 7)
y = rand(Poisson(7), 100)
# Assign grouping variables (5 groups of 20)
x = repeat([1, 2, 3, 4, 5],
inner = 20,
outer = 1)
df = hcat(y, x)
df = DataFrame(df)
colnames = ["y","x"]
rename!(df, Symbol.(colnames))
# Run model
modPoisson = fit(GeneralizedLinearModel,
@formula(y ~ 1 + x),
df,
Poisson(),
GLM.LogLink())
```
However, I keep getting the follow error when I run the function:
`ERROR: DomainError with -10.0: log1p will only return a complex result if called with a complex argument. Try log1p(Complex(x)).`
But, when I run the code outside of the function call (see below), I get the expected output and no error.
```julia
pred = GLM.predict(modPoisson)
response_y = df.y
dv = sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Poisson(), response_y, pred))
```
`julia> dv = sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Poisson(), response_y, pred)) 100-element Array{Float64,1}: 0.9061633419878498 1.5493622873332258 0.9061633419878498 0.217165948855438 -0.9383362226732821 -0.9383362226732821 ⋮ 0.07437940612764843 0.07437940612764843 0.8013566375722013 0.07437940612764843 1.4748443746174682 -0.7265187364238832`
I have tried adding `complex()` and `Complex()` to the terms in the code to create `dv`, with no success.
Any pointers or advice would be much appreciated. I am a Julia-newbie (\< 3 months experience), so I am sure I am missing something quite simple.
---
### Author:  [@sijo](https://discourse.julialang.org/u/sijo)
#### Post date: [March 23, 2021, 10:39am UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771/2 "2021-03-23T10:39:19Z")
Try to add `@show distType` after its definition in `calcDevResids`. It prints the following:
```julia
distType = Distributions.Poisson{Float64}
```
so the test `if distType == "Poisson"` will fail: `distType` is not the string “Poisson” but the type `Distributions.Poisson{Float64}`. Instead, you can write `if distType <: Poisson`.
A couple of remarks:
- The lone `return` looks weird.
- In `else distType == "Binomial"` the comparison is useless since there’s no `if`.
- Instead of
```julia
df = hcat(y, x)
df = DataFrame(df)
colnames = ["y","x"]
rename!(df, Symbol.(colnames))
```
you can write
```julia
DataFrame(; y, x)
```
---
### Author:  [@hendri54](https://discourse.julialang.org/u/hendri54)
#### Post date: [March 23, 2021, 12:31pm UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771/3 "2021-03-23T12:31:54Z")
Instead of
```julia
distType = typeof(model).parameters[1].parameters[1].parameters[2]
```
I would prefer the more robust
```julia
resp_type(m :: GeneralizedLinearModel{T1, T2}) where {T1, T2} = T1;
distr_type(::Type{GLM.GlmResp{T1, T2, T3}}) where {T1,T2,T3} = T2;
distType2 = distr_type(resp_type(model.model));
```
---
### Author:  [@hendri54](https://discourse.julialang.org/u/hendri54)
#### Post date: [March 23, 2021, 12:41pm UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771/4 "2021-03-23T12:41:42Z")
Then it would be nice to define
```julia
deviance_residual(D, response_y, pred) =
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(D, response_y, pred))
```
If one could then construct the distribution from its type somehow (from the return arg of `distr_type` above), the entire `if ... else` loop could be replaced with
```julia
distType = distr_type(resp_type(model.model));
# Construct the distribution object - how?
D = ??
return deviance_residual(D, response_y, pred)
```
However, I am not sure how to call the constructor `Poisson()` from the data type `Poisson{Float64}` that `distr_type` returns.
---
### Author:  [@guysutton](https://discourse.julialang.org/u/guysutton)
#### Post date: [March 23, 2021, 3:15pm UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771/5 "2021-03-23T15:15:07Z")
@hendri54 Thanks for your suggestions - I am going to need some time to work through your suggestions and wrap my head around the code.
@sijo. This seemed to solve the problem - Thanks very much. However, it created another problem later on in the workflow. When I pass the example model to the updated function, I get a new error message. Can you maybe advise as to what I am missing? Full updated code below:
```julia
##########################
# WRITE UPDATED FUNCTION
##########################
using Test
using Distributions
using Random
using DataFrames
using StatsBase
using GLM
using Statistics
using Compose
# Define function to calculate deviance residuals
function calcDevResids(;model, df)
###################################################
# Section #1: Extract response vector and model residuals
###################################################
# Extract fitted/predicted values from model object
pred = GLM.predict(model)
# Extract vector containing response variable
response_y = df.y
###################################################
# Section #2: Calculate deviance residuals
# - requires conditional statements to
# loop through different data distributions
###################################################
# Extract statistical distribution specified in model
distType = typeof(model).parameters[1].parameters[1].parameters[2]
# Use if-elseif-else statements to calculate deviance residuals depending
# on which errors were specified during model specification
# If model was specified with Poisson errors
if distType <: Poisson
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Poisson(), response_y, pred))
# If model was specified with negative binomial errors
elseif distType <: NegativeBinomial
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.NegativeBinomial(), response_y, pred))
# If model was specified with Gamma errors
elseif distType <: Gamma
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Gamma(), response_y, pred))
# If model was specified with Bernoulli errors
elseif distType <: Bernoulli
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Bernoulli(), response_y, pred))
# If model was specified with Binomial errors
elseif distType <: Binomial
sign.(response_y .- pred) .* sqrt.(GLM.devresid.(Distributions.Binomial(), response_y, pred))
end
end
###########################
# MAKE SOME EXAMPLE DATA
###########################
# Set reproducible seed
Random.seed!(123)
# Simulate n = 100 data points from a Poisson distribution (lambda = 7)
y = rand(Poisson(7), 100)
# Assign grouping variables (5 groups of 20)
x = repeat([1, 2, 3, 4, 5],
inner = 20,
outer = 1)
df = DataFrame(; y, x)
# Run model
modPoisson = fit(GeneralizedLinearModel,
@formula(y ~ 1 + x),
df,
Poisson(),
GLM.LogLink())
```
Applying the updated function returns the new error. I assume that I am not declaring the condition arguments in the if-elseif statements properly (i.e. the distType \<: Poisson/Bernoulli/Binomial/Gamma parts will all return the same error)? I didn’t think I could/should declare the conditions passed to `distType <: x`?
```julia
##################################################
# REPRODUCE NEW ERROR WHEN APPLYING FUNCTION
##################################################
calcDevResids(model = modPoisson, df = df)
```
This returns the following error.
`ERROR: UndefVarError: Poisson not defined`
---
### Author:  [@sijo](https://discourse.julialang.org/u/sijo)
#### Post date: [March 23, 2021, 3:47pm UTC](https://discourse.julialang.org/t/domain-error-when-running-user-defined-function-but-not-when-code-is-run-outside-function/57771/6 "2021-03-23T15:47:05Z")
For the conditions, they look fine to me (assuming `Gamma`, etc. are the correct types for the respective cases. I would add a print statement in each if/else to check everything works as you expect.
The code actually runs without error for me (on Julia 1.6-rc3). The `Poisson` symbol should be defined since you use Distributions (even without that, it should be defined (reexported) by GLM). So maybe there’s something wrong in your session… Did you try executing this code in a fresh Julia session?
In any case it would be useful to have the full error message.