# Confidence Interval for certain Contrast and model residual

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”?

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 you would use `stderror` to get the standard errors of your coefficients.

2 Likes

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.

A have a dataset 12248_2014_9661_MOESM1_ESM.txt from here:

and i try to validate Julia code for BE studies:

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.

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

2 Likes

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?

1 Like

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.

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.

1 Like

Hi!

It works, i have:

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

but in R i have:

``````> summary(lmobj )\$sigma^2
 0.006395846
``````

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);`

Great that it got resolved so quickly.

In described model i try to make:

``````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

`````` Anova(ols)
``````

and have same result. (Pacage by marcpabst)

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

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

Yes!