Search code examples
rinteractionquantilequartileggeffects

ggpredict and quartiles: predictions with only one data point?


I am using ggefects::ggpredict() to plot a regression model including two-way interactions: categorical:continuous, continuous:continuous, continuous:continuous (2nd degree polynomial). My first choice is to plot the interaction between continuous variables using quartiles, including minimum and maximum value, by setting terms = c("var.cont1[all]", "var.cont2[quart]"), as shown in this ggeffects vignette. In this vigentte, the minimum and maximum value of the variable of interest have many observations.

However, in my data I have a single observation for both mininimum and maximum, and ggpredict will plot predictions just fine, although with very wide confidence intervals for the max value. In such cases, will ggpredict() plot a prediction based on a single data point? I want to be sure I understand waht's going on. Many thanks in advance!

I have provided a reproducible example below, in which I force the data to have a single observation for the min and max values. The data comes from this ggeffects vignette.

library(ggeffects)
#> Warning: package 'ggeffects' was built under R version 4.2.3
library(tidyverse)
#> Warning: package 'ggplot2' was built under R version 4.2.2
#> Warning: package 'tibble' was built under R version 4.2.3
#> Warning: package 'tidyr' was built under R version 4.2.3
#> Warning: package 'purrr' was built under R version 4.2.3
#> Warning: package 'dplyr' was built under R version 4.2.3
#> Warning: package 'stringr' was built under R version 4.2.3

data(efc)

summary(efc$c12hour)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>    4.00   10.00   20.00   42.40   42.75  168.00       6

c(sum(!is.na(efc$c12hour) & efc$c12hour == 4), sum(!is.na(efc$c12hour) & efc$c12hour == 168))
#> [1] 24 86

efc2 <- efc[!is.na(efc$c12hour), ]
efc2 <- arrange(efc2, c12hour)
efc2$c12hour[1] <- 3
efc2$c12hour[nrow(efc2)] <- 169
summary(efc2$c12hour)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>    3.00   10.00   20.00   42.40   42.75  169.00
table(efc2$c12hour)
#> 
#>   3   4   5   6   7   8   9  10  11  12  14  15  16  17  18  20  21  22  24  25 
#>   1  23  29  60  32  42  11  70   6  23  34  43  12   4  16  58  12  11  20  28 
#>  26  27  28  30  35  36  39  40  42  43  45  48  49  50  55  56  59  60  62  65 
#>   1   2  28  40  24   1   1  35   9   1   4   4   2  22   1   6   1   7   1   3 
#>  70  77  80  84  85  89  90  91 100 105 110 118 119 120 125 126 128 130 140 148 
#>  13   3   5   6   2   1   2   3  13   1   4   1   1   8   1   1   1   1   6   1 
#> 150 154 160 161 162 168 169 
#>   4   1   6   2   1  85   1

fit <- lm(barthtot ~ c12hour * c161sex + neg_c_7, data = efc2)
plot(ggpredict(fit, terms = c("c161sex", "c12hour[quart]")))

sessionInfo()
#> R version 4.2.0 (2022-04-22 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 22000)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
#> [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.utf8    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] forcats_0.5.1     stringr_1.5.0     dplyr_1.1.2       purrr_1.0.1      
#>  [5] readr_2.1.2       tidyr_1.3.0       tibble_3.2.1      ggplot2_3.4.0    
#>  [9] tidyverse_1.3.1   ggeffects_1.3.1.5
#> 
#> loaded via a namespace (and not attached):
#>  [1] lubridate_1.8.0  assertthat_0.2.1 digest_0.6.29    utf8_1.2.2      
#>  [5] R6_2.5.1         cellranger_1.1.0 backports_1.4.1  reprex_2.0.2    
#>  [9] evaluate_0.15    httr_1.4.3       highr_0.9        pillar_1.9.0    
#> [13] rlang_1.1.1      readxl_1.4.0     rstudioapi_0.13  rmarkdown_2.14  
#> [17] labeling_0.4.2   munsell_0.5.0    broom_0.8.0      compiler_4.2.0  
#> [21] modelr_0.1.8     xfun_0.40        pkgconfig_2.0.3  htmltools_0.5.6 
#> [25] insight_0.19.5   tidyselect_1.2.0 fansi_1.0.3      crayon_1.5.1    
#> [29] tzdb_0.3.0       dbplyr_2.1.1     withr_2.5.0      grid_4.2.0      
#> [33] jsonlite_1.8.0   gtable_0.3.0     lifecycle_1.0.3  DBI_1.1.2       
#> [37] magrittr_2.0.3   scales_1.2.1     datawizard_0.7.0 cli_3.6.1       
#> [41] stringi_1.7.6    farver_2.1.0     fs_1.5.2         snakecase_0.11.0
#> [45] xml2_1.3.3       ellipsis_0.3.2   generics_0.1.2   vctrs_0.6.3     
#> [49] sjlabelled_1.2.0 tools_4.2.0      glue_1.6.2       hms_1.1.1       
#> [53] fastmap_1.1.0    yaml_2.3.5       colorspace_2.0-3 rvest_1.0.2     
#> [57] knitr_1.39       haven_2.5.0

Created on 2023-10-05 with reprex v2.0.2


Solution

  • Internally, ggpredict() call predict() (and ggemmeans() calls emmeans(), ggeffect() calls effect()). The predictions are based on the model's coefficients - the confidence intervals (resp. how narrow / wide these are) depends on the amount of data.

    Maybe this vignette clarifies how predicted values relate to regression coefficients.

    In your particular example, I assume that due to having only one data point at the "tails", confidence intervals are wide, indicating that based on the few data, there's a larger range of plausible values for your predictions (higher "uncertainty", if you like).