Logistic Regression in R

In binary logistic regression, there is a single binary dependent variable, coded by an indicator variable. For example, if we represent a response as 1 and non-response as 0, then the corresponding probability of response, can be between 0 (certainly not a response) and 1 (certainly a response) - hence the labeling !

The logistic model models the log-odds of an event as a linear combination of one or more independent variables (explanatory variables). If we observed $(y_i, x_i),$ where $y_i$ is a Bernoulli variable and $x_i$ a vector of explanatory variables, the model for $\pi_i = P(y_i=1)$ is

\[ \text{logit}(\pi_i)= \log\left\{ \frac{\pi_i}{1-\pi_i}\right\} = \beta_0 + \beta x_i, i = 1,\ldots,n \]

The model is especially useful in case-control studies and leads to the effect of risk factors by odds ratios.

Example: Lung Cancer Data

Data source: Loprinzi CL. Laurie JA. Wieand HS. Krook JE. Novotny PJ. Kugler JW. Bartel J. Law M. Bateman M. Klatt NE. et al. Prospective evaluation of prognostic variables from patient-completed questionnaires. North Central Cancer Treatment Group. Journal of Clinical Oncology. 12(3):601-7, 1994.

wgt_catn consists of: 1= patients a weight loss of zero or less, 0= patients with a weight loss of more than zero

trt01pn consists of 1= active treatment, 0 = placebo

Model Fit

We analyze the event of weight gain (or staying the same weight) in lung cancer patients in dependency of treatment (active or placebo), age, sex, ECOG performance score and calories consumed at meals. One of the most important things to remember is to ensure you tell R what your event is and what treatment comparison you are doing Active / Placebo or Placebo/Active! The easiest way to do this is to have event (or non-reference treatment) as 1, and non-event (reference treatment) as 0.

Below we are using wt_catn (0,1) and trt01pn (1,2) and sex (1,2). Let’s see what happens !

head(lung)

lung<-read.csv("../data/lung_cancer.csv")
head(lung)

  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss trt01p
1    3  306      2  74   1       1       90       100     1175      NA Active
2    3  455      2  68   1       0       90        90     1225      15 Active
3    3 1010      1  56   1       0       90        90       NA      15 Active
4    5  210      2  57   1       1       90        60     1150      11 Active
5    1  883      2  60   1       0      100        90       NA       0 Active
6   12 1022      1  74   1       1       50        80      513       0 Active
  trt01pn dose_mg dose_id wt_cat wt_catn cnsr
1       1      10       1   <NA>      NA    0
2       1      10       1   loss       0    0
3       1      10       1   loss       0    1
4       1      10       1   loss       0    0
5       1      10       1   gain       1    0
6       1      10       1   gain       1    1

m1 <- glm(wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
          data = lung, 
          family = binomial(link="logit"))
summary(m1)


Call:
glm(formula = wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
    family = binomial(link = "logit"), data = lung)

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.8623998  1.7675776  -2.185   0.0289 *
trt01pn      0.3887667  0.3781566   1.028   0.3039  
age          0.0122549  0.0211553   0.579   0.5624  
sex          0.8321005  0.3743793   2.223   0.0262 *
ph.ecog     -0.3763592  0.2638322  -1.427   0.1537  
meal.cal     0.0008500  0.0004486   1.895   0.0581 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.36  on 169  degrees of freedom
Residual deviance: 190.46  on 164  degrees of freedom
  (58 observations deleted due to missingness)
AIC: 202.46

Number of Fisher Scoring iterations: 4

R by default sets the first category as the baseline category, hence trt01pn =1 and sex =1 are the baseline level, and other levels of the variable are contrasted to this level. This is using contr.treatment option (more information on this later!). The estimate for those variables (0.3887 and 0.8321) are the increase in the log-odds of being of receiving treatment 2 vs 1, and of being of sex=2 vs 1. The exponential of the estimate is the odds ratio. For example, exp(0.3887)=1.475, hence Treatment 2 is 1.475 times as likely to have weight gain compared to Treatment 1.

