# Confidence Interval for certain Contrast and model residual

**URL:** <https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339>\
**Category:** Statistics\
**Created:** [February 1, 2019, 1:34am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339 "2019-02-01T01:34:47Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 1, 2019, 1:34am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/1 "2019-02-01T01:34:47Z")

</div>

I have a model like this:

`ols = lm(@formula(LnVar ~ Seq+Per+Trt+Subj), df)`

I try to get confidence interval for Trt contrast, LSM and model residual. How to obtain it?

In R you just do `confint(lmobj, c("Trt"), level=0.9)` and `summary(lmobj )$sigma^2`

And I can not find any working example for this. And model residual i can not find at all. How to make this kind of calculations. Can i do it “from the box”?

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [February 1, 2019, 6:42am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/2 "2019-02-01T06:42:41Z")

</div>

Please post a minimum working example of your code that allows others to reproduce what you’re doing.

In particular, it’s not clear which package you are using. In [GLM.jl](http://juliastats.github.io/GLM.jl/stable/manual/) you would use `stderror` to get the standard errors of your coefficients.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [February 1, 2019, 7:59am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/3 "2019-02-01T07:59:20Z")

</div>

You can get the confidence interval of the coefficients with `confint(ols)`, and the residuals with `residuals(ols)`. If that doesn’t work or you’re having problems post back here.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 1, 2019, 4:12pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/4 "2019-02-01T16:12:59Z")

</div>

A have a dataset 12248\_2014\_9661\_MOESM1\_ESM.txt from here:

> **[Reference datasets for 2-treatment, 2-sequence, 2-period bioequivalence...](https://pubmed.ncbi.nlm.nih.gov/25212768/)**
>
> It is difficult to validate statistical software used to assess bioequivalence since very few datasets with known results are in the public domain, and the few that are published are of moderate size and balanced. The purpose of this paper is...

and i try to validate Julia code for BE studies:

df = readtable(“12248\_2014\_9661\_MOESM1\_ESM.txt”, header = true, separator = ‘\t’)  
df.LnVar = log.(df.Var)  
df.Subj = string.(df.Subj)  
ols = lm(@formula(LnVar ~ Seq+Per+Trt+Subj), df)  
cint = confint(ols, 0.9)  
print(“Lower:”)  
println(exp(cint[4,1]))  
print(“Upper:”)  
println(exp(cint[4,2]))

I suppose to get 90% CI for Trt LSM.  
For this dataset i have:

Lower:0.9060577396300195  
Upper:0.997880935149158

If i run R code:

lmobj ← lm(log(Var)~Trt+Per+Seq+Subj, data=pkdata)  
lnCI ← confint(lmobj, c(“TrtT”), level=0.9)  
exp(lnCI)

I have another results:  
TrtT 0.9076208 0.9961624

In R i get residual variance by call:

summary(lmobj )$sigma^2

How to get it from GLM in Julia i do no know.

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [February 2, 2019, 9:42am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/5 "2019-02-02T09:42:40Z")

</div>

Hi @PharmCat it seems like it’s worth opening an issue on GLM to discuss this, which will catch the attention on the developers. In particular, R and Julia seem to deal differently with subj 8, which is abandoned in R but a value is estimated in Julia. This changes the residual degrees of freedom from 16 to 15. I don’t think that’s the full story though, as simply skipping subj 8 from the analysis does not realign the results.  
I don’t think there is a function in Julia to optain the residual variance (there isn’t one in R either, as your example shows).  
A quick note on presentation - it’s best to surround code by a block of triple backticks, and to include the `using CSV, DataFrames, GLM, StatsModels` part of the code. Also, `CSV.read("12248_2014_9661_MOESM1_ESM.txt", delim = '\t')` is the preferred syntax for reading DataFrames today.

Finally, here’s a direct link to the data set [https://static-content.springer.com/esm/art%3A10.1208%2Fs12248-014-9661-0/MediaObjects/12248\_2014\_9661\_MOESM1\_ESM.txt](https://static-content.springer.com/esm/art%3A10.1208%2Fs12248-014-9661-0/MediaObjects/12248_2014_9661_MOESM1_ESM.txt)

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 2, 2019, 10:27am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/6 "2019-02-02T10:27:42Z")

</div>

Hello! Thank you for explanation! In this example design matrix is singular and R use QR decomposition with pivoting to get coefficients. And in this example Subj in nested in sequence, but R calculate df without direct settings. I So, should i make an issue on github?

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [February 2, 2019, 10:39am UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/7 "2019-02-02T10:39:39Z")

</div>

Yes, I think so. It would be useful in the issue to mention that this dataset is from a collection specifically invented to present edge cases for assessing the robustness of statistical software implementations.

---

<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:** [February 2, 2019, 3:37pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/8 "2019-02-02T15:37:24Z")

</div>

You can get the residual variance a.k.a. dispersion parameter using the unexported function `dispersion`: `GLM.dispersion(lmobj.model, true)`. We should probably export it (and it already has a docstring): please file an issue in GLM.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 2, 2019, 7:44pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/9 "2019-02-02T19:44:04Z")

</div>

Hi!

> [@nalimilan](#):
>
> GLM.dispersion(lmobj.model, true)

It works, i have:

```julia
ulia> GLM.dispersion(ols.model, true)
0.006822235851285752

```

but in R i have:

```julia
> summary(lmobj )$sigma^2
[1] 0.006395846

```

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 3, 2019, 4:26pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/10 "2019-02-03T16:26:03Z")

</div>

For complience with R allowrankdeficient should be true:

`ols = lm(@formula(LnVar ~ Seq+Per+Trt+Subj), df, true);`

and for each categorical factor:

`categorical!(df, :Per);`

---

<div class="post-metadata">

**Author:** ![mkborregaard](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkborregaard/32/556_2.png) [@mkborregaard](https://discourse.julialang.org/u/mkborregaard)\
**Post date:** [February 3, 2019, 8:15pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/11 "2019-02-03T20:15:01Z")

</div>

Great that it got resolved so quickly.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 20, 2019, 4:22pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/12 "2019-02-20T16:22:23Z")

</div>

In described model i try to make:

```julia
ols = lm(@formula(LnVar ~ Seq+Per+Trt+Subj), df, true)
ols1 = lm(@formula(LnVar ~ Seq+Subj), df, true)
ols2 = lm(@formula(LnVar ~ Subj), df, true)
fts = ftest(ols1.model, ols2.model)

```

and i have:

LoadError: ArgumentError: FDist: the condition ν1 \> zero(ν1) && ν2 \> zero(ν2) is not satisfied.

And i try to make

```julia
 Anova(ols)

```

and have same result. (Pacage by [marcpabst](https://github.com/marcpabst))

Is a way to get Anova type III for crossover data?

---

<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:** [February 20, 2019, 5:45pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/13 "2019-02-20T17:45:06Z")

</div>

You need `ftest(ols2.model, ols1.model)`, right?

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 20, 2019, 6:13pm UTC](https://discourse.julialang.org/t/confidence-interval-for-certain-contrast-and-model-residual/20339/14 "2019-02-20T18:13:44Z")

</div>

Yes!