stan_glm model with no intercept gives unrealistic large residual variance

[f ko]

Dears
thanks for rstanarm! This is wonderful - exactly what we (ecologists) need! Really great.

I am sorry to bother you with maybe a very stupid question. I am interested in the difference between males and females in mean wing length of Snowfinches.
I  have only few data, but in addition I have for each group a mean with a SE from the literature that I would like to use as prior information.
When I fit two indepenent models for each group  and use a normal prior for the intercept I get a reasonable posterior distribution. However, when I fit the model to both groups separately excluding the intercept so that the two parameters correspond to the group means, the resulting residual variance is way too large (see r-code below). I could not find the reason for this. In my opinion the sigma in the combined analyses should be something between 2.5 and 5.2 which are the group-specific residual standard deviations. Because I am interested in the difference between the two means, I need to fit the combined model to obtain the joint posterior distribution of the two means.  
I am happy about any hint on what I do wrong. Thanks so much in advance!
Fränzi


# some data
group <- c("M", "M", "M", "W", "W", "W", "W", "W", "W")
y <- c(125, 122, 121, 119, 119.5, 116, 121, 121, 117.5)
dat <- data.frame(group=group, y=y)


# values from the literature (m SE)
mm <- 123.8
semm <- 0.381

mw <- 116.7
semw <- 0.387


# fit a model to the data without prior information
mod <- lm(y~group-1)
mod
# Coefficients:
#   groupM  groupW 
# 122.7   119.0
summary(mod)$sigma
#[1] 2.005944


# include prior information for the two means using rstanarm
library(rstanarm)
stanmod <- stan_glm(y~group-1, data=dat,
                    prior = normal(c(mm, mw), scale=c(semm, semw)),
                    prior_intercept = NULL)
print(stanmod)
# Estimates:
#   Median MAD_SD
# groupM 123.8    1.9
# groupW 116.7    2.0
# sigma  118.0   27.3    THIS SIGMA IS OVBIOUSLY WAY TOO LARGE

stanmod <- stan_glm(y~1, data=dat[dat$group=="M",],
                    prior = NULL,
                    prior_intercept = normal(mm, scale=semm))

print(stanmod)
# Estimates:
#   Median MAD_SD
# (Intercept) 123.7    0.4
# sigma         2.5    1.1  THIS SIGMA SEEMS REASONABLE



stanmod <- stan_glm(y~1, data=dat[dat$group=="W",],
                    prior = NULL,
                    prior_intercept = normal(mm, scale=semm))

print(stanmod)
# Estimates:
#   Median MAD_SD
# (Intercept) 123.6    0.4
# sigma         5.2    1.5   THIS SIGMA SEEMS REASONABLE

[Ben Goodrich]

On Saturday, March 5, 2016 at 4:53:11 AM UTC-5, f ko wrote:

Dears
thanks for rstanarm! This is wonderful - exactly what we (ecologists) need! Really great.

What you really need is rstanarm 2.9.0-3, which fixed a bug in the first rstanarm involving models without intercepts. Just do

update.packages("rstanarm")

Ben

[Paul Buerkner]

I may be mistaken but it appears to me that the bug is still there even with rstanarm 2.9.0-3 (CRAN version). Using the sample data that Fränzi provided, I get the same very large estimates for sigma.

Paul

[Ben Goodrich]

On Sunday, March 6, 2016 at 6:32:31 AM UTC-5, Paul Buerkner wrote:

I may be mistaken but it appears to me that the bug is still there even with rstanarm 2.9.0-3 (CRAN version). Using the sample data that Fränzi provided, I get the same very large estimates for sigma.

I get

> print(stanmod)
stan_glm(formula = y ~ group - 1, data = dat, prior = normal(c(mm,
    mw), scale = c(semm, semw)), prior_intercept = NULL)

Estimates:
       Median MAD_SD
groupM 123.0    1.0
groupW 118.6    0.9
sigma    2.2    0.6

Sample avg. posterior predictive
distribution of y (X = xbar):
         Median MAD_SD
mean_PPD 120.1    1.0

What does

t(file.info(system.file("libs", "rstanarm.so", package = "rstanarm")))

say for you?

Ben

[Paul Buerkner]

t(file.info(system.file("libs", "rstanarm.so", package = "rstanarm")))

is NA. When I write

t(file.info(system.file("libs", package = "rstanarm")))