The intercept represents the baseline log odds of the outcome when all predictor variables are set to zero. In the above model, we have variables treatment and sex, each coded as 1 and 2. Currently R is not treating these as binary factors, instead R thinks they are a continuous variable, and hence the intercept is where treatment=0 and sex=0. In the context of our model, this doesn’t make sense as you can’t have a zero gender and zero treatment (but it may not matter as we rarely look at the intercept term anyway!)

However, if you want to have an interpretable intercept term (and if you want to match SAS output for intercept!), then it’s important to ensure any factors in your model as fitted as such in R. You can do this using: lung$trt01pn<-as.factor(lung$trt01pn) or by changing the variable to a 0 and 1. Note if you are using the same contr.treatment option, then this only affects the intercept estimate not the variable estimates.

Modelling factors correctly and Interpretation of the Intercept

So far we’ve learnt that it is good practice to always have binary and categorical variables set as factors, and that we should specify what method we want R to use for doing any contrasts (e.g. contr.treatment). You can specify different contrast methods for each variable by including them in a list in the model. It can sometimes be hard to see what contrast method R is using, so best to always specify this in your model (e.g. contrasts = list(trt01pn = "contr.treatment", sex="contr.sum"). It is helpful to view what the contrasts look like before you select which to use.

A factor with 2 levels (2 treatments) using contr.treatment would set the first level to be the reference (0) and contrast the second level to the baseline.

A factor with 4 levels (4 treatments) using contr.treatment would set the first level to be the reference (0), you would see 3 parameters in the model with 3 estimates corresponding to the increase in log-odds attributable to contrasting the second, third or fourth level to the baseline.

contr.treatment(2)

  2
1 0
2 1

contr.treatment(4)

Below we apply contr.treatment to our trt01pn and sex variables. Note: how the variables trt01pn2 and sex2 are now shown in the output, this is indicating that these rows relates to treatment=2 and sex=2. See the next section for how to interpret the estimates and change this contr. option.

# Good practice to have binary variables (eg. treatment) identified as a factor 
# And to specify what contrast option you are using ! more on this below
lung$trt01pn<-as.factor(lung$trt01pn)
lung$sex<-as.factor(lung$sex)

m1 <- glm(wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
          data = lung, 
          family = binomial(link="logit"),
          contrasts = list(trt01pn = "contr.treatment", sex = "contr.treatment"))
summary(m1)


Call:
glm(formula = wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
    family = binomial(link = "logit"), data = lung, contrasts = list(trt01pn = "contr.treatment", 
        sex = "contr.treatment"))

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -2.6415326  1.5140191  -1.745   0.0810 .
trt01pn2     0.3887667  0.3781566   1.028   0.3039  
age          0.0122549  0.0211553   0.579   0.5624  
sex2         0.8321005  0.3743793   2.223   0.0262 *
ph.ecog     -0.3763592  0.2638322  -1.427   0.1537  
meal.cal     0.0008500  0.0004486   1.895   0.0581 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.36  on 169  degrees of freedom
Residual deviance: 190.46  on 164  degrees of freedom
  (58 observations deleted due to missingness)
AIC: 202.46

Number of Fisher Scoring iterations: 4

Model parameter estimates

The model summary contains the parameter estimates $\beta_j$ for each explanatory variable $x_j$, corresponding to the log-odds for the response variable to take the value $1$, conditional on all other explanatory variables remaining constant. For better interpretation, we can exponentiate these estimates, to obtain estimates for the odds instead and provide 95% confidence intervals.

How R parameterizes your variables in the model

contr.treatment - [default] sets the first level to be the reference and contrasts each other level with the baseline. For example, based on the model shown below. Log-odds for Treatment 2 = -2.64153 + 0.38876 =-2.25277 , Log-odds for Treatment 1=-2.64153
contr.sum - sets the last level to be the reference and compares the mean of the dependent variable for a given level to the overall mean of the dependent variable. For treatment 1, the intercept + trt01pn1 estimate. For Treatment 2, the intercept - trt01pn1 estimate.
Log-odds for Treatment 2 = -2.44715 - - 0.19438 = -2.25277, Log-odds for Treatment 1=-2.44715 - 0.19438 = -2.64153
use contr.sum(2) to show the contrast you are using and hence how to interpret your estimates.
exponential (ratio log-odds) = odds ratio. eg. -2.25277/ -2.64153 = 0.85283 Treatment 2 is 0.85 times as likely (eg. 15% less likely) to have weight gain compared to Treatment 1.

ma <- glm(wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
          data = lung, 
          family = binomial(link="logit"),
          contrasts = list(trt01pn = "contr.treatment", sex = "contr.treatment"))
summary(ma)


Call:
glm(formula = wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
    family = binomial(link = "logit"), data = lung, contrasts = list(trt01pn = "contr.treatment", 
        sex = "contr.treatment"))

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -2.6415326  1.5140191  -1.745   0.0810 .
trt01pn2     0.3887667  0.3781566   1.028   0.3039  
age          0.0122549  0.0211553   0.579   0.5624  
sex2         0.8321005  0.3743793   2.223   0.0262 *
ph.ecog     -0.3763592  0.2638322  -1.427   0.1537  
meal.cal     0.0008500  0.0004486   1.895   0.0581 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.36  on 169  degrees of freedom
Residual deviance: 190.46  on 164  degrees of freedom
  (58 observations deleted due to missingness)
AIC: 202.46

Number of Fisher Scoring iterations: 4

contr.treatment(2)

  2
1 0
2 1

mb <- glm(wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
          data = lung, 
          family = binomial(link="logit"),
          contrasts = list(trt01pn = "contr.sum", sex = "contr.treatment"))
summary(mb)


Call:
glm(formula = wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
    family = binomial(link = "logit"), data = lung, contrasts = list(trt01pn = "contr.sum", 
        sex = "contr.treatment"))

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -2.4471493  1.4929863  -1.639   0.1012  
trt01pn1    -0.1943833  0.1890783  -1.028   0.3039  
age          0.0122549  0.0211553   0.579   0.5624  
sex2         0.8321005  0.3743793   2.223   0.0262 *
ph.ecog     -0.3763592  0.2638322  -1.427   0.1537  
meal.cal     0.0008500  0.0004486   1.895   0.0581 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.36  on 169  degrees of freedom
Residual deviance: 190.46  on 164  degrees of freedom
  (58 observations deleted due to missingness)
AIC: 202.46

Number of Fisher Scoring iterations: 4

contr.sum(2)

  [,1]
1    1
2   -1

Calculation of confidence intervals

Using the above model, you can output the estimates and confidence intervals using coef() and confint() and exponential back transforming using exp(). NOTE: that there are two types of confidence intervals that you can calculate. Function confint.default gives the Wald confidence limits, which is the default option in SAS PROC LOGISTIC procedure; whereas confint gives the profile-likelihood limits.

# model coefficients summary
summary(m1)$coefficients

                 Estimate   Std. Error    z value   Pr(>|z|)
(Intercept) -2.6415326252 1.5140190574 -1.7447156 0.08103439
trt01pn2     0.3887666677 0.3781565596  1.0280574 0.30392281
age          0.0122549015 0.0211552875  0.5792831 0.56239813
sex2         0.8321005169 0.3743792762  2.2226137 0.02624186
ph.ecog     -0.3763592487 0.2638321918 -1.4265100 0.15372119
meal.cal     0.0008499918 0.0004486401  1.8945961 0.05814593

## Wald confidence limits
cbind(est = exp(coef(m1)), 
          exp(confint.default(m1)))

                   est       2.5 %   97.5 %
(Intercept) 0.07125198 0.003664896 1.385263
trt01pn2    1.47516031 0.702994248 3.095470
age         1.01233030 0.971213751 1.055188
sex2        2.29814096 1.103327506 4.786840
ph.ecog     0.68635572 0.409236994 1.151128
meal.cal    1.00085035 0.999970674 1.001731

## profile-likelihood limits
cbind(est = exp(coef(m1)), 
          exp(confint(m1)))

Waiting for profiling to be done...

                   est       2.5 %   97.5 %
(Intercept) 0.07125198 0.003312288 1.302587
trt01pn2    1.47516031 0.696210092 3.085089
age         1.01233030 0.971670916 1.056194
sex2        2.29814096 1.107651762 4.836770
ph.ecog     0.68635572 0.405659156 1.147452
meal.cal    1.00085035 0.999978126 1.001761

Comparing 2 models

To compare two logistic models, the residual deviances (-2 * log likelihoods) are compared against a $\chi^2$-distribution with degrees of freedom calculated using the difference in the two models’ parameters. Below, the only difference is the inclusion/exclusion of age in the model, hence we test using $\chi^2$ with 1 df. Here testing at the 5% level.

m2 <- glm(wt_catn ~ trt01pn + sex + ph.ecog + meal.cal, 
                    data = lung, 
                    family = binomial(link="logit"),
                    contrasts = list(trt01pn = "contr.treatment"))
summary(m2)


Call:
glm(formula = wt_catn ~ trt01pn + sex + ph.ecog + meal.cal, family = binomial(link = "logit"), 
    data = lung, contrasts = list(trt01pn = "contr.treatment"))

Coefficients:
              Estimate Std. Error z value Pr(>|z|)   
(Intercept) -1.8350766  0.5810727  -3.158  0.00159 **
trt01pn2     0.3681563  0.3761976   0.979  0.32777   
sex2         0.7919227  0.3674936   2.155  0.03117 * 
ph.ecog     -0.3312031  0.2527586  -1.310  0.19008   
meal.cal     0.0007896  0.0004378   1.803  0.07131 . 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.36  on 169  degrees of freedom
Residual deviance: 190.80  on 165  degrees of freedom
  (58 observations deleted due to missingness)
AIC: 200.8

Number of Fisher Scoring iterations: 4

anova(m1, m2, test = "LRT")

Analysis of Deviance Table

Model 1: wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal
Model 2: wt_catn ~ trt01pn + sex + ph.ecog + meal.cal
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1       164     190.46                     
2       165     190.80 -1 -0.33867   0.5606

Stackexchange here has a good article describing this method and the difference between comparing 2 models using the likelihood ratio tests versus using wald tests and Pr>chisq (from the maximum likelihood estimate). Note: anova(m1, m2, test = "Chisq") and using test="LRT" as above are synonymous in this context.

Predicting likelihood of response for new patients

Predictions from the model for the log-odds of a patient with new data to experience a weight loss are derived using predict():

# new female, symptomatic but completely ambulatory patient consuming 2500 calories
new_pt <- data.frame(trt01pn=1, age=48, sex=2, ph.ecog=1, meal.cal=2500)
new_pt$trt01pn<-as.factor(new_pt$trt01pn)
new_pt$sex<-as.factor(new_pt$sex)
predict(m1, new_pt, type = "response")

       1 
0.628882

Creating Treatment Contrasts for 2 treatments

{emmeans}

Here we will use {emmeans} to output the log-odds of weight gain for treatment 1 and treatment 2.

NOTE as per the output, these are on the logit scale, you need to exponentiate to get the odds (or use the type=“response” option).

The treatment comparison can also be output on the log-odds or back transformed scale as shown below.

m3 <- glm(wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
          data = lung, 
          family = binomial(link="logit"),
          contrasts = list(trt01pn = "contr.treatment"))
summary(m3)


Call:
glm(formula = wt_catn ~ trt01pn + age + sex + ph.ecog + meal.cal, 
    family = binomial(link = "logit"), data = lung, contrasts = list(trt01pn = "contr.treatment"))

Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept) -2.6415326  1.5140191  -1.745   0.0810 .
trt01pn2     0.3887667  0.3781566   1.028   0.3039  
age          0.0122549  0.0211553   0.579   0.5624  
sex2         0.8321005  0.3743793   2.223   0.0262 *
ph.ecog     -0.3763592  0.2638322  -1.427   0.1537  
meal.cal     0.0008500  0.0004486   1.895   0.0581 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 202.36  on 169  degrees of freedom
Residual deviance: 190.46  on 164  degrees of freedom
  (58 observations deleted due to missingness)
