Why is the quantile function not working for this dplyr function?

Question

I'm working through Faraway's 2016 book Extending the Linear Model with R and have encountered an issue with the code that I don't know how to fix. Here is the relevant syntax leading up to the error:

#### Load Data & Libraries ####
library(faraway)
library(tidyverse)
data(wcgs)

#### Add Variables ####
wcgs$y <- ifelse(wcgs$chd == "no",0,1) # create binary response from chd 
wcgs$bmi <- with(wcgs,
                 703*wcgs$weight/(wcgs$height^2)) # create BMI variable

#### Create GLM Model ####
lmod <- glm(chd ~ height + cigs,
            family = binomial, 
            wcgs)

#### Mutate Data ####
wcgs <- mutate(wcgs,
               residuals=residuals(lmod), 
               linpred=predict(lmod)) # create residuals/pred values

And this is the part where the error arises (the third line which includes a mutate function:

#### Error Code (Last Line) ####
wcgsm <- na.omit(wcgs) # omit NA values
wcgsm <- mutate(wcgsm,
                predprob=predict(lmod,
                                 type="response")) # make pred data
gdf <- group_by(wcgsm, 
                cut(linpred,
                    breaks=unique(quantile(linpred, 
                                           (1:100)/101)))) # bin NA

Which gives me this error:

Error in `group_by()`:
! Problem adding computed columns.
Caused by error in `mutate()`:
! Problem while computing `..1 = cut(linpred, breaks = unique(quantile(linpred,
  (1:100)/101)))`.
✖ `..1` must be size 3140 or 1, not 3154.

I dont understand what this error means. When I run dim(wcgs), I get there are 3154 rows, and when I run dim(na.omit(wcgs)) I get 3140 rows. The only thing I can think of is that the predicted model values dont line up with the new na.omit data, but I'm not sure now how to work around that given the rest of this chapter uses this data manipulation.

Why are you blaming `quantile`? The error message says the problem is with `group_by`. — IRTFM, Jul 16 '22 at 04:54
I can see that now, but it doesn't seem to look much different from a previous command they made earlier in the book, which looks like this: `gdf <- group_by(wcgs, cut(linpred, breaks=unique(quantile(linpred, (1:100)/101))))`. The only difference seems to be the NA values being removed. — Shawn Hemelstrand, Jul 16 '22 at 05:14
Any chance you've created a variable called `linpred` at some point? — Hobo, Jul 16 '22 at 05:47
I'm getting a similar error in the `predict` line. It was solved by making the `newdata` argument explicit since it defaults to the data set the model was fitted to. Try `predprob = predict(lmod, newdata = wcgsm, type="response")` and see if the next line still throws an error. — Rui Barradas, Jul 16 '22 at 05:48
I just tried your `newdata=wcgsm` part and it still gave me the same error. — Shawn Hemelstrand, Jul 16 '22 at 05:55
Can you restart from a fresh R session and see if my suggestion solves the problem? I'm still getting the error in the prediction step. I'll post as an answer the full code I'm running. — Rui Barradas, Jul 16 '22 at 05:58
I just restarted it and now it suddenly works. But please post your answer as well so I can accept it given it fixed my issue. — Shawn Hemelstrand, Jul 16 '22 at 06:01
@ShawnHemelstrand sorry - I meant a separate variable (outside the mutate). But that obviously wasn't the issue - it was just a thought. Glad you've found a solution — Hobo, Jul 17 '22 at 00:44

score 1 · Accepted Answer · answered Jul 16 '22 at 06:13

predict methods for R's modeling functions always predict from the original data set the models were fitted to. To have a new data set, in this case a subset of the data wcgs, argument newdata must be explicitly set.
The error in the predict line at the bottom is therefore expected behavior.

#### Load Data & Libraries ####
suppressPackageStartupMessages({
  #library(faraway)
  library(dplyr)
})
data(wcgs, package = "faraway")

#### Add Variables ####
wcgs$y <- as.integer(wcgs$chd == "yes") # create binary response from chd 
wcgs$bmi <- with(wcgs, 703*weight/(height^2)) # create BMI variable

#### Create GLM Model ####
lmod <- glm(chd ~ height + cigs, family = binomial, data = wcgs)

#### Mutate Data ####
# create residuals/pred values
wcgs <- mutate(wcgs,
               residuals = residuals(lmod), 
               linpred = predict(lmod)) 

wcgsm <- na.omit(wcgs) # omit NA values
wcgsm <- mutate(wcgsm,
                predprob = predict(lmod, type="response")) # make pred data
#> Error in `mutate()`:
#> ! Problem while computing `predprob = predict(lmod, type = "response")`.
#> ✖ `predprob` must be size 3140 or 1, not 3154.

^{Created on 2022-07-16 by the reprex package (v2.0.1)}

See where the error comes from.

predprob_all <- predict(lmod, type = "response")
predprob_na.omit <- predict(lmod, newdata = wcgsm, type = "response")

length(predprob_all)
#> [1] 3154
length(predprob_na.omit)
#> [1] 3140

^{Created on 2022-07-16 by the reprex package (v2.0.1)}

These lengths are the values in the error message, once again, as expected.

There is also the problem of the quantiles in cut(., breaks) not spanning the entire range of linpred. Values outside the quantiles' range will become NA. This is solved with the two endpoints of the breaks vector.
And I have given a name to the binned vector.
The following code works and, I believe, does what is needed.

wcgsm <- na.omit(wcgs) # omit NA values
wcgsm <- mutate(wcgsm,
                predprob = predict(lmod, newdata = wcgsm, type="response")) # make pred data

breaks <- c(-Inf, 
            unique(quantile(wcgsm$linpred, (1:100)/101)),
             Inf)
gdf <- group_by(wcgsm, 
                bins = cut(linpred, breaks = breaks)) # bin NA

anyNA(gdf$bins)
#> [1] FALSE

^{Created on 2022-07-16 by the reprex package (v2.0.1)}

Thank you for the detailed explanation and your hard work. That did the trick. — Shawn Hemelstrand, Jul 16 '22 at 06:18

Why is the quantile function not working for this dplyr function?

1 Answers1