Confidence Interval for certain Contrast and model residual

Specific Domains Statistics
1

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

2

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
3

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.

4

A have a dataset 12248_2014_9661_MOESM1_ESM.txt from here:

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.

5

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
6

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
7

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.

8

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
9

Hi!

It works, i have:

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

but in R i have:

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

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

11

Great that it got resolved so quickly.

12

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?

13

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

14

Yes!