AIC: 202.46

Number of Fisher Scoring iterations: 4

contr.treatment(2)

  2
1 0
2 1

#log-odds for each treatment
lsm<-emmeans(m3,"trt01pn")
lsm

 trt01pn emmean    SE  df asymp.LCL asymp.UCL
 1       -1.034 0.225 Inf     -1.47   -0.5935
 2       -0.645 0.302 Inf     -1.24   -0.0529

Results are averaged over the levels of: sex 
Results are given on the logit (not the response) scale. 
Confidence level used: 0.95

#log-odds ratios (treatment comparison): This does all pairwise comparisons
# However as seen below, this is TRT 1 - TRT 2
pairs(lsm)

 contrast            estimate    SE  df z.ratio p.value
 trt01pn1 - trt01pn2   -0.389 0.378 Inf  -1.028  0.3039

Results are averaged over the levels of: sex 
Results are given on the log odds ratio (not the response) scale.

# the below creates tests and CI's prior to back transformation (ratios of geometric means)
pairs(lsm, type="response")

 contrast            odds.ratio    SE  df null z.ratio p.value
 trt01pn1 / trt01pn2      0.678 0.256 Inf    1  -1.028  0.3039

Results are averaged over the levels of: sex 
Tests are performed on the log odds ratio scale

# see coefficients of the linear functions
coef(pairs(lsm))

         trt01pn c.1
