Hello all,

I was using limma for assessing the deferentially expressed genes in my microarray experiment. However, during the analysis process I've encountered the following error

Error in chol.default(Zinfo) : 
  the leading minor of order 11 is not positive definite
Calls: arrayWeights -> chol -> chol.default
Execution halted

I guess this error is due to a larger design matrix. Here it is

design = model.matrix(~0+factor(c(1,2,3,1,2,3,7,8,9,10,11,12)))
colnames(design) = c("CDH", "CDN", "CDP", "CDPCIS", "CDPFU", "CDPUN", "CDNCIS", "CDNFU", "CDNUN")
aw = arrayWeights(data.filt,design)
fit = lmFit(data.filt,design,weights=aw) 

I would be grateful for your help.

Thank you.

This error is caused by the presence of several groups that contain only a single array. arrayWeights works by identifying noisy arrays with large residuals from the linear model, and downweighting them accordingly. However, for the groups with only one array, the residual for that array will always be zero, i.e., the linear model fits it perfectly for all genes. I would guess that this results in too many zeroes in whatever matrix is being used in the Cholesky decomposition, resulting in the observed error about the matrix not being positive definite.

There's no point computing weights for the arrays in single-array groups, as the choice of weight won't affect the fit. Each of these arrays is always fitted perfectly and doesn't contribute to the variance estimate. So, the simplest workaround would be to estimate weights for only the arrays in groups 1-3, where the choice of weight does affect the variance estimate; and set the other weights to unity (or whatever non-zero value you desire).

Perhaps arrayWeights could also do this automatically, i.e., return unity weights for arrays that don't contribute to the residual d.f. of the fit. This might be more elegant than dropping an error, as it would not impede progress through to the rest of the analysis.

There's no point in using arrayWeights() unless you have groups with at least 3 replicates. With fewer than 3 replicates, it is impossible to identify unusual cases because it is impossible for any single case to stand out.

Just omit the arrayWeights() step and proceed with the analysis as normal. As Aaron says, groups with just n=1 replicate are not ideal, but limma will do the best that can be done, which is to use the replicates that are available.