# Auxiliary GLM Statistics

**URL:** <https://discourse.julialang.org/t/auxiliary-glm-statistics/14722>\
**Category:** Statistics\
**Tags:** glm\
**Created:** [September 9, 2018, 2:37am UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722 "2018-09-09T02:37:28Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![iwelch](https://avatars.discourse-cdn.com/v4/letter/i/8c91f0/32.png) [@iwelch](https://discourse.julialang.org/u/iwelch)\
**Post date:** [September 9, 2018, 2:37am UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/1 "2018-09-09T02:37:28Z")

</div>

Julia 1.0. GLM 1.0.0.

### My Basic Example

```julia-auto
julia> using GLM, StatsBase

julia> y=[1:10;] ; X= hcat( fill(1,10), y.^2, y.^3 ); lm1= GLM.lm(X, y);

```

### Basic Questions

1. is the sigma (i.e., `sqrt(deviance(lm)/dof_residual(lm)`) already programmed into julia?

2. is the f-statistic of the model (i.e., whether the model beats a constant) already available somewhere ?

My motivation is to try to write a function that prints regression output in a different format that I prefer. At first, I thought I would write a `show` function, but then I figured that the spirit of julia is for the user not to replace the basics, but to write different functions…or maybe add a particular kind of keyword that dispatches to my function?

### Less Important Question on ftest

1. how does ftest work? The example from the docs does not work:

```julia-auto
julia> using DataFrames, GLM, StatsBase

julia> dat = DataFrame(Treatment=[1, 1, 1, 2, 2, 2, 1, 1, 1, 2, 2, 2.],
                                Result=[1.1, 1.2, 1, 2.2, 1.9, 2, .9, 1, 1, 2.2, 2, 2],
                                Other=[1, 1, 2, 1, 2, 1, 3, 1, 1, 2, 2, 1]);

julia> using StatsModels

julia> mod = lm(@formula(Result ~ 1 + Treatment), dat);
ERROR: cannot assign variable Base.mod from module Main
Stacktrace:
 [1] top-level scope at none:0

julia> mod = lm(@formula(Result ~ 1 + Treatment), dat);
ERROR: cannot assign variable Base.mod from module Main
Stacktrace:
 [1] top-level scope at none:0

```

so, I tried to make up my own simpler example and without DataFrame:

```julia-auto
julia> using GLM, StatsBase

julia> y=[1:10;] ; X= hcat( fill(1,10), y.^2, y.^3 );

julia> y=[1:10;] ; XX= hcat( fill(1,10), y.^2, y.^3, y.^4 );

julia> lm1= GLM.lm(X, y); lm2= GLM.lm(XX, y);

julia> ftest(lm2, lm1)
Error showing value of type GLM.FTestResult{2}:
ERROR: MethodError: no method matching Array{String,2}(::Int64, ::Int64)
Closest candidates are:
  Array{String,2}(::UndefInitializer, ::Int64, ::Int64) where T at boot.jl:396
  Array{String,2}(::UndefInitializer, ::Int64...) where {T, N} at boot.jl:400
  Array{String,2}(::UndefInitializer, ::Integer, ::Integer) where T at sysimg.jl:143
  ...
Stacktrace:
 [1] show(::IOContext{REPL.Terminals.TTYTerminal}, ::GLM.FTestResult{2}) at /Users/ivo/.julia/packages/GLM/8TwVM/src/ftest.jl:110

```

### Prediction?

in R, one gets a predicted value like this. How do I do it in this simple example?

```julia-auto
julia> using DataFrames, GLM

julia> y=[1:10;] ; X= hcat( fill(1,10), y.^2, y.^3 ); lm1= GLM.lm(X, y);

julia> predict( lm1, [1.0, 2.0, 3.0])
ERROR: MethodError: no method matching predict(::LinearModel{LmResp{Array{Float64,1}},DensePredChol{Float64,LinearAlgebra.Cholesky{Float64,Array{Float64,2}}}}, ::Array{Float64,1})

```

### pointers appreciated.

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dominiquemakowski/32/51410_2.png) [@DominiqueMakowski](https://discourse.julialang.org/u/DominiqueMakowski)\
**Post date:** [September 10, 2018, 6:14am UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/2 "2018-09-10T06:14:04Z")

</div>

To me this is illustrative of the paradox currently existing in Julia: you can do cutting-edge super crazy complex stuff but there is currently no straightforward function to perform a Pearson’s correlation test 😅 I believe the implementation of these “basic” level statistics (including GLMs fitting) will be necessary for Julia’s adoption to really explode (there are fields with millions of researchers almost exclusively (but intensively) using LMs or ANOVAs).

With that being said, after trying to implement some [helper functions for GLM’s](https://github.com/neuropsychology/Psycho.jl), I similarly found it fairly uneasy to deal with the structures and extract things from the models ☹

I believe this is gonna change with time as they are refactoring things (such as the [formula system](https://github.com/JuliaStats/StatsModels.jl/pull/71)) but we might have to be a bit more patient.

Sorry for not answering your question!

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [September 10, 2018, 7:23am UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/3 "2018-09-10T07:23:44Z")

</div>

> [@DominiqueMakowski](#):
>
> there is currently no straightforward function to perform a Pearson’s correlation test 😅

This is a particularly unfortunate example, since there are [at least 5 ways to perform frequentist inference on Pearson’s \rho](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Inference), depending on the assumptions you want to make.

If you believe that there is a “straightforward” way to do something just because a particular language/framework that you have been using has that implemented as the default, you may expect to be disappointed with the Julia package ecosystem, as it is common to implement things in a modular way, giving the user a lot of choices. If, however, you understand what you are doing, you may find the latter very convenient.

---

<div class="post-metadata">

**Author:** ![mcreel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcreel/32/30088_2.png) [@mcreel](https://discourse.julialang.org/u/mcreel)\
**Post date:** [September 10, 2018, 7:33am UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/4 "2018-09-10T07:33:23Z")

</div>

It’s `predict(lm1, [1 2 3])`. The input you supply must be the same type as that used to fit the model. In this case, X is an array of integers, so you need to use integers for the values you supply for the out of sample fit.

More info, from ?predict  
Form the predicted response of model obj. An object with new covariate values newX can be supplied, which should have the **same type and structure** as that used to fit obj; e.g. for a GLM it would generally be a DataFrame with the same variable names as the original predictors.

Regarding methods and functions, I suspect that the statistics/data science environment will become more complete and perhaps somewhat more homogeneous now that Julia 1.0 is here. I will be looking to use DataFrames, StatsModels, etc. more in my own work.

---

<div class="post-metadata">

**Author:** ![DominiqueMakowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dominiquemakowski/32/51410_2.png) [@DominiqueMakowski](https://discourse.julialang.org/u/DominiqueMakowski)\
**Post date:** [September 10, 2018, 7:53am UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/5 "2018-09-10T07:53:39Z")

</div>

> [@Tamas\_Papp](#):
>
> This is a particularly unfortunate example, since there are [at least 5 ways to perform frequentist inference on Pearson’s ρ\rho](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Inference), depending on the assumptions you want to make.

You’re absolutely right! By “straightforward” I meant R-like with R-like defaults (the language which OP was referring to), which clearly might not be the only / best / recommended way of doing things. However, a “default” implementation is usually how people (at least in my field) start with statistics and programming, by running `cor.test(x, y)` and then (hopefully) wondering: “wait a sec, what does this parameter / function _actually_ do?”

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [September 10, 2018, 2:51pm UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/6 "2018-09-10T14:51:54Z")

</div>

> [@iwelch](#):
>
> how does ftest work? The example from the docs does not work:

Good catch. See [Fix ftest printing by nalimilan · Pull Request #252 · JuliaStats/GLM.jl · GitHub](https://github.com/JuliaStats/GLM.jl/pull/252).

---

<div class="post-metadata">

**Author:** ![iwelch](https://avatars.discourse-cdn.com/v4/letter/i/8c91f0/32.png) [@iwelch](https://discourse.julialang.org/u/iwelch)\
**Post date:** [September 10, 2018, 6:53pm UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/7 "2018-09-10T18:53:50Z")

</div>

thx, nalimilan. now I get the same error in both cases.

```julia
julia> ft = ftest(model.model, nullmodel.model)
Error showing value of type GLM.FTestResult{2}:
ERROR: MethodError: no method matching Array{String,2}(::Int64, ::Int64)
Closest candidates are:
  Array{String,2}(::UndefInitializer, ::Int64, ::Int64) where T at boot.jl:396
  Array{String,2}(::UndefInitializer, ::Int64...) where {T, N} at boot.jl:400
  Array{String,2}(::UndefInitializer, ::Integer, ::Integer) where T at sysimg.jl:143
  ...

```

regards, /iaw

PS: (thanks mcreel for the predict info. and is the regression stderr of the resids and overall F-test really not provided? I am not asking inference. I am asking about basic stats.)

---

<div class="post-metadata">

**Author:** ![compleat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/compleat/32/8958_2.png) [@compleat](https://discourse.julialang.org/u/compleat)\
**Post date:** [January 27, 2021, 6:11pm UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/8 "2021-01-27T18:11:38Z")

</div>

Hi All, I am a latecomer to this discussion (maybe it is all obsolete now), but I found a handy workaround for ftest, given that the ‘r2’ function is there:

```julia
ftest(lmodel)=(1/(1-r2(lmodel))-1)*(length(predict(lmodel))- 
    length(coef(lmodel)))/(length(coef(lmodel))-1)

```

It would be cool if this could be implemented, as I use this interface for my introductory stats classes.

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [January 27, 2021, 6:24pm UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/9 "2021-01-27T18:24:45Z")

</div>

Feel free to make a pull request to implement a single-argument `ftest` method. But it should probably return the same kind of table as the multiple argument method (not a single number).

---

<div class="post-metadata">

**Author:** ![compleat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/compleat/32/8958_2.png) [@compleat](https://discourse.julialang.org/u/compleat)\
**Post date:** [January 27, 2021, 6:35pm UTC](https://discourse.julialang.org/t/auxiliary-glm-statistics/14722/10 "2021-01-27T18:35:30Z")

</div>

Thanks. I tried to, but unfortunately I don’t know how to do this (I got bogged down with ‘branch comparisons’). I take your point about the name … maybe call it ‘fstat’ because it is a standard diagnostic, or else make it an allowable method for ftest)