trt01pn1       1   1
trt01pn2       2  -1

# Output treatment contrasts 2 vs 1 and 95% CIs, the type="response" option back transforms the results
trtdiff<-contrast(lsm,"poly")
trtdiff

 contrast estimate    SE  df z.ratio p.value
 linear      0.389 0.378 Inf   1.028  0.3039

Results are averaged over the levels of: sex 
Results are given on the log odds ratio (not the response) scale.

confint(trtdiff, type="response")

 contrast odds.ratio    SE  df asymp.LCL asymp.UCL
 linear         1.48 0.558 Inf     0.703       3.1

Results are averaged over the levels of: sex 
Confidence level used: 0.95 
Intervals are back-transformed from the log odds ratio scale

In Summary: Treatment 2 is on average 1.48 times as likely to have weight gain compared to treatment 1, however this is not statistically significant (95% Confidence interval = 0.703-3.100, p-value= 0.3039).

Creating Treatment Contrasts for 2 or more treatments

{emmeans} can also be used to do specific contrasts, instead of investigating treatment (active vs placebo), suppose we now want to look at 2 dose groups vs placebo.

We have a 3 level treatment variable (dose_id), where 1=10mg, 2=20mg doses for active treatment and 3= placebo which is 0mg. We want to test the null hypothesis that 0.5*dose10mg + 0.5*dose20mg - placebo = 0. That there is not difference between the doses.

