# Modelling sigma

Hi,

I would like to decompose sigma into sigma_method_old, sigma_method_new and sigma_leftover but the model cannot estimate the parameters of simulated data.

I would like to allow a variable with 2-levels to have each own variance term (something like unequal variances in anova). I’m looking to do it in rstan but here I’m using brms which is a bit simpler to me. According to bmrs and these examples I could use the sigma ~ formula but I’d like to keep,

1. the model error term and
2. the 2 variances.

I’ve come with this model

prior_bm <- c(
prior(normal(10, .1), class = Intercept),
prior(exponential(2), class = sd),
prior(exponential(2), class = sd,  dpar = sigma),
prior(normal(0, .5),  class = Intercept, dpar = sigma))

bm <- brm(
bf(value ~ 1 + (1|year) + (1|subject),
sigma ~ 1 + (0+method|method)),
prior = prior_bm,
data = d1, cores = 4, chains = 4, seed = 2802)


Which is giving what I was looking for

1. Total sigma (sigma_Intercept)
2. A variance estimation per level (sd(sigma_methodnew) & sd(sigma_methodold)

But since I’ve simulated the data I know that this estimation for the old method is double the real.
The simulated data have the following sd

• new method: 0.6,
• old method 0.3.
> exp(-.034) # sigma_Intercept - close to real
[1] 0.9665715
> exp(-.034- 0.36) # sigma_Intercept - sd(sigma_methodnew) to keep the error from new method
[1] 0.6743541
> exp(-.034 - 0.45) # sigma_Intercept - sd(sigma_methodold) - double the real
[1] 0.6163132


Is the model I’m using correct? Any ideas what migth be wrong?

 Family: gaussian
Links: mu = identity; sigma = log
Formula: value ~ 1 + (1 | year) + (1 | subject)
sigma ~ 1 + (0 + method | method)
Data: d1 (Number of observations: 424)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000

Group-Level Effects:
~subject (Number of levels: 212)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.09      0.06     0.00     0.23 1.01      656     1445

~year (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.05      0.04     0.00     0.14 1.00     1563     1358

~method (Number of levels: 2)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(sigma_methodnew)                      0.36      0.35     0.01     1.32 1.00     1538     1462
sd(sigma_methodold)                      0.45      0.37     0.02     1.46 1.00     1208     1241
cor(sigma_methodnew,sigma_methodold)    -0.01      0.57    -0.95     0.95 1.00     2429     1786

Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept           9.99      0.03     9.92    10.06 1.00     2284     2244
sigma_Intercept    -0.34      0.22    -0.75     0.15 1.00      972     1450

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).


Here is the data I’v simulated
The within sd in new is double compared to old.

d1<- structure(list(year = c("A", "A", "A", "A", "A", "A", "A", "A",
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A",
"A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "B", "B", "B",
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
"B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B",
"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
"C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C",
"C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D", "D", "D",
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D",
"D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D", "D",
"D", "D", "D", "D", "D", "D", "E", "E", "E", "E", "E", "E", "E",
"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E",
"E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "E", "F", "F",
"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F",
"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F",
"F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F", "F",
"F", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
"G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G", "G",
"G", "G", "G", "G", "H", "H", "H", "H", "H", "H", "H", "H", "H",
"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
"H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H", "H",
"H", "H", "H", "H", "H", "H", "H", "I", "I", "I", "I", "I", "I",
"I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I",
"I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I", "I",
"I", "I", "I", "I", "I", "I", "I", "I", "J", "J", "J", "J", "J",
"J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J",
"J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J",
"J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J", "J"
), subject = c(1L, 1L, 2L, 2L, 3L, 3L, 4L, 4L, 5L, 5L, 6L, 6L,
7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 12L, 12L, 13L, 13L,
14L, 14L, 15L, 15L, 16L, 16L, 17L, 17L, 18L, 18L, 19L, 19L, 20L,
20L, 21L, 21L, 22L, 22L, 45L, 45L, 46L, 46L, 47L, 47L, 48L, 48L,
49L, 49L, 50L, 50L, 51L, 51L, 52L, 52L, 53L, 53L, 54L, 54L, 55L,
55L, 56L, 56L, 57L, 57L, 58L, 58L, 59L, 59L, 60L, 60L, 61L, 61L,
62L, 62L, 63L, 63L, 64L, 64L, 65L, 65L, 87L, 87L, 88L, 88L, 89L,
89L, 90L, 90L, 91L, 91L, 92L, 92L, 93L, 93L, 94L, 94L, 95L, 95L,
96L, 96L, 97L, 97L, 98L, 98L, 99L, 99L, 100L, 100L, 101L, 101L,
102L, 102L, 103L, 103L, 104L, 104L, 105L, 105L, 106L, 106L, 107L,
107L, 108L, 108L, 131L, 131L, 132L, 132L, 133L, 133L, 134L, 134L,
135L, 135L, 136L, 136L, 137L, 137L, 138L, 138L, 139L, 139L, 140L,
140L, 141L, 141L, 142L, 142L, 143L, 143L, 144L, 144L, 145L, 145L,
146L, 146L, 147L, 147L, 148L, 148L, 149L, 149L, 150L, 150L, 171L,
171L, 172L, 172L, 173L, 173L, 174L, 174L, 175L, 175L, 176L, 176L,
177L, 177L, 178L, 178L, 179L, 179L, 180L, 180L, 181L, 181L, 182L,
182L, 183L, 183L, 184L, 184L, 185L, 185L, 186L, 186L, 187L, 187L,
188L, 188L, 189L, 189L, 190L, 190L, 191L, 191L, 192L, 192L, 215L,
215L, 216L, 216L, 217L, 217L, 218L, 218L, 219L, 219L, 220L, 220L,
221L, 221L, 222L, 222L, 223L, 223L, 224L, 224L, 225L, 225L, 226L,
226L, 227L, 227L, 228L, 228L, 229L, 229L, 230L, 230L, 231L, 231L,
232L, 232L, 233L, 233L, 234L, 234L, 235L, 235L, 257L, 257L, 258L,
258L, 259L, 259L, 260L, 260L, 261L, 261L, 262L, 262L, 263L, 263L,
264L, 264L, 265L, 265L, 266L, 266L, 267L, 267L, 268L, 268L, 269L,
269L, 270L, 270L, 271L, 271L, 272L, 272L, 273L, 273L, 274L, 274L,
275L, 275L, 276L, 276L, 277L, 277L, 299L, 299L, 300L, 300L, 301L,
301L, 302L, 302L, 303L, 303L, 304L, 304L, 305L, 305L, 306L, 306L,
307L, 307L, 308L, 308L, 309L, 309L, 310L, 310L, 311L, 311L, 312L,
312L, 313L, 313L, 314L, 314L, 315L, 315L, 316L, 316L, 317L, 317L,
318L, 318L, 319L, 319L, 341L, 341L, 342L, 342L, 343L, 343L, 344L,
344L, 345L, 345L, 346L, 346L, 347L, 347L, 348L, 348L, 349L, 349L,
350L, 350L, 351L, 351L, 352L, 352L, 353L, 353L, 354L, 354L, 355L,
355L, 356L, 356L, 357L, 357L, 358L, 358L, 359L, 359L, 360L, 360L,
381L, 381L, 382L, 382L, 383L, 383L, 384L, 384L, 385L, 385L, 386L,
386L, 387L, 387L, 388L, 388L, 389L, 389L, 390L, 390L, 391L, 391L,
392L, 392L, 393L, 393L, 394L, 394L, 395L, 395L, 396L, 396L, 397L,
397L, 398L, 398L, 399L, 399L, 400L, 400L, 401L, 401L, 402L, 402L
), method = c("old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new", "old", "new", "old", "new", "old", "new",
"old", "new", "old", "new", "old", "new", "old", "new", "old",
"new", "old", "new"), value = c(9.64666806231916, 11.1406419173765,
10.4255105848025, 11.2457821661948, 9.07356103467091, 9.45910433300655,
9.68186990403523, 10.2287331839032, 9.86930975551746, 10.3965714296506,
10.9113399620658, 8.24846280586141, 8.98225815774507, 11.1336520796119,
10.3070993770525, 8.87715191233007, 10.4557673582883, 9.86343249579437,
9.37038028432163, 9.96782883358531, 9.32228071391576, 8.15960239924106,
10.200148530534, 9.72383720912501, 11.0241788453915, 10.1807015126005,
10.4099130279561, 10.6248177280729, 9.32825545794063, 10.8145414023864,
9.4600725188663, 9.66927819917436, 11.0237248719275, 8.49430572189433,
10.0372192827134, 9.98781204165761, 9.63540324398594, 9.89709696692786,
10.420399542048, 8.77030250314109, 9.98244042373313, 9.4129421747976,
11.1516531044121, 11.1669941814043, 9.4198436434474, 10.2996453494226,
10.0575786643415, 10.2678425179664, 11.0715050307537, 10.955053881425,
10.2307753004977, 11.5240462794419, 8.81795652195423, 10.0067491796938,
9.63492072421524, 11.0050036031399, 9.34314797915736, 9.72442585814951,
10.3680122635163, 9.10346731525654, 9.84791036611555, 9.53221752827408,
9.21810362387959, 9.56636553222821, 10.1718365692315, 9.88400974088326,
8.97336435444177, 8.92291905088711, 10.4393708780108, 10.5516172114199,
10.222692730273, 10.0779855422247, 10.0979092385019, 9.9295180029111,
10.36451663369, 9.23265126579389, 8.99059350228469, 11.1283448645178,
10.1031624272184, 11.6512149781172, 9.65814308003735, 10.2768123521153,
10.1403423365961, 10.0017093202053, 9.33295084519743, 9.81862429720658,
10.1888378420736, 9.98633824809675, 10.0134483304201, 10.5505687917062,
9.83310695654405, 8.98543124563584, 9.27453009436458, 9.65924585256207,
9.30765342355441, 10.2171508316655, 9.59157524457701, 12.5263240009562,
9.73970722723746, 10.6099417224675, 9.41976450587833, 9.34554958035681,
10.8922029377097, 8.5619131726814, 10.1287061485434, 8.74918813209084,
10.8407248266056, 8.6767637625822, 10.4139972830072, 9.77209072884671,
9.3877371836339, 10.8259383887078, 9.5936064956609, 11.1581888822931,
10.2033795417513, 9.29865409831781, 9.48470597337375, 11.1973280371938,
9.84234399073567, 10.4182996372002, 10.2955906689641, 9.95320408037742,
10.1333719400826, 10.9229796391423, 10.1217875718174, 8.78671889189415,
10.5398452627003, 9.77380840868652, 9.85425878561668, 9.5531601336851,
9.67439520630761, 8.86947503420209, 9.88734689277898, 9.45095393829788,
9.37647883084962, 9.16911939624334, 9.56072034314916, 9.71126797566575,
9.75437147381779, 9.20271114799561, 10.6311987251933, 11.4680106702452,
9.63525804354694, 10.5332882036465, 9.87219391749674, 9.56933005069624,
9.72752465375751, 10.8878740847704, 10.1155128640606, 8.85576213831347,
9.66185749135493, 9.40362690246215, 10.9858567907881, 9.39746722480769,
9.71201627563607, 9.10564777922846, 8.85156844105756, 9.72419010004377,
9.73397872618711, 9.91489429221665, 9.556990940533, 11.5294059234401,
9.51656745317379, 9.7823717719027, 9.24831749797161, 11.1513977967016,
10.5330848851481, 9.71255512098891, 10.4418147403912, 11.6467420876022,
10.1680310997984, 9.80179368592425, 8.79523810131724, 9.6745850797733,
9.45608277288618, 10.7545639686315, 10.241305377294, 10.6087812808745,
11.4399003647411, 10.2296100794102, 10.5772408532628, 10.1015936844673,
10.0854824440185, 10.4256594425489, 9.81181066715515, 9.91628154544913,
10.5498161987714, 10.682072920475, 10.8312293701424, 10.5393773268653,
9.77765614044698, 9.67397098959081, 10.3030742758939, 9.45047408829272,
10.0801269528226, 9.58011786025431, 10.4331370404103, 9.20745935177094,
9.88551800404376, 9.68684680565471, 9.91319447213167, 9.75917319275408,
11.1238137707432, 10.556274664099, 10.3428054099433, 11.162674093232,
9.98380821078587, 10.2833087183452, 9.62851531531102, 8.42315754582768,
10.9237084040946, 9.60708060464494, 10.163533706097, 11.2762324363098,
10.6287372568558, 9.4691608645923, 10.5536660586234, 9.87029571605979,
9.91242282528335, 9.83820147675914, 9.93458656946169, 9.95033563568217,
9.75405723773715, 9.82487237460319, 10.1186320742338, 8.93629947443907,
10.8576624893521, 9.87202812935888, 10.3831902930352, 10.1920404066772,
9.5239514935371, 10.1698361124872, 9.63501931196244, 9.53434614358507,
10.9760422590876, 9.60394544440441, 10.0047619344575, 10.1474452498884,
10.2974314681013, 9.85179508843413, 9.84137108944241, 9.98013033579613,
11.0225296651751, 9.47406390418915, 10.8124614176175, 8.99243382684779,
10.0105822617517, 10.6768402756995, 9.8315583622288, 9.88555633739469,
9.25171311568785, 9.95376597949405, 9.97043183957485, 9.12008196734006,
10.2745538847873, 9.82381130225172, 9.36589007648649, 7.97974355540044,
10.9155195613281, 10.202304537148, 10.7215561563741, 8.46916191644844,
9.42098383221351, 10.2525138849946, 10.268546817668, 10.0951210882669,
10.4787506356114, 10.2231259271882, 9.48696504919052, 10.4217994558522,
9.22665448639598, 10.3202844245279, 8.77540023589803, 9.41379627008984,
10.611245796184, 10.291118982885, 10.1658334320981, 9.5097704294992,
9.80446942411316, 10.5187236924982, 10.728009962582, 9.80821864037385,
10.0408325161359, 9.65551599705282, 10.205052576579, 9.71199729976287,
10.8489948382187, 8.97538644201378, 10.5315320681001, 8.62855942629956,
10.6511914672025, 9.97265173885171, 9.9479622741747, 9.0461870197572,
10.8996252380421, 12.3599734384568, 10.7791131494668, 9.39083006285456,
9.27343378596625, 9.64671395390188, 10.4250180551466, 10.3860730834155,
9.59299070604241, 9.93941576249559, 10.2315399709258, 10.7241856094573,
10.3536466055067, 9.21930881655912, 10.0124320245946, 10.5817949676844,
9.78077184503572, 8.94773928584931, 8.95771617121254, 11.0943992108845,
9.67286494998046, 9.69010680414319, 8.55625741738904, 10.0996584351527,
10.0978820815018, 10.695949749899, 9.4659208384556, 9.00974731384131,
9.35904538363268, 11.1631409422, 10.0978949764737, 8.22147742199733,
10.1788503685018, 8.307445068494, 10.1030066323082, 9.54394535542576,
9.99365295701485, 9.62585483758286, 10.1725522249741, 9.86325708038225,
11.332094595289, 11.1220831427406, 10.2526884307178, 10.5079954715359,
9.76303390595084, 9.91645247809355, 10.4419322530252, 10.0010456173615,
9.14546634978685, 9.75386688753041, 9.93992493230075, 10.5397130560331,
11.1733004798793, 10.0440934469521, 9.82860126768683, 9.48840181330592,
10.2332175830011, 10.6422073305935, 10.2157323020761, 9.18660757619801,
10.5565745733978, 10.2652979905713, 10.6686993540023, 10.2757483111122,
9.28953270967829, 9.72926277144741, 10.6020489916687, 9.34008491043955,
10.2546876223525, 11.1845897000729, 10.3551925468931, 10.4353987795316,
9.33630309115289, 9.63047123935044, 9.41686283537383, 10.1720903474441,
9.62554868148692, 10.7024339853934, 10.0124089313601, 10.5701103985974,
9.62082523615912, 9.43601587467301, 10.0140427243452, 10.6577124935894,
10.1538626021251, 10.7650376693023, 10.3876485903888, 9.72499394983051,
10.5032974102856, 11.5942293340956, 9.56370482367918, 9.72817131853183,
9.61816990625114, 10.2500031094947, 9.97875025067029, 10.1505112693768,
9.19773878395889, 9.95359869051835, 10.2845730429639, 9.07792694868719,
9.99379614370394, 10.7058778542982, 9.75102309138852, 10.785528291343,
10.2539560650092, 10.0632774458996, 8.79411523093736, 9.15717026297178,
9.25564105724771, 9.01994559093434, 10.5636776737627, 10.4205734402798,
9.35132573872114, 9.67958072052341, 9.60768056642342, 11.357059786111,
10.0105427173501, 9.65443890298642, 9.5020194287417, 10.6696256729768,
9.2764884535406, 10.9867174910803, 10.3769784753715, 9.19910692367298,
9.68779970153859, 10.3200914799168, 10.6194858775346, 10.5610676478136,
10.2200784805803, 9.17732774137462)), row.names = c(NA, -424L
), class = c("tbl_df", "tbl", "data.frame"))


If I understand your modelling aims correctly, you want the residual variance to vary for the two methods? If that is the case, then you should simply predict sigma by method:

sigma ~ 0 + method


Can you clarify what you mean by wanting to keep the error term? Allowing unequal variance means that you have a separe resudual variance for each method and there is nothing left over to keep.

When it comes to the model you provided:

sigma ~ 1 + (0+method|method))


I’m not sure whether this has a meaningful interpretation. This means that sigma has an overall intercept, and then and a separate random effect for each method, that in turns varies within method. In other words, you are saying that the effect of methodnew on sigma is 0, but that it varies by whether the method is new or old. This, doesn’t quite make sense.

An additional issue is that your estimate for sigma is -0.34, but you are using -0.034 in your exp(…) statements. In those you are adding an intercept and subtracting from it the sd on the differences of each method from each method. This a bit like apples and oranges :)

I think the very simple approach above should achieve what you want. But if I have misunderstood your intention, it would help if you share your data simulation script - this should make it clear what you want to model and we can see how to adapt the formula to achieve it.

In essence I’d like to estimate the variance of method_new and method_old. On my head could be, u_{method= new} \sim N(0, \sigma_{method = new}) and u_{method= old} \sim N(0, \sigma_{method = old}) so the

sigma ~ 0 + method


could fit but I would like an overall error estimation as well, \epsilon_{year, subject, method}.