'singular.ok' logic is reversed in 'stan_lm'

rvlenth commented 4 years ago

Summary:

It appears that in stan_lm(), singular.ok does the opposite of what it is supposed to do. I think the logic is backwards somewhere.

Description:

In the example below, I deleted a cell from the data which creates a singularity. With singular = TRUE, stan_lm() exits with an error while singular = FALSE, a coefficient is discarded. This is the opposite of what is done in lm().

Reproducible Steps:

> packageVersion("rstanarm")
[1] ‘2.19.2’

> ### stan_lm tests ###
> stan_lmT <- stan_lm(breaks ~ wool*tension, data = warpbreaks, 
+     subset = -(9:18), seed = 1776, prior = R2(0.5), singular.ok = TRUE)
Error in backsolve(R, diag(K)) : 
  singular matrix in 'backsolve'. First zero in diagonal [5]

> stan_lmF <- stan_lm(breaks ~ wool*tension, data = warpbreaks, 
+     subset = -(9:18), seed = 1776, prior = R2(0.5), singular.ok = FALSE)
> coef(stan_lmF)
   (Intercept)          woolB       tensionM       tensionH woolB:tensionH 
    40.1165632    -11.8845024      0.5421783    -15.2374347      6.9444381

> ### stats::lm tests -- just the opposite ###
> stats_lmT <- lm(breaks ~ wool*tension, data = warpbreaks, 
+     subset = -(9:18), singular.ok = TRUE)
> coef(stats_lmT)
   (Intercept)          woolB       tensionM       tensionH woolB:tensionM woolB:tensionH 
    41.7500000    -13.5277778      0.5555556    -17.1944444             NA      7.7500000 

> stats_lmF <- lm(breaks ~ wool*tension, data = warpbreaks, 
+     subset = -(9:18), singular.ok = FALSE)
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
  singular fit encountered

> ### stan_aov has no singular.ok but discards variables when singular ###
> stan_aov_ <- stan_aov(breaks ~ wool*tension, data = warpbreaks, 
+     subset = -(9:18), seed = 1776, prior = R2(0.5))
> coef(stan_aov_)
   (Intercept)          woolB       tensionM       tensionH woolB:tensionH 
    40.1165632    -11.8845024      0.5421783    -15.2374347      6.9444381

RStanARM Version:

2.19.2

R Version:

3.6.1

Operating System:

Windows 10 v 1903

rvlenth commented 4 years ago

Just a note that the singularity issue in stan_lm still does not work correctly.

> lm.mod = stan_lm(breaks ~ wool*tension, data = warpbreaks[1:44,])
Error in backsolve(R, diag(K)) : 
  singular matrix in 'backsolve'. First zero in diagonal [5]

# Try explicitly specifying singular.ok...
> lm.mod = stan_lm(breaks ~ wool*tension, data = warpbreaks[1:44,], singular.ok = TRUE)
Error in backsolve(R, diag(K)) : 
  singular matrix in 'backsolve'. First zero in diagonal [5]

> lm.mod2 = stan_lm(breaks ~ wool*tension, data = warpbreaks[1:44,], singular.ok = FALSE)
Error in stopifnot(is.numeric(location)) : 'location' must be specified

Meanwhile, using stan_glm() does handle this right:

> glm.mod = stan_glm(breaks ~ wool*tension, data = warpbreaks[1:44,])
Warning message:
In center_x(x, sparse) : Dropped empty interaction levels: woolB:tensionH

> library(emmeans)
> emmeans(glm.mod, ~ wool * tension)
 wool tension emmean lower.HPD upper.HPD
 A    L         43.6      35.8      51.1
 B    L         28.7      19.9      36.4
 A    M         24.5      16.7      32.7
 B    M         28.5      19.9      36.8
 A    H         24.8      16.6      33.1
 B    H       nonEst        NA        NA

Point estimate displayed: median 
HPD interval probability: 0.95

I'm unclear when this issue was closed. If I did that, it was a mistake. But in any case, thanks for your attention.

jgabry commented 4 years ago

Thanks for following up about this and sorry no one fixed this earlier! I think you're right that the logic is reversed. I think

https://github.com/stan-dev/rstanarm/blob/887dbdd00770baac27b17ce32c111fde3f94f3fe/R/stan_lm.fit.R#L63

should say singular.ok instead of !singular.ok. I'll push a commit fixing this in a few minutes.

rvlenth commented 4 years ago

Thanks. I agree that this will appear to fix it; but I'll test it.

jgabry commented 4 years ago