lung2<-lung |> 
    mutate(dose_id2 = as.factor(lung$dose_id))

m3 <- glm (wt_catn ~ dose_id2 + age + sex + ph.ecog + meal.cal, 
                     data = lung2, 
                     family = binomial(link="logit"),
                     contrasts = list(dose_id2 = "contr.treatment"))
#log-odds for each treatment
lsm3<-emmeans(m3,"dose_id2")
lsm3

 dose_id2 emmean    SE  df asymp.LCL asymp.UCL
 1        -0.903 0.296 Inf     -1.48   -0.3228
 2        -1.201 0.348 Inf     -1.88   -0.5201
 3        -0.643 0.302 Inf     -1.23   -0.0518

Results are averaged over the levels of: sex 
Results are given on the logit (not the response) scale. 
Confidence level used: 0.95

contrast (lsm3,list(AveDose_vs_pbo=c(0.5,0.5,-1)))

 contrast       estimate   SE  df z.ratio p.value
 AveDose_vs_pbo    -0.41 0.38 Inf  -1.078  0.2813

Results are averaged over the levels of: sex 
Results are given on the log odds ratio (not the response) scale.

confint(contrast (lsm3,list(AveDose_vs_pbo=c(0.5,0.5,-1))))

 contrast       estimate   SE  df asymp.LCL asymp.UCL
 AveDose_vs_pbo    -0.41 0.38 Inf     -1.15     0.335

