0
$\begingroup$

The code

confint(model, level = 1 -0.15/length(coef(model)))

was recommended by https://stat.ethz.ch/pipermail/r-help/2005-March/067570.html

However, I have two questions:

  1. How do I avoid including the intercept? For example, here is my output:

    > confint(ACFL.m1, level = 1 -0.15/length(coef(ACFL.m1)))
                   1.87 %    98.12 %
    (Intercept)   1.5838374  3.2846837
    PDCD          8.6257799 12.3346769
    PMIX        -17.6815440  6.8408247
    RTCL          0.2888126  0.6987631
    

    I only want 85% CI for my predictors and not the intercept?

  2. In the above example, level = 1-0.15/4 = 0.9625. I know, to get 85% CI for each predictor, the individual CI has to be > 85%, but how do I understand 0.9625% for each summing to as an 85% simultaneous CI?

$\endgroup$
  • $\begingroup$ I don't understand your last question. To get an interval of a certain width use confint(fit, level=mylevel). Is that what you're asking? $\endgroup$ – AdamO Mar 26 '18 at 18:30
  • $\begingroup$ Based on the formula above, my 85% family-wise CI = 1-0.15/4 = 0.9625% for each variable. I now understand how the lower (1.87 =(1-0.9625)/2) and upper bound values (98.12= ((1-0.9625)/2)+0.9625) came about. I had made a fundamental error in assuming that the 96.25% was a sum of CIs for each variable. Now I understand that each variable is tested at 96.25% so that the family-wise error is 85%. I have left the question as it is so others with similar questions can benefit. Thank you. $\endgroup$ – Guphadi Mar 27 '18 at 17:48
  • $\begingroup$ could you include your last comment as a paragraph in your answer below? $\endgroup$ – Ben Bolker Mar 28 '18 at 16:44
1
$\begingroup$

The answer to question 1:

I learned that the package multcomp allows calculation of simultaneous CIs without the intercept. Details about the package are in https://cran.r-project.org/web/packages/multcomp/multcomp.pdf

Its' default is the single-step method and is considered to be more powerful than the Bonferroni method (Bretz et al. 2011). The results of the default method is shown below:

> K <- cbind(0, diag(length(coef(ACFL.m1))-1))
> rownames(K) <- names(coef(ACFL.m1)) [-1]
> ACLF.mc <- glht(ACFL.m1, linfct = K)
> confint(ACLF.mc, level = 0.85)

Simultaneous Confidence Intervals

Fit: lm(formula = PINDX ~ PDCD + PMIX + RTCL, data = ACFL)

Quantile = 1.9021
85% family-wise confidence level


Linear Hypotheses:
           Estimate   lwr      upr     
PDCD == 0  10.4802   8.7910  12.1695
PMIX == 0  -5.4204 -16.5895   5.7488
RTCL == 0   0.4938   0.3071   0.6805

To get the Bonferroni 85% CI, I had to first run the Bonferroni Test

ACFL.bf <- summary(ACLF.mc, test= adjusted(type = "bonferroni")) 

Then calculate the confidence interval

confint(ACFL.bf, level = 0.85)

I got the following output:

 Simultaneous Confidence Intervals

Fit: lm(formula = PINDX ~ PDCD + PMIX + RTCL, data = ACFL)

Quantile = 1.9021
85% family-wise confidence level


Linear Hypotheses:
           Estimate   lwr      upr     
PDCD == 0  10.4802   8.7909  12.1696
PMIX == 0  -5.4204 -16.5898   5.7491
RTCL == 0   0.4938   0.3071   0.6805

Good reference for the multcomp package:

Bretz F, Hothorn T, Westfall PH (2011) Multiple comparisons using R. CRC Press, Boca Raton, FL

Hothorn T (2017) Additional multcomp examples. Accessed online form the R project website https://cran.r-project.org/web/packages/multcomp/vignettes/multcomp-examples.pdf

Hothorn T, Bretz F, Westfall P (2008) Simultaneous inference in general parametric models. Biom J 50:346–363

The answer to question 2:

Based on the formula above, my 85% family-wise CI = 1-0.15/4 = 0.9625% for each variable. I now understand how the lower (1.87 =(1-0.9625)/2) and upper bound values (98.12= ((1-0.9625)/2)+0.9625) came about. I had made a fundamental error in assuming that the 96.25% was a sum of CIs for each variable. I should not have been confused between the 'Bonferroni' and the 'sequential Bonferroni' (Holm, 1979). Now I understand that each variable is tested at 96.25% so that the family-wise error is 85%. I have left the question as it is so others with similar questions can benefit.

$\endgroup$
  • $\begingroup$ To appropriate a quote from a famous painter whom I admire, "You are a professional statistician, you don't need no stinkin' package". To get the right interval width, just use 3 instead of 4 to get the family-wise corrected alpha level for the individual tests, and do the intervals based on that. You're done. The output is nice, but the conceptual shortcoming was not realizing that the intercept didn't need a test at all. $\endgroup$ – AdamO Mar 28 '18 at 16:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.