Thanks. I tested it before closing the issue and it worked well. If it doesn't fix it for you let me know and I'll dig deeper.

rvlenth commented 4 years ago

Great -- I figured you had. That may have to do for now, as I was unable to build it for myself. I think it was actually the StanHeaders package I couldn't build. Here are the last few lines of a very long stream from the build process:

C:/Users/rlenth/Documents/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:
C:/Users/rlenth/Documents/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:13: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]
 static void set_zero_all_adjoints() {
             ^~~~~~~~~~~~~~~~~~~~~

cc1plus.exe: out of memory allocating 14907328 bytes
make: *** [C:/PROGRA~1/R/R-40~1.1/etc/i386/Makeconf:227: stan_files/jm.o] Error 1
rm stan_files/jm.cc stan_files/binomial.cc stan_files/lm.cc stan_files/continuous.cc stan_files/bernoulli.cc stan_files/polr.cc stan_files/count.cc stan_files/mvmer.cc
ERROR: compilation failed for package 'rstanarm'
* removing 'C:/Users/rlenth/Documents/R/win-library/4.0/rstanarm'
* restoring previous 'C:/Users/rlenth/Documents/R/win-library/4.0/rstanarm'
Warning messages:
1: In utils::untar(tarfile, ...) :
  ‘tar.exe -xf "C:\Users\rlenth\AppData\Local\Temp\Rtmpsv7k0D\file1ce4740d1102.tar.gz" -C "C:/Users/rlenth/AppData/Local/Temp/Rtmpsv7k0D/remotes1ce42479d7a"’ returned error code 1
2: In i.p(...) :
  installation of package ‘C:/Users/rlenth/AppData/Local/Temp/Rtmpsv7k0D/file1ce4460b74da/rstanarm_2.19.3.tar.gz’ had non-zero exit status

This may relate to a memory-allocation bug in R 4.0.1 that Peter Dahlgaard mentioned in a June 9 e-mail to r-announce.

bgoodri commented 4 years ago

That is usually due to it picking up RTools35 rather than RTools40.

On Wed, Jun 10, 2020 at 2:50 PM Russell V. Lenth notifications@github.com wrote:

Great -- I figured you had. That may have to do for now, as I was unable to build it for myself. I think it was actually the StanHeaders package I couldn't build. Here are the last few lines of a very long stream from the build process:

C:/Users/rlenth/Documents/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp: At global scope:

C:/Users/rlenth/Documents/R/win-library/4.0/StanHeaders/include/stan/math/rev/core/set_zero_all_adjoints.hpp:14:13: warning: 'void stan::math::set_zero_all_adjoints()' defined but not used [-Wunused-function]

static void set_zero_all_adjoints() {
         ^~~~~~~~~~~~~~~~~~~~~
cc1plus.exe: out of memory allocating 14907328 bytes

make: *** [C:/PROGRA~1/R/R-40~1.1/etc/i386/Makeconf:227: stan_files/jm.o] Error 1

rm stan_files/jm.cc stan_files/binomial.cc stan_files/lm.cc stan_files/continuous.cc stan_files/bernoulli.cc stan_files/polr.cc stan_files/count.cc stan_files/mvmer.cc

ERROR: compilation failed for package 'rstanarm'

removing 'C:/Users/rlenth/Documents/R/win-library/4.0/rstanarm'

restoring previous 'C:/Users/rlenth/Documents/R/win-library/4.0/rstanarm'

Warning messages:

1: In utils::untar(tarfile, ...) :

‘tar.exe -xf "C:\Users\rlenth\AppData\Local\Temp\Rtmpsv7k0D\file1ce4740d1102.tar.gz" -C "C:/Users/rlenth/AppData/Local/Temp/Rtmpsv7k0D/remotes1ce42479d7a"’ returned error code 1

2: In i.p(...) :

installation of package ‘C:/Users/rlenth/AppData/Local/Temp/Rtmpsv7k0D/file1ce4460b74da/rstanarm_2.19.3.tar.gz’ had non-zero exit status

This may relate to a memory-allocation bug in R 4.0.1 that Peter Dahlgaard mentioned in a June 9 e-mail to r-announce.

— You are receiving this because you were assigned. Reply to this email directly, view it on GitHub https://github.com/stan-dev/rstanarm/issues/402#issuecomment-642193508, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZ2XKTCYSK3M7P3RHWBAGDRV7IYPANCNFSM4JJMAXMA .

rvlenth commented 4 years ago

I do have rtools40 installed. I renamed the old rtools directory just to make sure; but re-trying the build did not reveal errors about not finding stuff, just the same errors.

stan-dev / rstanarm