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,
- the model error term and
- 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
- Total sigma (
sigma_Intercept
) - 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"))
Thanks for your time.