For example, geneA ~ family + cur_env*pre_env.

• cur_env and pre_env are two categorical variables and both include CT and HD two levels.
  • I want to compare the fold change between CTCT (CT at pre_env and CT at cur_env) and HDHD (HD at pre_env and HD at cur_env).

I am interested in understanding how it does be compared behind the limma code and manually extract information of CTCT and HDHD and compute the fold change of their mu coefficients.

At the moment, I am using GAMLSS to replicate limma results by model interaction terms. My problem are

• How am I able to run similar contrast and obtain the coefficients and p value for CTCT and HDHD?
• Is it possible to construct reduced model to know whether combinations of two variables (CTCT, HDHD, CTHD, HDCT) have effects on response variables by taking interactions into account?

# design matrix
cur_env <- c(rep("CT", 4), rep("HD", 4))
pre_env <- c(rep("CT", 2), rep("HD", 2), rep("CT", 2), rep("HD", 2))
family <- c(rep("1", 3), rep("2", 2), rep("3", 3))

cur_env <- factor(cur_env, levels=c("CT","HD"))
pre_env <- factor(pre_env, levels=c("CT","HD"))
family  <- factor(family)

design <- model.matrix(~ family + cur_env*pre_env)

# make up coefficients
# (Intercept) family1 family2 family3  cur_envHD pre_envHD cur_envHD:pre_envHD 
 3.2       2.4       2.6      0.09      -0.14      0.09      -0.11

# design
  (Intercept) family2 family3 cur_envHD pre_envHD cur_envHD:pre_envHD
1           1       0       0         0         0                   0
2           1       0       0         0         0                   0
3           1       0       0         0         1                   0
4           1       1       0         0         1                   0
5           1       1       0         1         0                   0
6           1       0       1         1         0                   0
7           1       0       1         1         1                   1
8           1       0       1         1         1                   1
attr(,"assign")
[1] 0 1 1 2 3 4
attr(,"contrasts")
attr(,"contrasts")$family
[1] "contr.treatment"

attr(,"contrasts")$cur_env
[1] "contr.treatment"

attr(,"contrasts")$pre_env
[1] "contr.treatment"

The limma Users Guide explains how to test for interactions or to test for contrasts very easily. limma can work with any design matrix, but the easiest way by far is to combine pre_env and cur_env into one factor, and then use contrasts.fit to form the comparisons you want explicitly.

There is no need to form "reduced models" in limma. It is much simpler and more direct, you just form the comparisons you need directly.

At the moment, you don't seem to have tested any hypotheses in either package, limma or GAMLSS, so I don't understand what you mean by replicating results. You don't show any code from either package. GAMLSS is intended for joint modelling of the mean and scale parameter, but you don't mention why you think that might be relevant here. The design matrix you have formed doesn't seem to me to relate to your questions of interest at all closely.