Refitting ash using mixsqp

[william-denault]

Hello,

I am trying to use mixsqp to solve some weighted EM problems. However, I am quite puzzled by the results and I noticed that I was unable to use mixsqp to recover the weight found by ash. Please a minimal example below of the problem. I think I specify the matrix L correctly. I would be glad if you could let me know what am I missing here.

set.seed(1)
library(ashr)
library(mixsqp)

### Generate the data under the nul
n_coef <- 10
sd_hat <- runif(n_coef , 0.1, 0.1)
beta_hat <- rnorm(n_coef , mean = rep(0,n_coef ), sd= sd_hat)

m <- ash(beta_hat, sd_hat, mixcompdist="normal")
m$fitted_g$sd
m$fitted_g$pi 

# Trying to recover the pi vector using mixsqp

L <- matrix( NA, nrow=  length(beta_hat), ncol=length(m$fitted_g$pi))

for (k in 1:dim(L)[2]){
  L[, k ] <- c(dnorm( beta_hat, mean=rep(0, length(beta_hat)),
                      sd = sqrt(m$fitted_g$sd[k]^2 + sd_hat^2),
                      log = log))  
}

t_pi <- mixsqp(L , log=TRUE )$x
t_pi
m$fitted_g$pi

# Piggybacking on ashr code
sdmat <-sqrt(outer(sd_hat^2,m$fitted_g$sd^2,"+"))
matrix_llik = (dnorm(outer(beta_hat,m$fitted_g$mean,FUN="-")/sdmat,log=TRUE) -log(sdmat )) # an n by k matrix

image(matrix_llik - L) #it seems that the L matrix was correctly specified
mixsqp(matrix_llik, log=TRUE )$x
t_pi #previous estimates
m$fitted_g$pi

[pcarbo]

@william-denault A few small comments on sharing code (for the next time): (1) use the GitHub Markdown features to format your code so that it is easier to read; (2) please provide the output you get so I can compare with what I get; (3) add a sessionInfo() call so I know what version of R and package versions you are using.

[pcarbo]

Try this:

library(ashr)
library(mixsqp)
set.seed(1)

# Generate the data under the null.
n_coef <- 10
sd_hat <- rep(0.1,n_coef)
beta_hat <- rnorm(n_coef,mean = rep(0,n_coef),sd = sd_hat)

m <- ash(beta_hat,sd_hat,mixcompdist = "normal",nullweight = 1,outputlevel = 3)

# Trying to recover the pi vector using mixsqp.
L <- matrix(as.numeric(NA),length(beta_hat),length(m$fitted_g$pi))

for (k in 1:ncol(L))
  L[,k] <- dnorm(beta_hat,
                 mean = rep(0,length(beta_hat)),
                 sd = sqrt(m$fitted_g$sd[k]^2 + sd_hat^2),log = TRUE)

t_pi <- mixsqp(L,log = TRUE,control = list(eps = 1e-6,numiter.em = 20))$x
plot(m$fitted_g$pi,t_pi,pch = 20)
abline(a = 0,b = 1,col = "skyblue",lty = "dotted")

When I ran this I got the same mixture proportions. There are a few subtleties here, resulting in different maximum-likelihood estimates of the mixture proportions in small examples such as yours.

Also outputlevel = 3 isn't needed here but useful for debugging.

[william-denault]

Thank you very much for your feedback. I was having a hard time understanding my "mistake". Is there any other option/ control I should specify to some given value to get "robust" estimate was using small examples (which is actually quite critical in my setting).

[pcarbo]

William, it may be better to discuss this on Slack so I understand your problem better. You might find that it isn't all that important to get a "robust" estimate of the mixture weights (pi), but let's discuss.