# standard figure size and generate clean outputknitr::opts_chunk$set(fig.height=4, fig.width=6, warning=FALSE, message=FALSE, cache=TRUE, error=T, echo=T)library(here)
here() starts at C:/Users/dan.killian/Documents/StatisticalRethinking
source(here("scripts/prep.R"))
Warning: package 'arm' was built under R version 4.1.3
Loading required package: MASS
Warning: package 'MASS' was built under R version 4.1.3
Loading required package: Matrix
Warning: package 'Matrix' was built under R version 4.1.3
Loading required package: lme4
Warning: package 'lme4' was built under R version 4.1.3
arm (Version 1.13-1, built: 2022-8-25)
Working directory is C:/Users/dan.killian/Documents/StatisticalRethinking/scripts
Warning: package 'BMA' was built under R version 4.1.3
Loading required package: survival
Loading required package: leaps
Loading required package: robustbase
Warning: package 'robustbase' was built under R version 4.1.2
Attaching package: 'robustbase'
The following object is masked from 'package:survival':
heart
Loading required package: inline
Loading required package: rrcov
Warning: package 'rrcov' was built under R version 4.1.2
Scalable Robust Estimators with High Breakdown Point (version 1.6-2)
Attaching package: 'rrcov'
The following objects are masked from 'package:lme4':
getData, isSingular
Warning: package 'brms' was built under R version 4.1.3
Loading required package: Rcpp
Warning: package 'Rcpp' was built under R version 4.1.3
Attaching package: 'Rcpp'
The following object is masked from 'package:inline':
registerPlugin
Loading 'brms' package (version 2.18.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: 'brms'
The following object is masked from 'package:robustbase':
epilepsy
The following object is masked from 'package:survival':
kidney
The following object is masked from 'package:lme4':
ngrps
The following object is masked from 'package:stats':
ar
Warning: package 'corrplot' was built under R version 4.1.2
corrplot 0.92 loaded
Attaching package: 'corrplot'
The following object is masked from 'package:arm':
corrplot
dummies-1.5.6 provided by Decision Patterns
Attaching package: 'dummies'
The following object is masked from 'package:lme4':
dummy
Warning: package 'DescTools' was built under R version 4.1.3
Attaching package: 'DescTools'
The following object is masked from 'package:rrcov':
Cov
Warning: package 'estimatr' was built under R version 4.1.3
Warning: package 'extrafont' was built under R version 4.1.3
Registering fonts with R
Warning: package 'forestplot' was built under R version 4.1.3
Loading required package: grid
Loading required package: checkmate
Warning: package 'checkmate' was built under R version 4.1.3
Loading required package: abind
Attaching package: 'janitor'
The following objects are masked from 'package:stats':
chisq.test, fisher.test
Warning: package 'memisc' was built under R version 4.1.3
Loading required package: lattice
Warning: package 'lattice' was built under R version 4.1.2
Attaching package: 'memisc'
The following objects are masked from 'package:DescTools':
%nin%, Cor, Cov, Mean, Median, Range, Var
The following object is masked from 'package:rrcov':
Cov
The following object is masked from 'package:Matrix':
as.array
The following objects are masked from 'package:stats':
contr.sum, contr.treatment, contrasts
The following object is masked from 'package:base':
as.array
Loading required package: rstan
Warning: package 'rstan' was built under R version 4.1.3
Loading required package: StanHeaders
Loading required package: ggplot2
Warning: package 'ggplot2' was built under R version 4.1.3
Attaching package: 'ggplot2'
The following object is masked from 'package:memisc':
syms
rstan (Version 2.21.7, GitRev: 2e1f913d3ca3)
For execution on a local, multicore CPU with excess RAM we recommend calling
options(mc.cores = parallel::detectCores()).
To avoid recompilation of unchanged Stan programs, we recommend calling
rstan_options(auto_write = TRUE)
Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
Attaching package: 'rstan'
The following object is masked from 'package:arm':
traceplot
Loading required package: parallel
rethinking (Version 2.13)
Attaching package: 'rethinking'
The following objects are masked from 'package:brms':
LOO, stancode, WAIC
The following objects are masked from 'package:arm':
logit, sim, standardize
The following object is masked from 'package:stats':
rstudent
Warning: package 'tidyr' was built under R version 4.1.3
Attaching package: 'tidyr'
The following object is masked from 'package:rstan':
extract
The following object is masked from 'package:reshape2':
smiths
The following objects are masked from 'package:Matrix':
expand, pack, unpack
Warning: package 'broom' was built under R version 4.1.3
Warning: package 'haven' was built under R version 4.1.3
Warning: package 'HH' was built under R version 4.1.3
Loading required package: latticeExtra
Warning: package 'latticeExtra' was built under R version 4.1.3
Attaching package: 'latticeExtra'
The following object is masked from 'package:ggplot2':
layer
Loading required package: multcomp
Warning: package 'multcomp' was built under R version 4.1.3
Loading required package: mvtnorm
Loading required package: TH.data
Warning: package 'TH.data' was built under R version 4.1.2
Attaching package: 'TH.data'
The following object is masked from 'package:MASS':
geyser
Loading required package: gridExtra
Attaching package: 'HH'
The following object is masked from 'package:rethinking':
logit
The following object is masked from 'package:DescTools':
OddsRatio
The following object is masked from 'package:brms':
mmc
The following object is masked from 'package:arm':
logit
Warning: package 'Hmisc' was built under R version 4.1.3
Loading required package: Formula
Attaching package: 'Hmisc'
The following objects are masked from 'package:memisc':
%nin%, html, Mean
The following objects are masked from 'package:DescTools':
%nin%, Label, Mean, Quantile
The following objects are masked from 'package:base':
format.pval, units
Warning: package 'lubridate' was built under R version 4.1.2
Attaching package: 'lubridate'
The following object is masked from 'package:HH':
interval
The following objects are masked from 'package:memisc':
as.interval, is.interval
The following objects are masked from 'package:base':
date, intersect, setdiff, union
Warning: package 'knitr' was built under R version 4.1.3
Attaching package: 'margins'
The following object is masked from 'package:brms':
marginal_effects
Warning: package 'magrittr' was built under R version 4.1.3
Attaching package: 'magrittr'
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:rstan':
extract
Attaching package: 'plotrix'
The following object is masked from 'package:arm':
rescale
Warning: package 'scales' was built under R version 4.1.3
Attaching package: 'scales'
The following object is masked from 'package:plotrix':
rescale
The following object is masked from 'package:memisc':
percent
The following object is masked from 'package:arm':
rescale
Attaching package: 'survey'
The following object is masked from 'package:Hmisc':
deff
The following object is masked from 'package:graphics':
dotchart
Warning: package 'srvyr' was built under R version 4.1.2
Attaching package: 'srvyr'
The following object is masked from 'package:Hmisc':
summarize
The following objects are masked from 'package:memisc':
collect, rename, syms
The following object is masked from 'package:MASS':
select
The following object is masked from 'package:stats':
filter
Warning: package 'foreign' was built under R version 4.1.3
Warning: package 'car' was built under R version 4.1.3
Loading required package: carData
Warning: package 'carData' was built under R version 4.1.2
Attaching package: 'car'
The following objects are masked from 'package:HH':
logit, vif
The following object is masked from 'package:rethinking':
logit
The following object is masked from 'package:memisc':
recode
The following object is masked from 'package:DescTools':
Recode
The following object is masked from 'package:arm':
logit
Warning: package 'ICC' was built under R version 4.1.3
Warning: package 'openxlsx' was built under R version 4.1.3
Warning: package 'readr' was built under R version 4.1.3
Attaching package: 'readr'
The following object is masked from 'package:scales':
col_factor
Warning: package 'readxl' was built under R version 4.1.3
Warning: package 'labelled' was built under R version 4.1.3
Warning: package 'texreg' was built under R version 4.1.2
Version: 1.38.5
Date: 2022-03-03
Author: Philip Leifeld (University of Essex)
Consider submitting praise using the praise or praise_interactive functions.
Please cite the JSS article in your publications -- see citation("texreg").
Attaching package: 'texreg'
The following object is masked from 'package:magrittr':
extract
The following object is masked from 'package:tidyr':
extract
The following object is masked from 'package:rstan':
extract
Warning: package 'sjmisc' was built under R version 4.1.2
Attaching package: 'sjmisc'
The following objects are masked from 'package:labelled':
to_character, to_factor
The following object is masked from 'package:Hmisc':
%nin%
The following object is masked from 'package:tidyr':
replace_na
The following object is masked from 'package:memisc':
%nin%
The following objects are masked from 'package:janitor':
remove_empty_cols, remove_empty_rows
The following object is masked from 'package:DescTools':
%nin%
Warning: package 'psych' was built under R version 4.1.3
Attaching package: 'psych'
The following object is masked from 'package:car':
logit
The following objects are masked from 'package:scales':
alpha, rescale
The following object is masked from 'package:plotrix':
rescale
The following object is masked from 'package:Hmisc':
describe
The following object is masked from 'package:HH':
logit
The following objects are masked from 'package:rethinking':
logistic, logit, sim
The following object is masked from 'package:rstan':
lookup
The following objects are masked from 'package:ggplot2':
%+%, alpha
The following objects are masked from 'package:DescTools':
AUC, ICC, SD
The following object is masked from 'package:brms':
cs
The following objects are masked from 'package:arm':
logit, rescale, sim
Warning: package 'dplyr' was built under R version 4.1.3
Attaching package: 'dplyr'
The following object is masked from 'package:car':
recode
The following objects are masked from 'package:Hmisc':
src, summarize
The following object is masked from 'package:gridExtra':
combine
The following objects are masked from 'package:memisc':
collect, recode, rename, syms
The following object is masked from 'package:MASS':
select
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Warning: package 'tidybayes' was built under R version 4.1.2
Attaching package: 'tidybayes'
The following object is masked from 'package:DescTools':
Mode
The following objects are masked from 'package:brms':
dstudent_t, pstudent_t, qstudent_t, rstudent_t
Warning: package 'viridis' was built under R version 4.1.2
Loading required package: viridisLite
Warning: package 'viridisLite' was built under R version 4.1.3
Attaching package: 'viridis'
The following object is masked from 'package:scales':
viridis_pal
Warning: package 'ggridges' was built under R version 4.1.3
Attaching package: 'ggridges'
The following objects are masked from 'package:tidybayes':
scale_point_color_continuous, scale_point_color_discrete,
scale_point_colour_continuous, scale_point_colour_discrete,
scale_point_fill_continuous, scale_point_fill_discrete,
scale_point_size_continuous
Warning: package 'DT' was built under R version 4.1.3
Warning: package 'workflowr' was built under R version 4.1.2
This is workflowr version 1.7.0
Run ?workflowr for help getting started
Warning: package 'pointblank' was built under R version 4.1.3
Warning: package 'cmdstanr' was built under R version 4.1.3
This is cmdstanr version 0.5.3
- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
A newer version of CmdStan is available. See ?install_cmdstan() to install it.
To disable this check set option or environment variable CMDSTANR_NO_VER_CHECK=TRUE.
Warning: package 'patchwork' was built under R version 4.1.3
Attaching package: 'patchwork'
The following object is masked from 'package:MASS':
area
Warning: package 'tidyverse' was built under R version 4.1.3
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
Found more than one class "atomicVector" in cache; using the first, from namespace 'Matrix'
Also defined by 'Rmpfr'
-- Attaching packages --------------------------------------- tidyverse 1.3.2 --
v tibble 3.1.8 v stringr 1.4.1
v purrr 0.3.5 v forcats 0.5.2
Warning: package 'tibble' was built under R version 4.1.3
Warning: package 'purrr' was built under R version 4.1.3
Warning: package 'stringr' was built under R version 4.1.3
Warning: package 'forcats' was built under R version 4.1.3
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x purrr::%@%() masks memisc::%@%()
x psych::%+%() masks ggplot2::%+%()
x tibble::add_case() masks sjmisc::add_case()
x psych::alpha() masks scales::alpha(), ggplot2::alpha()
x lubridate::as.difftime() masks base::as.difftime()
x lubridate::as.interval() masks memisc::as.interval()
x readr::col_factor() masks scales::col_factor()
x dplyr::collect() masks srvyr::collect(), memisc::collect()
x dplyr::combine() masks gridExtra::combine()
x lubridate::date() masks base::date()
x purrr::discard() masks scales::discard()
x tidyr::expand() masks Matrix::expand()
x texreg::extract() masks magrittr::extract(), tidyr::extract(), rstan::extract()
x dplyr::filter() masks srvyr::filter(), stats::filter()
x lubridate::intersect() masks base::intersect()
x lubridate::interval() masks HH::interval()
x lubridate::is.interval() masks memisc::is.interval()
x purrr::is_empty() masks sjmisc::is_empty()
x dplyr::lag() masks stats::lag()
x latticeExtra::layer() masks ggplot2::layer()
x purrr::map() masks rethinking::map()
x tidyr::pack() masks Matrix::pack()
x dplyr::recode() masks car::recode(), memisc::recode()
x dplyr::rename() masks srvyr::rename(), memisc::rename()
x sjmisc::replace_na() masks tidyr::replace_na()
x dplyr::select() masks srvyr::select(), MASS::select()
x purrr::set_names() masks magrittr::set_names()
x lubridate::setdiff() masks base::setdiff()
x purrr::some() masks car::some()
x dplyr::src() masks Hmisc::src()
x dplyr::summarize() masks srvyr::summarize(), Hmisc::summarize()
x dplyr::syms() masks srvyr::syms(), ggplot2::syms(), memisc::syms()
x purrr::transpose() masks HH::transpose()
x lubridate::union() masks base::union()
x tidyr::unpack() masks Matrix::unpack()
x tibble::view() masks memisc::view()
Conditional Manatees
Building an interaction
Let’s look at the relationship between terrain ruggedness and economic performance. Generally speaking, one might think that the less rugged the terrain, the easier it is to conduct economic activity and therefore the greater economic performance of a country.
To examine this relationship, we’ll look at terrain ruggedness, measured as the meters of elevation difference per grid point, and the log of country GDP in the year 2000.
## R code 8.1 ## prepare datadata(rugged)d <- rugged %>%filter(complete.cases(rgdppc_2000))rm(rugged)d <- d %>%mutate(log_gdp =log(d$rgdppc_2000), # make log version of outcomelog_gdp_std = log_gdp /mean(log_gdp), # text (proportion of mean)log_gdp_sc =scale(log_gdp), # mean 0 sd 1rugged_std = rugged /max(rugged), # text version (0-1)rugged_sc =scale(rugged)) # mean 0 sd 1str(d)
'data.frame': 170 obs. of 56 variables:
$ isocode : Factor w/ 234 levels "ABW","AFG","AGO",..: 3 5 8 9 10 12 13 14 15 16 ...
$ isonum : int 24 8 784 32 51 28 36 40 31 108 ...
$ country : Factor w/ 234 levels "Afghanistan",..: 6 2 219 9 10 8 12 13 14 34 ...
$ rugged : num 0.858 3.427 0.769 0.775 2.688 ...
$ rugged_popw : num 0.714 1.597 0.316 0.22 0.934 ...
$ rugged_slope : num 2.27 10.45 2.11 2.27 8.18 ...
$ rugged_lsd : num 0.228 1.006 0.191 0.226 0.799 ...
$ rugged_pc : num 4.91 62.13 6.14 9.41 50.56 ...
$ land_area : num 124670 2740 8360 273669 2820 ...
$ lat : num -12.3 41.1 23.9 -35.4 40.3 ...
$ lon : num 17.6 20.1 54.3 -65.2 44.9 ...
$ soil : num 26.7 68.1 0 35.7 30.1 ...
$ desert : num 0.425 0 77.28 0 0 ...
$ tropical : num 44.3 0 0 0 0 ...
$ dist_coast : num 0.428 0.048 0.065 0.352 0.348 0.001 0.336 0.242 0.584 0.906 ...
$ near_coast : num 13.2 94.7 75.7 13 0 ...
$ gemstones : int 47756 0 0 0 0 0 264154 0 0 0 ...
$ rgdppc_2000 : num 1795 3703 20604 12174 2422 ...
$ rgdppc_1950_m : num 1052 1001 15798 4987 NA ...
$ rgdppc_1975_m : num 1073 2289 25465 8122 NA ...
$ rgdppc_2000_m : num 765 2741 17568 8544 4565 ...
$ rgdppc_1950_2000_m : num 1107 1932 20120 6927 NA ...
$ q_rule_law : num -1.567 -0.82 0.913 0.033 -0.453 ...
$ cont_africa : int 1 0 0 0 0 0 0 0 0 1 ...
$ cont_asia : int 0 0 1 0 1 0 0 0 1 0 ...
$ cont_europe : int 0 1 0 0 0 0 0 1 0 0 ...
$ cont_oceania : int 0 0 0 0 0 0 1 0 0 0 ...
$ cont_north_america : int 0 0 0 0 0 1 0 0 0 0 ...
$ cont_south_america : int 0 0 0 1 0 0 0 0 0 0 ...
$ legor_gbr : int 0 0 1 0 0 1 1 0 0 0 ...
$ legor_fra : int 1 0 0 1 0 0 0 0 0 1 ...
$ legor_soc : int 0 1 0 0 1 0 0 0 1 0 ...
$ legor_deu : int 0 0 0 0 0 0 0 1 0 0 ...
$ legor_sca : int 0 0 0 0 0 0 0 0 0 0 ...
$ colony_esp : int 0 0 0 1 0 0 0 0 0 0 ...
$ colony_gbr : int 0 0 1 0 0 1 1 0 0 0 ...
$ colony_fra : int 0 0 0 0 0 0 0 0 0 0 ...
$ colony_prt : int 1 0 0 0 0 0 0 0 0 0 ...
$ colony_oeu : int 0 0 0 0 0 0 0 0 0 1 ...
$ africa_region_n : int 0 0 0 0 0 0 0 0 0 0 ...
$ africa_region_s : int 0 0 0 0 0 0 0 0 0 0 ...
$ africa_region_w : int 0 0 0 0 0 0 0 0 0 0 ...
$ africa_region_e : int 0 0 0 0 0 0 0 0 0 1 ...
$ africa_region_c : int 1 0 0 0 0 0 0 0 0 0 ...
$ slave_exports : num 3610000 0 0 0 0 ...
$ dist_slavemkt_atlantic: num 5.67 NA NA NA NA ...
$ dist_slavemkt_indian : num 6.98 NA NA NA NA ...
$ dist_slavemkt_saharan : num 4.93 NA NA NA NA ...
$ dist_slavemkt_redsea : num 3.87 NA NA NA NA ...
$ pop_1400 : int 1223208 200000 19200 276632 105743 747 200000 1250000 226250 239693 ...
$ european_descent : num 2 100 0 89.9 0.5 ...
$ log_gdp : num 7.49 8.22 9.93 9.41 7.79 ...
$ log_gdp_std : num 0.88 0.965 1.166 1.104 0.915 ...
$ log_gdp_sc : num [1:170, 1] -0.878 -0.257 1.214 0.763 -0.621 ...
..- attr(*, "scaled:center")= num 8.52
..- attr(*, "scaled:scale")= num 1.17
$ rugged_std : num 0.138 0.553 0.124 0.125 0.433 ...
$ rugged_sc : num [1:170, 1] -0.407 1.792 -0.483 -0.478 1.159 ...
..- attr(*, "scaled:center")= num 1.33
..- attr(*, "scaled:scale")= num 1.17
The hypothesized negative relationship between ruggedness and economic performance holds for all non-Africa countries. In Africa, there is a positive relationship between ruggedness and economic performance. How could this be?
## R code 8.2m8.1<-quap(alist( log_gdp_std ~dnorm( mu , sigma ) , mu <- a + b*( rugged_std -0.215 ) , # deviation from mean a ~dnorm( 1 , 1 ) , b ~dnorm( 0 , 1 ) , sigma ~dexp( 1 ) ) , data=dd )
Error in quap(alist(log_gdp_std ~ dnorm(mu, sigma), mu <- a + b * (rugged_std - : object 'dd' not found
Let’s look at the priors we created for this model.
## R code 8.3set.seed(7)prior <-extract.prior( m8.1 )
Error in h(simpleError(msg, call)): error in evaluating the argument 'fit' in selecting a method for function 'extract.prior': object 'm8.1' not found
plot(density(prior$a))
Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object of type 'closure' is not subsettable
Error in quap(alist(log_gdp_std ~ dnorm(mu, sigma), mu <- a[cid] + b * : object 'dd' not found
## R code 8.9compare( m8.1 , m8.2 )
Error in compare(m8.1, m8.2): object 'm8.1' not found
## R code 8.10precis( m8.2 , depth=2 )
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'precis': object 'm8.2' not found
## R code 8.11post <-extract.samples(m8.2)
Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'extract.samples': object 'm8.2' not found
diff_a1_a2 <- post$a[,1] - post$a[,2]
Error in eval(expr, envir, enclos): object 'post' not found
PI( diff_a1_a2 )
Error in quantile(samples, probs = c(a, 1 - a)): object 'diff_a1_a2' not found
## R code 8.12rugged.seq <-seq( from=-0.1 , to=1.1 , length.out=30 )# compute mu over samples, fixing cid=2 and then cid=1mu.NotAfrica <-link( m8.2 ,data=data.frame( cid=2 , rugged_std=rugged.seq ) )
Error in h(simpleError(msg, call)): error in evaluating the argument 'fit' in selecting a method for function 'link': object 'm8.2' not found
Error in compare(m8.1, m8.2, m8.3, func = PSIS): object 'm8.1' not found
## R code 8.16plot( PSIS( m8.3 , pointwise=TRUE )$k )
Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': error in evaluating the argument 'object' in selecting a method for function 'PSIS': object 'm8.3' not found
## R code 8.17# plot Africa - cid=1d.A1 <- dd[ dd$cid==1 , ]
Error in eval(expr, envir, enclos): object 'dd' not found
plot( d.A1$rugged_std , d.A1$log_gdp_std , pch=16 , col=rangi2 ,xlab="ruggedness (standardized)" , ylab="log GDP (as proportion of mean)" ,xlim=c(0,1) )
Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'd.A1' not found
mu <-link( m8.3 , data=data.frame( cid=1 , rugged_std=rugged_seq ) )
Error in h(simpleError(msg, call)): error in evaluating the argument 'fit' in selecting a method for function 'link': object 'm8.3' not found
mu_mean <-apply( mu , 2 , mean )
Error in apply(mu, 2, mean): object 'mu' not found
mu_ci <-apply( mu , 2 , PI , prob=0.97 )
Error in apply(mu, 2, PI, prob = 0.97): object 'mu' not found
lines( rugged_seq , mu_mean , lwd=2 )
Error in xy.coords(x, y): object 'mu_mean' not found
