stan_aov is completely messed up

[bgoodri]

Summary:

example(stan_aov) produces bogus output

Description:

example(stan_aov) yields

stan_aov
 family:       gaussian [identity]
 formula:      yield ~ block + N * P * K
 observations: 24
 predictors:   13
------
            Median        MAD_SD       
(Intercept)  5.490000e+01  9.000000e-01
block1      -1.302753e+14  4.523190e+15
block2      -4.342511e+13  1.507730e+15
block3      -2.171256e+13  7.538650e+14
block4       3.908260e+13  1.356957e+15
block5       2.605507e+13  9.046380e+14
N1           1.900000e+00  8.000000e-01
P1          -4.000000e-01  8.000000e-01
K1          -1.300000e+00  7.000000e-01
N1:P1       -6.000000e-01  7.000000e-01
N1:K1       -8.000000e-01  8.000000e-01
P1:K1        1.000000e-01  8.000000e-01
N1:P1:K1     1.302753e+14  4.523190e+15

This is probably due to changing the singular.ok logic. If you just do aov, there are only 12 coefficients estimated, rather than 13.

Reproducible Steps:

Run example(stan_aov, package = "rstanarm"))

RStanARM Version:

2.21.2

R Version:

4.1

Operating System:

Linux

[jgabry]

Hmm, that's not good. Are you saying the change to the singular.ok logic was wrong (it seemed correct, right?) or that stan_aov is doing the wrong thing with the correct singular.ok logic?

[bgoodri]

stan_aov was doing the correct thing before the change to singular.ok but possibly for the wrong reason

…

On Tue, Jul 21, 2020 at 4:59 PM Jonah Gabry ***@***.***> wrote: Hmm, that's not good. Are you saying the change to the singular.ok logic was wrong (it seemed correct, right?) or that stan_aov is doing the wrong thing with the correct singular.ok logic? — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#448 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZ2XKQWAKYSFVCXW7AGIA3R4X6ULANCNFSM4PEAAEVA> .

[jgabry]

I guess this can probably needs to be changed

rstanarm/R/stan_aov.R

Line 56 in 9c6d1cb

lmcall$singular.ok <- FALSE

to singular.ok <- TRUE and it will at least emulate the old behavior.

Just tried that and I get this, which looks better:

stan_aov
 family:       gaussian [identity]
 formula:      yield ~ block + N * P * K
 observations: 24
 predictors:   12
------
            Median MAD_SD
(Intercept) 52.7    2.6  
block2       2.4    2.5  
block3       4.7    2.6  
block4      -2.7    2.5  
block5      -2.4    2.5  
block6       1.6    2.5  
N1           6.8    2.6  
P1           0.3    2.5  
K1          -1.3    2.5  
N1:P1       -2.7    2.9  
N1:K1       -3.2    2.9  
P1:K1        0.4    2.8  

Auxiliary parameter(s):
              Median MAD_SD
R2            0.5    0.1   
log-fit_ratio 0.0    0.1   
sigma         4.2    0.7   

ANOVA-like table:
              Median MAD_SD
Mean Sq block 43.3   21.2  
Mean Sq N     48.1   28.0  
Mean Sq P     13.7   11.3  
Mean Sq K     28.7   19.5  
Mean Sq N:P   12.5   16.5  
Mean Sq N:K   16.6   21.5  
Mean Sq P:K    5.4    7.5  

If you're ok with that then let me know and I can push or just go ahead and make the change yourself.

[bgoodri]

Push it

…

On Tue, Jul 21, 2020 at 11:41 PM Jonah Gabry ***@***.***> wrote: I guess this can probably needs to be changed https://github.com/stan-dev/rstanarm/blob/9c6d1cb505e893c2d58ba57e63b09287908ad32d/R/stan_aov.R#L56 to singular.ok <- TRUE and it will at least emulate the old behavior. Just tried that and I get this, which looks better: stan_aov family: gaussian [identity] formula: yield ~ block + N * P * K observations: 24 predictors: 12 ------ Median MAD_SD (Intercept) 52.7 2.6 block2 2.4 2.5 block3 4.7 2.6 block4 -2.7 2.5 block5 -2.4 2.5 block6 1.6 2.5 N1 6.8 2.6 P1 0.3 2.5 K1 -1.3 2.5 N1:P1 -2.7 2.9 N1:K1 -3.2 2.9 P1:K1 0.4 2.8 Auxiliary parameter(s): Median MAD_SD R2 0.5 0.1 log-fit_ratio 0.0 0.1 sigma 4.2 0.7 ANOVA-like table: Median MAD_SD Mean Sq block 43.3 21.2 Mean Sq N 48.1 28.0 Mean Sq P 13.7 11.3 Mean Sq K 28.7 19.5 Mean Sq N:P 12.5 16.5 Mean Sq N:K 16.6 21.5 Mean Sq P:K 5.4 7.5 If you're ok with that then let me know and I can push or just go ahead and make the change yourself. — You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#448 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAZ2XKVVCXZYN6O4K6QP2HTR4ZNXTANCNFSM4PEAAEVA> .

[jgabry]

@bgoodri Unfortunately I think this warrants submitting a bugfix rstanarm to CRAN ASAP. They'll allow it if it's a major bug fix even though we released recently.

[jgabry]

@bgoodri Any progress on an rstanarm release? stan_aov has been messed up for 5 months now! We just got a stan_aov question on the forums

https://discourse.mc-stan.org/t/stan-aov-is-confusing-me/19841

but installing rstanarm from github isn't a general purpose solution for everyone since that requires the toolchain, compilation, etc.

[jgabry]

Also @bgoodri if you have a chance can you answer the stan_aov question on the forums I linked to in my previous post? It's about more than the bug. Thanks!