I get
     C:/Users/Paul/Documents/R/win-library/3.2/rstanarm
size  "0"                                              
isdir "TRUE"                                           
mode  "777"                                            
mtime "2016-03-06 12:28:50"                            
ctime "2016-03-06 12:28:50"                            
atime "2016-03-06 12:28:50"                            
exe   "no"
 

I have updated rstanarm and also tried to remove it completely and install it again. Still the same. Maybe I am missing some important thing. If it works for the OP, I think it doesn't really matter if it works for me.

[Ben Goodrich]

On Sunday, March 6, 2016 at 1:10:23 PM UTC-5, Paul Buerkner wrote:

t(file.info(system.file("libs", "rstanarm.so", package = "rstanarm")))

is NA.

I think it would be rstanarm.dll on Windows, or in general

t(file.info(system.file("libs", paste0("rstanarm", .Platform$dynlib.ext), package = "rstanarm")))

Ben

[Paul Buerkner]

Still NA. Tomorrow, I will check it on another Windows machine (both the model and the file.info).

Paul

[Paul Buerkner]

After testing it on another Windows machine (with no previous rstanarm installation), I can confirm that the latest CRAN version (2.9.0-3) still produces those inflated sigma estimates. Unfortunately, the file.info is also NA on the other machine.
I can't test it with the current dev version from github at the moment, because my Windows refuses to create the DLL (likely a problem on my side though). Based on your model results, I think it's save to assume that the dev version would work fine also on Windows. Again, not a big deal for me and since the OP didn't respond again, your suggestions may have solved the issue for him.

Paul

[f ko]

Dear Ben and Paul
thanks for treating my question. And sorry, for my late response.
I updated rstanarm and I still get the inflated sigma. Thus it did not help.

And writing the code suggested by Ben I get a lot of NAs

> t(file.info(system.file("libs", paste0("rstanarm", .Platform$dynlib.ext), package = "rstanarm")))

        
size  NA
isdir NA
mode  NA
mtime NA
ctime NA
atime NA
exe   NA

I also work on a Windows machine.
Any hint on what I could try is welcome!
Thanks a lot
Fränzi

[Ben Goodrich]

On Monday, March 7, 2016 at 9:18:41 AM UTC-5, f ko wrote:

thanks for treating my question. And sorry, for my late response.
I updated rstanarm and I still get the inflated sigma. Thus it did not help.

It works for me on Windows, which it should because the intercept-handling code is not platform specific. The only thing I can think of is that you are using an older version of R (like 3.1.x) and it is downloading 2.9.0-1 rather than 2.9.0-3.

Try these:

packageVersion("rstanarm")

ip <- installed.packages()
ip["rstanarm",]

t(file.info(system.file("libs", .Platform$r_arch, package = "rstanarm")))

Ben

[Paul Buerkner]

Here's the output:

packageVersion("rstanarm")

[1] ‘2.9.0.3’

ip <- installed.packages()
ip["rstanarm",]

                                                                                                                                   Package
                                                                                                                                "rstanarm"
                                                                                                                                   LibPath
                                                                                               "C:/Users/Paul/Documents/R/win-library/3.2"
                                                                                                                                   Version
                                                                                                                                 "2.9.0-3"
                                                                                                                                  Priority
                                                                                                                                        NA
                                                                                                                                   Depends
                                                                                                 "R (>= 3.2.0), Rcpp (>= 0.11.0), methods"
                                                                                                                                   Imports
"ggplot2 (>= 2.0.0), lme4 (>= 1.1-8), loo (>= 0.1.4), Matrix,\nnlme (>= 3.1-124), rstan (>= 2.9.0-3), shinystan (>= 2.1.0),\nstats, utils"
                                                                                                                                 LinkingTo
                                                "StanHeaders (>= 2.9.0), rstan (>= 2.9.0-3), BH (>= 1.58.0),\nRcpp (>= 0.11.0), RcppEigen"
                                                                                                                                  Suggests
                                                            "arm, gamm4, gridExtra, HSAUR3, KernSmooth, knitr, MASS,\nrmarkdown, testthat"
                                                                                                                                  Enhances
                                                                                                                                        NA
                                                                                                                                   License
                                                                                                                              "GPL (>= 3)"
                                                                                                                           License_is_FOSS
                                                                                                                                        NA
                                                                                                                     License_restricts_use
                                                                                                                                        NA
                                                                                                                                   OS_type
                                                                                                                                        NA
                                                                                                                                    MD5sum
                                                                                                                                        NA
                                                                                                                          NeedsCompilation
                                                                                                                                     "yes"
                                                                                                                                     Built
                                                                                                                                   "3.2.3"

