# Computing confidence limits for $\beta$ estimates on the reference level of a categorical variable

tl;dr I have some $\beta$ estimates from a logistic regression model. One of the predictors is a categorical variable with effects coding applied. I want to compute confidence limits for all levels of that variable, including the one with all contrasts set to -1. I am accustomed to using the profile methods of MASS::confint.glm(), but that function doesn't calculate limits for the reference level. I'm attempting to recreate the output of confint.glm() using a Z-approximation, but seem to be failing. I'm not sure where I'm going off the rails exactly, but suspect it has something to do with my selection of the appropriate $N$ to use for calculating the CI width. In fact, I've come fairly close to the profile limits by using $N=1$, but I would like to be sure that I'm making statistically-valid choices, not just finding numbers that match.

So my questions are:

1. Is $N=1$ the way to go here? If so, why is that the case? My intuition is that $N$ might be the number of cases falling into the particular category for which the estimate is being made.
2. If $N=1$ is not the way to go, where am I going wrong?

First, set up the environment and create some play data:

### Setup

library(MASS)         # for confint.glm()
library(plyr)         # for mapvalues()
library(tidyverse)    # for the forward-pipe and friends

set.seed(1138)
N <- 2000                                   # establish sample size
alpha <- 0.05                               # define desired confidence level
fruits <- c("papaya", "pomegranate",
"soursop", "guava", "durian")   # possible fruits
animals <- c("tribble", "ewok")             # experimental species (Empire- and Federation-approved)
eaten <- c("no", "yes")                     # possible treatments
location <- c("Tatooine", "Endor")          # environmental variable

# Craft the data
.df <- data.frame(ill = rbinom(n=N, size=1, prob=0.267),
fruit = fruits[ round( runif(n=N, min=1, max=length(fruits) ) ) ],
animal = animals[ round( runif(n=N, min=1, max=length(animals) ) ) ],
eaten = eaten[ round( runif(n=N, min=1, max=length(eaten) ) ) ],
planet = location[ round( runif(n=N, min=1, max=length(location) ) ) ] )