Results are averaged over the levels of: sex 
Results are given on the log odds ratio (not the response) scale. 
Confidence level used: 0.95

Here we found that on average there is -0.41 times the risk of weight gain on active vs placebo but that this is not statisticallly significantly different (95% CI -1.15 to 0.335, p-value = 0.2813).

See the emmeans vignette on creating bespoke contrasts here.

{gmodels}

{gmodels} is an alternative package to create contrasts instead of <emmeans>, you can use the fit.contrast() function from the gmodels package. The same result is obtained as <emmeans>.

lung2<-lung |> 
    mutate(dose_id2 = as.factor(lung$dose_id))

m3 <- glm (wt_catn ~ dose_id2 + age + sex + ph.ecog + meal.cal, 
                     data = lung2, 
                     family = binomial(link="logit"),
                     contrasts = list(dose_id2 = "contr.treatment"))

fit.contrast (m3,'dose_id2',c(0.5,0.5,-1),conf.int=0.95)

                            Estimate Std. Error   z value  Pr(>|z|)  lower CI
dose_id2 c=( 0.5 0.5 -1 ) -0.4096323  0.3801683 -1.077502 0.2812558 -1.160322
                           upper CI
dose_id2 c=( 0.5 0.5 -1 ) 0.3410574

Reference

#\| echo: false

#List all the packages needed

si <- sessioninfo::session_info(c('dplyr','emmeans','gmodels')) 
si

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       macOS Sequoia 15.6
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       Europe/London
 date     2025-08-20
 pandoc   3.4 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 ! package      * version    date (UTC) lib source
 P cli            3.6.3      2024-06-21 [?] CRAN (R 4.4.0)
 P dplyr        * 1.1.4      2023-11-17 [?] CRAN (R 4.4.0)
 P emmeans      * 1.10.4     2024-08-21 [?] CRAN (R 4.4.1)
 P estimability   1.5.1      2024-05-12 [?] CRAN (R 4.4.0)
 P fansi          1.0.6      2023-12-08 [?] CRAN (R 4.4.0)
 P gdata          3.0.1      2024-10-22 [?] RSPM
 P generics       0.1.3      2022-07-05 [?] CRAN (R 4.4.0)
 P glue           1.8.0      2024-09-30 [?] CRAN (R 4.4.1)
 P gmodels      * 2.19.1     2024-03-06 [?] RSPM
 P gtools         3.9.5      2023-11-20 [?] RSPM
 P lifecycle      1.0.4      2023-11-07 [?] CRAN (R 4.4.0)
 P magrittr       2.0.3      2022-03-30 [?] CRAN (R 4.4.0)
   MASS           7.3-61     2024-06-13 [2] CRAN (R 4.4.2)
 P mvtnorm        1.3-1      2024-09-03 [?] CRAN (R 4.4.1)
   numDeriv       2016.8-1.1 2019-06-06 [1] CRAN (R 4.4.0)
 P pillar         1.9.0      2023-03-22 [?] CRAN (R 4.4.0)
 P pkgconfig      2.0.3      2019-09-22 [?] CRAN (R 4.4.0)
 P R6             2.5.1      2021-08-19 [?] CRAN (R 4.4.0)
 P rlang          1.1.4      2024-06-04 [?] CRAN (R 4.4.0)
 P tibble         3.2.1      2023-03-20 [?] CRAN (R 4.4.0)
 P tidyselect     1.2.1      2024-03-11 [?] CRAN (R 4.4.0)
 P utf8           1.2.4      2023-10-22 [?] CRAN (R 4.4.0)
 P vctrs          0.6.5      2023-12-01 [?] CRAN (R 4.4.0)
   withr          3.0.1      2024-07-31 [1] CRAN (R 4.4.0)

 [1] /Users/christinafillmore/Documents/GitHub/CAMIS/renv/library/macos/R-4.4/aarch64-apple-darwin20
 [2] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library

 P ── Loaded and on-disk path mismatch.

──────────────────────────────────────────────────────────────────────────────