t(file.info(system.file("libs", .Platform$r_arch, package = "rstanarm")))

      C:/Users/Paul/Documents/R/win-library/3.2/rstanarm/libs/x64

size  "0"                                                       
isdir "TRUE"                                                    
mode  "777"                                                     

mtime "2016-03-07 12:31:21"                                     
ctime "2016-03-07 12:31:21"                                     
atime "2016-03-07 12:31:21"                                     
exe   "no"  

[Ben Goodrich]

On Monday, March 7, 2016 at 3:46:13 PM UTC-5, Paul Buerkner wrote:

Here's the output:

This is very odd. If you do rstan::get_stancode() on the $stanfit element of the list that is returned by the model with the crazy variance does it have code like

if (has_intercept == 1) {
  // stuff
}
else { // this was added in 2.9.0-3
  eta <- eta + dot_product(xbar, beta);
}

Ben

[Paul Buerkner]

yep:

 

vector[N] eta;  // linear predictor
  if (K > 0) eta <- X * beta;
  else eta <- rep_vector(0.0, N);
  if (has_offset == 1) eta <- eta + offset;
  if (t > 0) eta <- eta + csr_matrix_times_vector(N, q, w, v, u, b);
  if (has_intercept == 1) {
    if (family == 1 || link == 2) eta <- eta + gamma[1];
    else eta <- eta - min(eta) + gamma[1];
  }
  else {
  // correction to eta if model has no intercept (because X is centered)

  eta <- eta + dot_product(xbar, beta);
  }

It feels like I have made some sort of (very) stupid mistake, but I currently don't see where...

[f ko]

I use R 3.2.3 but indeed it has downloaded rstanarm 2.9.0-1 when I updated rstanarm from within RStudio. Now I updated rstanarm not using RStudio and it downloaded rstanarm 2.9.0-3 AND IT WORKS AND GIVES THE SAME RESULT AS BEN HAS GOT!  Wonderful!
Thanks so much for this help!
Fränzi

[Bob Carpenter]

Hi, Fränzi. Ben --- is this just slow mirroring on
Rstudio's part? Is it possible to set a different mirror
so you can get more up to date CRAN pacakges?

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Ben Goodrich]

On Tuesday, March 8, 2016 at 4:14:41 PM UTC-5, Bob Carpenter wrote:

Hi, Fränzi.  Ben --- is this just slow mirroring on
Rstudio's part?  Is it possible to set a different mirror
so you can get more up to date CRAN pacakges?

I don't know. I wouldn't think any mirror would be weeks behind the mothership. I would just always use https://cloud.r-project.org as the repository which will redirect to whatever is the closest reliable mirror that has the latest version of the package.

Ben

[Paul Buerkner]

Ok, I found out why the results differed. When I replicate the script of the OP exactly (including priors) I get the same  results as Ben.

However, when I just call

stanmod2 <- stan_glm(y~group-1, data=dat)

I get some weired sigma estimates again. Also, the other parameters seem to be a bit off:

> print(stanmod2)
stan_glm(formula = y ~ group - 1, data = dat)

Estimates:
       Median MAD_SD
groupM   5.1   13.0
groupW  10.0   12.9
sigma  110.3   26.0

Sample avg. posterior predictive
distribution of y (X = xbar):
         Median MAD_SD

mean_PPD  9.8   37.7 

[Ben Goodrich]

On Saturday, March 12, 2016 at 12:08:20 PM UTC-5, Paul Buerkner wrote:

Ok, I found out why the results differed. When I replicate the script of the OP exactly (including priors) I get the same  results as Ben.

However, when I just call

stanmod2 <- stan_glm(y~group-1, data=dat)

I get some weired sigma estimates again. Also, the other parameters seem to be a bit off:

I believe this is due to the default priors of stan_glm being quite inconsistent with the data-generating process assumed by this model. If you nullify all the priors, you get something similar to the MLEs

stan_glm(y~group-1, data=dat, prior = NULL, prior_intercept = NULL, prior_ops = prior_options(scaled = FALSE))

Ben