# Use effects coding/sum-to-zero contrasts on .df$fruit .df <- within( .df, { contrasts(fruit) <- contr.sum( levels(fruit) ) })  The data then look as follows—a binary response (ill), three binary predictors (animal, eaten, and planet), and one categorical predictor (fruit): ## ill fruit animal eaten planet ## 1 0 guava tribble yes Endor ## 2 0 papaya tribble no Tatooine ## 3 0 papaya tribble yes Endor ## 4 0 pomegranate ewok yes Tatooine ## 5 0 guava ewok yes Endor ## 6 0 soursop ewok yes Tatooine  Now I define the model... im1 <- glm( ill ~ fruit + animal + eaten + planet, data = .df, family = binomial(link="logit") )  ### Profile limits ...and calculate the odds ratios and 95% confidence limits using the profile method of MASS::confint.glm(). These are what I'm trying to replicate. im1.OR_profile <- exp( cbind( OR = coef(im1), # exponentiate the parameter-estimate coefficients... confint(im1) ) ) %>% # ...together with the conf limits to get odds-ratio estimates as.data.frame.matrix( . ) %>% mutate( Parameter = rownames( . ) ) %>% filter( grepl( "Intercept|fruit", Parameter ) ) %>% # for simplicity, I'm restricting this to the variable under sum-to-zero contrasts select( Parameter, everything() )  That gives me the following estimates and profile confidence limits: ## Parameter OR 2.5 % 97.5 % ## 1 (Intercept) 0.3082855 0.2492946 0.3794741 ## 2 fruit1 1.2874891 1.0121429 1.6273995 ## 3 fruit2 0.9159739 0.7543939 1.1079502 ## 4 fruit3 0.8640049 0.6672064 1.1066674 ## 5 fruit4 0.8450893 0.6919834 1.0273585  But it only gives me limits for four of the five levels of fruit. It's leaving out fruit5, the all–-1 contrasts level. I can calculate the$\beta$estimate for that level, but to get confidence limits I then need to use a Z approximation. So I try that: ### Z-approximation limits: beta <- coef(im1) # extract parameter estimates from the model beta[["fruit5"]] <- -sum( beta[ grepl("^fruit[1-4]", names(beta)) ] ) # compute beta estimate for reference level of 'fruit' vcmat <- vcov(im1) # extract variance-covariance matrix from model Rinds <- grepl("^fruit", rownames(vcmat)) # find relevant rows of vcmat (corresponding to levels of sum-to-zero factor) Cinds <- grepl("^fruit", colnames(vcmat)) # find relevant cols of vcmat (corresponding to levels sum-to-zero factor) vcmat_tiny <- vcmat[ Rinds, Cinds ] # restrict vcmat to relevant submatrix  This gives me the following variance-covariance (sub)matrix in vcmat_tiny: ## fruit1 fruit2 fruit3 fruit4 ## fruit1 0.014639932 -0.003104047 -0.005432924 -0.003282766 ## fruit2 -0.003104047 0.009599889 -0.003754821 -0.001598346 ## fruit3 -0.005432924 -0.003754821 0.016606235 -0.003944584 ## fruit4 -0.003282766 -0.001598346 -0.003944584 0.010146859  Next I compute the$\beta$estimate and standard deviation of the estimate for fruit5: beta.s2 <- diag(vcmat) # extract variances for each parameter from vcmat beta.s2[["fruit5"]] <- sum(vcmat_tiny) # compute variance of reference level as # Var(f1 + f2 + f3 + f4) = Var(f1) + Var(f2) + Var(f3) + Var(f4) + # 2Cov(f1,f2) + 2Cov(f1,f3) + 2Cov(f1,f4) + 2Cov(f2,f3) + # 2Cov(f2,f4) + 2Cov(f3,f4) beta.s <- sqrt(beta.s2) # calculate SD of parameter estimates as square root of variances beta <- beta[ grepl("Intercept|fruit", names(beta)) ] # for simplicity, I will limit this to the factor beta.s <- beta.s[ grepl("Intercept|fruit", names(beta.s)) ] # under sum-to-zero contrasts  This is where I start to really lose confidence in my calculations. I know (read: I think I know) that the half-width of the confidence interval is computed as$\bar{x}±z\frac{s}{\sqrt{n}}$. I have$\bar{x}$—that's the vector of$\beta$coefficients; I just calculated their standard deviations,$s$.$z$is an easy calculation. That leaves$n$. Now, is this$n=N$? That is, the size of the entire sample, 2000? Or, does each level get its own$n$corresponding to the number of cases in that level in the data? In fact, based on what I've experienced so far while trying to replicate the output of MASS::confint(), it seems that neither of these is the case. Instead, I get closest to the values produced by that function's methods if I don't divide$s$at all (effectively,$n=1$). I'll go through all cases here. Population$N$: z <- qnorm(1 - (alpha / 2)) # compute the critical z quantile corresponding to the desired confidence level CIWidth.N <- z * beta.s / sqrt(N) # compute the half-width of the confidence interval (N=2000 for all parameters) CIlo.N <- beta - CIWidth.N # calculate lower limit of CI CIhi.N <- beta + CIWidth.N # calculate upper limit of CI im1.OR_popN <- exp( cbind( OR = beta, # put it all together and exponentiate 2.5 % = CIlo.N, 97.5 % = CIhi.N ) ) %>% as.data.frame.matrix( . ) %>% mutate( Parameter = rownames( . ) ) %>% select( Parameter, everything() )  Comparing these estimates + limits to those computed earlier I have: ## Parameter OR 2.5 % 97.5 % Parameter OR 2.5 % 97.5 % ## 1 (Intercept) 0.3082855 0.2492946 0.3794741 (Intercept) 0.3082855 0.3068414 0.3097363 ## 2 fruit1 1.2874891 1.0121429 1.6273995 fruit1 1.2874891 1.2806799 1.2943345 ## 3 fruit2 0.9159739 0.7543939 1.1079502 fruit2 0.9159739 0.9120491 0.9199156 ## 4 fruit3 0.8640049 0.6672064 1.1066674 fruit3 0.8640049 0.8591390 0.8688983 ## 5 fruit4 0.8450893 0.6919834 1.0273585 fruit4 0.8450893 0.8413667 0.8488284 ## 6 <NA> NA NA NA fruit5 1.1613271 1.1565738 1.1661000  Here, the original profile (confint()) estimates and limits are shown on the left; the corresponding Z-approximation ($n=N$) limits are on the right. Clearly, the latter intervals are far too small;$N$is too big. Maybe the more appropriate$n$to use is the$n_i$in each level of the fruit factor? Level$n_i$: n_i <- table(.df$fruit)                                    # determine the number of cases for each level of fruit
fruitlevels <- c(fruit1="durian", fruit2="guava",          # draft a map of factor levels --> contrast labels
fruit3="papaya", fruit4="pomegranate",
fruit5="soursop" )

