Hi fellow ComBat fans, Not sure if my previous message posted to the listserv because of the attachment. Now re-sending without attachment: I was having a problem with ComBat when using large batch sizes and the non-parametric parameter estimation. Almost all probe intensity values after ComBat correction would be NAs. I've traced the source of this problem and have a solution (switching from likelihood to log likelihood) described below. I hope you all find it helpful, or that you would let me know if you find any errors. Also, I can send the attachment which shows why I changed what I changed if anyone would like to take a look. ComBat correction ComBat was sometimes outputting NAs for all probe values using non-parametric parameter estimates. I traced this to the likelihood calculation in the int.eprior function. The original function is pasted below: # Monte Carlo integration function to find the nonparametric adjustments for a particular batch int.eprior <- function(sdat,g.hat,d.hat){ g.star <- d.star <- NULL r <- nrow(sdat) for(i in 1:r){ g <- g.hat[-i] d <- d.hat[-i] x <- sdat[i,!is.na(sdat[i,])] n <- length(x) j <- numeric(n)+1 dat <- matrix(as.numeric(x),length(g),n,byrow=T) resid2 <- (dat-g)^2 sum2 <- resid2%*%j LH <- 1/(2*pi*d)^(n/2)*exp(-sum2/(2*d)) LH[LH=="NaN"]=0 g.star <- c(g.star,sum(g*LH)/sum(LH)) d.star <- c(d.star,sum(d*LH)/sum(LH)) ##if(i%%1000==0){cat(i,'\n')} } adjust <- rbind(g.star,d.star) rownames(adjust) <- c("g.star","d.star") adjust } The likelihood calculation LH here goes to zero when sum2 is very large due to numerical precision errors. This causes the parameter estimates g.star and d.star to become NaN from a divide by zero error, and consequently for the output of the program to be NAs. This can happen when there are a lot of samples in a batch or when there are large residuals (or both). To correct this, can use the log-likelihood, which is less susceptible to numerical precision errors since it is not so small. I changed the function to use log-likelihood and the entire edited function is now pasted below: # Monte Carlo integration function to find the nonparametric adjustments for a particular batch int.eprior <- function(sdat,g.hat,d.hat){ #Get the number of probes r <- nrow(sdat); #The number of probes removing the probe of interest G = r - 1; #g.star and d.star are the output matrices g.star <- d.star <- matrix(NA,nrow=G,ncol=1); #Loop over each probe for(i in 1:r){ #Remove the beta value of the probe of interest g <- g.hat[-i]; #Remove the variance of the probe of interest d <- d.hat[-i]; #Isolate the value of this probe across all samples x <- sdat[i,!is.na(sdat[i,])]; #Number of samples in this batch n <- length(x) #A matrix of dimension gxn which contains the probe values repeated in each row dat <- matrix(as.numeric(x),length(g),n,byrow=T) #The residual between the probe values and the regression coefficient resid2 <- (dat-g)^2 #Calculate Log-likelihood using Gaussian LL = matrix(0,nrow=length(G),ncol=1); for (j in 1:n) { LL = LL + -0.5*log(2*pi*d)-(resid2[,j]/(2*d)); } #Order by the log-likelihood values to avoid Inf in parameter estimates orderind = order(LL,decreasing=TRUE); LL = LL[orderind]; g = g[orderind]; d = d[orderind]; #Determine the parameter estimates numerator.g = 1; denominator = 1; numerator.d = 1; for (k in 2:G) { numerator.g = numerator.g + (g[k]/g[1]) * exp(LL[k]-LL[1]); denominator = denominator + exp(LL[k]-LL[1]); numerator.d = numerator.d + (d[k]/d[1]) * exp(LL[k]-LL[1]); } g.star[i] = g[1]*(numerator.g/denominator); d.star[i] = d[1]*(numerator.d/denominator); } adjust <- rbind(g.star,d.star) rownames(adjust) <- c("g.star","d.star") adjust }