Hello,

I would like to test for differential methylation variability using missMethyl, but with a block design to take sample pairs/quartets into account.
If I wanted to test for differential (mean) methylation with my experimental design, I would use a code like this for example:

design <- model.matrix(~0+A1+A2+B1+B2)

blocks <- id
dupcor <- duplicateCorrelation(object=mval, design=design, block=blocks)
fit <- lmFit(mval, design, blocks, correlation=dupcor$consensus)

contrast.matrix <- makeContrasts(
  contrastA=A2-A1,
  contrastB=B2-B1,
  contrast1=B1-A1,
  contrast2=B2-A2,
  levels=design)

fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

topTable(fit2, coef="contrastA")
topTable(fit2, coef="contrastB")
topTable(fit2, coef="contrast1")
topTable(fit2, coef="contrast2")

This works fine.
However, when I try to "translate" this into code for missMethyl, the following happens:

> fitvar <- varFit(mval, design, blocks, correlation=dupcor$consensus)
Error in varFit(mval, design, blocks, correlation = dupcor$consensus) : 
  unused argument (correlation = dupcor$consensus)

I used blocks because I have samples from the same (or a subset of the same) individuals across all conditions and would like to take this into account. As the block design didn't work, I also tried to simply add the ids (which are factors) as covariates to the design, but then I got the following warning and no results (all NA):

> design <- model.matrix(~0+A1+A2+B1+B2+id)
> [...]
> fitvar <- varFit(mval, design=design)
Warning messages:
1: Partial NA coefficients for 466566 probe(s) 

Does anybody know how I could properly test for differential variability taking the information about the individuals into account?

Best,
Simone

> sessionInfo()
R version 3.2.1 (2015-06-18)
Platform: x86_64-unknown-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux Server 7.3 (Maipo)

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_GB.UTF-8   
 [6] LC_MESSAGES=en_GB.UTF-8    LC_PAPER=en_GB.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] limma_3.26.9     missMethyl_1.4.0 RSQLite_1.0.0    DBI_0.5-1   

I got the similar warning message when calling varFit under a contrast design matrix setting (multiple comparison groups).

> fitvar <- varFit(mval, design = design)
Warning messages:
1: Partial NA coefficients for 466566 probe(s)

The issue turns out to be a typo within the function varFit.default() in the source code DiffVar.R:

if(is.null(coef)) 
   coef <- c(1, ncol(design))
z <- getLeveneResiduals(data, design = design[,coef], type = type)

As you can see, the coefficients only take the first and the last comparison group. It should be coef <- 1:ncol(design).

I hope the authors can fix it in the next release but meanwhile you could just modify this line of code by yourself and re-install the package from the fixed source file. HTH.