names(n_i) <- mapvalues( names(n_i),                       # rename the table of cases according to that map
fruitlevels, names(fruitlevels) )

n_i[["(Intercept)"]] <- mean(n_i)                          # add an n for the 'average' case as the 'average'
n_i <- n_i[ match( names(beta.s), names(n_i) ) ]           #     of the ns????; rearrange to match beta.s

CIWidth.n_i <- z * beta.s / sqrt(n_i)                      # compute the CI half-width using level-specific ns

CIlo.n_i <- beta - CIWidth.n_i                             # calculate the lower confidence limit
CIhi.n_i <- beta + CIWidth.n_i                             # calculate the upper confidence limit

im1.OR_n_i <- exp( cbind( OR = beta,                       # put everything together and exponentiate
2.5 % = CIlo.n_i,
97.5 % = CIhi.n_i ) ) %>%
as.data.frame.matrix( . ) %>%
mutate( Parameter = rownames( . ) ) %>%
select( Parameter, everything() )


Again, I'll compare these to the limits estimated using profile methods:

##     Parameter        OR     2.5 %    97.5 %   Parameter        OR     2.5 %    97.5 %
## 1 (Intercept) 0.3082855 0.2492946 0.3794741 (Intercept) 0.3082855 0.3050658 0.3115392
## 2      fruit1 1.2874891 1.0121429 1.6273995      fruit1 1.2874891 1.2680903 1.3071846
## 3      fruit2 0.9159739 0.7543939 1.1079502      fruit2 0.9159739 0.9081644 0.9238505
## 4      fruit3 0.8640049 0.6672064 1.1066674      fruit3 0.8640049 0.8504203 0.8778065
## 5      fruit4 0.8450893 0.6919834 1.0273585      fruit4 0.8450893 0.8375471 0.8526995
## 6        <NA>        NA        NA        NA      fruit5 1.1613271 1.1519695 1.1707607


Still too narrow! Again, the confint() limits are on the left; the corresponding $n=n_i$ limits are on the right.

I guess in the extreme case, $n=1$. I'll try that.

Singular $n$:
CIWidth.1 <- z * beta.s / sqrt(1)                  # compute the half-width of the confidence interval using n=1

CI.lo1 <- beta - CIWidth.1                         # calculate the lower confidence limit
CI.hi1 <- beta + CIWidth.1                         # calculate the upper confidence limit

im1.OR_1 <- exp( cbind( OR = beta,                 # put everything together and exponentiate
2.5 % = CI.lo1,
97.5 % = CI.hi1 ) ) %>%
as.data.frame.matrix( . ) %>%
mutate( Parameter = rownames( . ) ) %>%
select( Parameter, everything() )


One last time! Compare:

##     Parameter        OR     2.5 %    97.5 %   Parameter        OR     2.5 %   97.5 %
## 1 (Intercept) 0.3082855 0.2492946 0.3794741 (Intercept) 0.3082855 0.2498973 0.380316
## 2      fruit1 1.2874891 1.0121429 1.6273995      fruit1 1.2874891 1.0156684 1.632056
## 3      fruit2 0.9159739 0.7543939 1.1079502      fruit2 0.9159739 0.7559328 1.109898
## 4      fruit3 0.8640049 0.6672064 1.1066674      fruit3 0.8640049 0.6711598 1.112260
## 5      fruit4 0.8450893 0.6919834 1.0273585      fruit4 0.8450893 0.6936808 1.029545
## 6        <NA>        NA        NA        NA      fruit5 1.1613271 0.9667092 1.395126


Much closer; now the CIs are just barely too wide, and appear to be slightly up-shifted.

To reiterate my questions:

1. Is $N=1$ the way to go here? If so, why is that the case? My intuition is that $N$ might be the number of cases falling into the particular category for which the estimate is being made.
2. If $N=1$ is not the way to go, where am I going wrong?