GLM model runs in interactive code but not when I use knitr

I was able to reproduce the problem you describe, but haven't yet fully understood what happens under the hood.
This Markdown chunck is interesting :

```{r}
scabiesrisk_OR <- exp(cbind(OR= coef(scabiesrisk), confint((scabiesrisk))))
scabiesrisk_summary <- summary(scabiesrisk)
scabiesrisk_summary <- cbind(scabiesrisk_OR, scabiesrisk_summary$coefficients)
scabiesrisk_summary
```

If I manually quickly execute the lines in the chunck one after another (ctrl+Enter x 4), sometimes I get two profiling messages:

Waiting for profiling to be done...
Waiting for profiling to be done...

In this case, summary(scabiesrisk) is a matrix:

> class(scabiesrisk_summary)
[1] "matrix" "array" 

If I manually slowly execute the lines in the chunk, I get only one profiling message:

Waiting for profiling to be done...

summary(scabiesrisk) is a summary.glm :

> class(scabiesrisk_summary)
[1] "summary.glm"

Looks like profiling is launched on a separate thread, and depending on whether it was finished or not, summary function doesn't have the same behaviour. If profiling is finished, it returns the expected summary.glm object, but if it isn't the case it launches another profiling and returns a matrix.
In particular, with a matrix scabiesrisk_summary$coefficients isn't available and I get in this situation the following error message:

Error in scabiesrisk_summary$coefficients : 
  $ operator is invalid for atomic vectors

This could possibly also happen while knitting : does knitting overhead make profiling slower so that the problem occurs?

With the workaround found here (use confint.defaultinstead of confint), I wasn't able to reproduce the above problem:

scabiesrisk_OR <- exp(cbind(OR= coef(scabiesrisk), confint.default((scabiesrisk))))
scabiesrisk_summary <- summary(scabiesrisk)
scabiesrisk_summary <- cbind(scabiesrisk_OR, scabiesrisk_summary$coefficients)
scabiesrisk_summary
                       OR      2.5 %    97.5 %   Estimate Std. Error
(Intercept)    0.09357141 0.06984512 0.1253575 -2.3690303  0.1492092
agegroups0-10  2.20016940 1.60953741 3.0075383  0.7885344  0.1594864
agegroups11-20 2.53291768 1.79985894 3.5645415  0.9293719  0.1743214
gendermale     1.44749159 1.13922803 1.8391682  0.3698321  0.1221866
house_cat6-10  1.30521927 1.02586104 1.6606512  0.2663710  0.1228792
house_cat10+   1.17003712 0.67405594 2.0309692  0.1570355  0.2813713
                   z value     Pr(>|z|)
(Intercept)    -15.8772359 9.110557e-57
agegroups0-10    4.9442116 7.645264e-07
agegroups11-20   5.3313714 9.747386e-08
gendermale       3.0267824 2.471718e-03
house_cat6-10    2.1677478 3.017788e-02
house_cat10+     0.5581076 5.767709e-01

So you could also probably try this in your case.

Contrary to confint.defaut which is a directly readable R function, confint is a S3 dispatch method (thanks @Ben Bolker for the internal references in comments), and I didn't yet investigate further what could explain this surprising behaviour.

Another option seems to save scabiesrisk_summary in another variable.
I tried hard but was never able to reproduce the problem after doing so :

```{r}
scabiesrisk_OR <- exp(cbind(OR= coef(scabiesrisk), confint((scabiesrisk))))
scabiesrisk_summary <- summary(scabiesrisk)
scabiesrisk_final <- cbind(scabiesrisk_OR, scabiesrisk_summary$coefficients)
scabiesrisk_final
```

I strongly suspect that you forgot to include library(tidyverse) in your script. If tidyverse is loaded, then your code works fine. If it's not:

• the step where you try to mutate() (and use %>%) fails, so the scabies variable is never created within the scabies data set
• glm(scabies ~ ...) then interprets the response variable scabies as being the whole data set, and complains that the response variable is "invalid type(list)".

For this reason it's good practice to avoid having variables within data frames that have the same name as the data frames themselves ...

Your data transformation steps can be cleaned up a little bit (as.factor() is redundant; you can do all of the transformations as steps within a single mutate() call; as.numeric(x=="yes") is a shorter way to turn a string into a 0/1 variable ...) If I were going to do a lot more of this I would write a custom mycut() function that took breakpoints and a desired reference level as input arguments, constructed custom labels, and did the releveling.

library(tidyverse)
scabies <- (read.csv(file = "S1-Dataset_CSV.csv") %>%
            mutate(agegroups <- cut(age, c(0,10,20,Inf),
                                    labels = c("0-10","11-20","21+"),
                                    include.lowest = TRUE),
                   agegroups = relevel(agegroups, ref = "21+"),
                   house_cat = cut(house_inhabitants, c(0,5,10,Inf),
                                   labels = c("0-5","6-10","10+"),
                                   include.lowest = TRUE),
                   house_cat = relevel(house_cat, ref = "0-5"),
                   scabies = as.numeric(scabies_infestation=="yes"),
                   impetigo = as.numeric(impetigo_active=="yes"))
)