n <- 1000
x <- sort(rnorm(n))
alpha <- -1
beta <- 0.3
obs_sd <- 1.2
y <- rnorm(n, alpha + beta * x, obs_sd)
library(greta)
#>
#> Attaching package: 'greta'
#> The following objects are masked from 'package:stats':
#>
#> binomial, cov2cor, poisson
#> The following objects are masked from 'package:base':
#>
#> %*%, apply, backsolve, beta, chol2inv, colMeans, colSums, diag,
#> eigen, forwardsolve, gamma, identity, rowMeans, rowSums, sweep,
#> tapply
y_ga <- as_data(y)
#> ℹ Initialising python and checking dependencies, this may take a moment.
#> ✔ Initialising python and checking dependencies ... done!
#>
alpha <- normal(0, 10)
beta <- normal(0, 10)
obs_sd <- normal(0, 1, truncation = c(0, Inf))
eta <- alpha + beta * x
distribution(y_ga) <- normal(eta, obs_sd)
m <- model(alpha, beta, obs_sd)
draws <- mcmc(m)
#> running 4 chains simultaneously on up to 8 CPU cores
#>
#> warmup 0/1000 | eta: ?s warmup == 50/1000 | eta: 19s | 9% bad warmup ==== 100/1000 | eta: 10s | 4% bad warmup ====== 150/1000 | eta: 7s | 3% bad warmup ======== 200/1000 | eta: 6s | 2% bad warmup ========== 250/1000 | eta: 5s | 2% bad warmup =========== 300/1000 | eta: 4s | 2% bad warmup ============= 350/1000 | eta: 4s | 1% bad warmup =============== 400/1000 | eta: 3s | 1% bad warmup ================= 450/1000 | eta: 3s | 1% bad warmup =================== 500/1000 | eta: 2s | <1% bad warmup ===================== 550/1000 | eta: 2s | <1% bad warmup ======================= 600/1000 | eta: 2s | <1% bad warmup ========================= 650/1000 | eta: 2s | <1% bad warmup =========================== 700/1000 | eta: 1s | <1% bad warmup ============================ 750/1000 | eta: 1s | <1% bad warmup ============================== 800/1000 | eta: 1s | <1% bad warmup ================================ 850/1000 | eta: 1s | <1% bad warmup ================================== 900/1000 | eta: 0s | <1% bad warmup ==================================== 950/1000 | eta: 0s | <1% bad warmup ====================================== 1000/1000 | eta: 0s | <1% bad
#> sampling 0/1000 | eta: ?s sampling == 50/1000 | eta: 3s sampling ==== 100/1000 | eta: 3s sampling ====== 150/1000 | eta: 3s sampling ======== 200/1000 | eta: 2s sampling ========== 250/1000 | eta: 2s sampling =========== 300/1000 | eta: 2s sampling ============= 350/1000 | eta: 2s sampling =============== 400/1000 | eta: 2s sampling ================= 450/1000 | eta: 2s sampling =================== 500/1000 | eta: 1s sampling ===================== 550/1000 | eta: 1s sampling ======================= 600/1000 | eta: 1s sampling ========================= 650/1000 | eta: 1s sampling =========================== 700/1000 | eta: 1s sampling ============================ 750/1000 | eta: 1s sampling ============================== 800/1000 | eta: 1s sampling ================================ 850/1000 | eta: 0s sampling ================================== 900/1000 | eta: 0s sampling ==================================== 950/1000 | eta: 0s sampling ====================================== 1000/1000 | eta: 0s
# not this
eta_sim <- calculate(eta, values = draws, nsim = 100)
# this
y_sim <- calculate(y_ga, values = draws, nsim = 100)
library(DHARMa)
#> This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
resids <- createDHARMa(simulatedResponse = t(y_sim$y_ga[, , 1]),
observedResponse = y)
#> No fitted predicted response provided, using the mean of the simulations
plot(resids)
# can also convert to expected normality
res <- resids$scaledResiduals
hist(res, breaks = 100)
hist(runif(n), breaks = 100)
res_normal <- qnorm(res)
hist(res_normal)
hist(qnorm(runif(1e5), 3, 2))
hist(rnorm(1e5, 3, 2))
Created on 2025-02-27 with reprex v2.1.1
Session info
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.4.2 (2024-10-31)
#> os macOS Sequoia 15.1
#> system aarch64, darwin20
#> ui X11
#> language (EN)
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Australia/Hobart
#> date 2025-02-27
#> pandoc 3.2.1 @ /opt/homebrew/bin/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> abind 1.4-8 2024-09-12 [1] CRAN (R 4.4.1)
#> backports 1.5.0 2024-05-23 [1] CRAN (R 4.4.0)
#> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.4.0)
#> boot 1.3-31 2024-08-28 [2] CRAN (R 4.4.2)
#> callr 3.7.6 2024-03-25 [1] CRAN (R 4.4.0)
#> cli 3.6.3 2024-06-21 [1] CRAN (R 4.4.0)
#> coda 0.19-4.1 2024-01-31 [1] CRAN (R 4.4.0)
#> codetools 0.2-20 2024-03-31 [2] CRAN (R 4.4.2)
#> crayon 1.5.3 2024-06-20 [1] CRAN (R 4.4.0)
#> curl 6.2.0 2025-01-23 [1] CRAN (R 4.4.1)
#> DHARMa * 0.4.7 2024-10-18 [1] CRAN (R 4.4.1)
#> digest 0.6.37 2024-08-19 [1] CRAN (R 4.4.1)
#> doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.4.0)
#> dplyr 1.1.4 2023-11-17 [1] CRAN (R 4.4.0)
#> evaluate 1.0.1 2024-10-10 [1] CRAN (R 4.4.1)
#> fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.4.0)
#> foreach 1.5.2 2022-02-02 [1] CRAN (R 4.4.0)
#> fs 1.6.5 2024-10-30 [1] CRAN (R 4.4.1)
#> future 1.34.0 2024-07-29 [1] CRAN (R 4.4.0)
#> gap 1.6 2024-08-27 [1] CRAN (R 4.4.1)
#> gap.datasets 0.0.6 2023-08-25 [1] CRAN (R 4.4.1)
#> generics 0.1.3 2022-07-05 [1] CRAN (R 4.4.0)
#> globals 0.16.3 2024-03-08 [1] CRAN (R 4.4.0)
#> glue 1.8.0 2024-09-30 [1] CRAN (R 4.4.1)
#> greta * 0.5.0.9000 2025-01-24 [1] local
#> hms 1.1.3 2023-03-21 [1] CRAN (R 4.4.0)
#> htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.4.0)
#> httpuv 1.6.15 2024-03-26 [1] CRAN (R 4.4.0)
#> iterators 1.0.14 2022-02-05 [1] CRAN (R 4.4.0)
#> jsonlite 1.8.9 2024-09-20 [1] CRAN (R 4.4.1)
#> knitr 1.49 2024-11-08 [1] CRAN (R 4.4.1)
#> later 1.4.1 2024-11-27 [1] CRAN (R 4.4.1)
#> lattice 0.22-6 2024-03-20 [2] CRAN (R 4.4.2)
#> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.4.0)
#> listenv 0.9.1 2024-01-29 [1] CRAN (R 4.4.0)
#> lme4 1.1-35.5 2024-07-03 [1] CRAN (R 4.4.0)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.4.0)
#> MASS 7.3-64 2025-01-04 [1] CRAN (R 4.4.1)
#> Matrix 1.7-1 2024-10-18 [2] CRAN (R 4.4.2)
#> mgcv 1.9-1 2023-12-21 [2] CRAN (R 4.4.2)
#> mime 0.12 2021-09-28 [1] CRAN (R 4.4.0)
#> minqa 1.2.8 2024-08-17 [1] CRAN (R 4.4.0)
#> nlme 3.1-166 2024-08-14 [2] CRAN (R 4.4.2)
#> nloptr 2.1.1 2024-06-25 [1] CRAN (R 4.4.0)
#> parallelly 1.41.0 2024-12-18 [1] CRAN (R 4.4.1)
#> pillar 1.10.1 2025-01-07 [1] CRAN (R 4.4.1)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.4.0)
#> plyr 1.8.9 2023-10-02 [1] CRAN (R 4.4.0)
#> png 0.1-8 2022-11-29 [1] CRAN (R 4.4.0)
#> prettyunits 1.2.0 2023-09-24 [1] CRAN (R 4.4.0)
#> processx 3.8.5 2025-01-08 [1] CRAN (R 4.4.1)
#> progress 1.2.3 2023-12-06 [1] CRAN (R 4.4.0)
#> promises 1.3.2 2024-11-28 [1] CRAN (R 4.4.1)
#> ps 1.8.1 2024-10-28 [1] CRAN (R 4.4.1)
#> qgam 1.3.4 2021-11-22 [1] CRAN (R 4.4.0)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.4.0)
#> rbibutils 2.3 2024-10-04 [1] CRAN (R 4.4.1)
#> Rcpp 1.0.13-1 2024-11-02 [1] CRAN (R 4.4.1)
#> Rdpack 2.6.2 2024-11-15 [1] CRAN (R 4.4.1)
#> reprex 2.1.1 2024-07-06 [1] CRAN (R 4.4.0)
#> reticulate 1.40.0 2024-11-15 [1] CRAN (R 4.4.1)
#> rlang 1.1.5 2025-01-17 [1] CRAN (R 4.4.1)
#> rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.4.1)
#> rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.4.1)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.4.0)
#> shiny 1.10.0 2024-12-14 [1] CRAN (R 4.4.1)
#> tensorflow 2.16.0 2024-04-15 [1] CRAN (R 4.4.0)
#> tfautograph 0.3.2 2021-09-17 [1] CRAN (R 4.4.0)
#> tfruns 1.5.3 2024-04-19 [1] CRAN (R 4.4.0)
#> tibble 3.2.1 2023-03-20 [1] CRAN (R 4.4.0)
#> tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.4.0)
#> vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.4.0)
#> whisker 0.4.1 2022-12-05 [1] CRAN (R 4.4.0)
#> withr 3.0.2 2024-10-28 [1] CRAN (R 4.4.1)
#> xfun 0.50.5 2025-01-15 [1] Github (yihui/xfun@116d689)
#> xml2 1.3.6 2023-12-04 [1] CRAN (R 4.4.0)
#> xtable 1.8-4 2019-04-21 [1] CRAN (R 4.4.0)
#> yaml 2.3.10 2024-07-26 [1] CRAN (R 4.4.0)
#>
#> [1] /Users/nick/Library/R/arm64/4.4/library
#> [2] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
#>
#> ─ Python configuration ───────────────────────────────────────────────────────
#> python: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/bin/python
#> libpython: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/libpython3.10.dylib
#> pythonhome: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2:/Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2
#> version: 3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:51:49) [Clang 16.0.6 ]
#> numpy: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/python3.10/site-packages/numpy
#> numpy_version: 1.26.4
#> tensorflow: /Users/nick/Library/r-miniconda-arm64/envs/greta-env-tf2/lib/python3.10/site-packages/tensorflow
#>
#> NOTE: Python version was forced by use_python() function
#>
#> ──────────────────────────────────────────────────────────────────────────────