I am comparing expression values in paired samples given a continuous variable. In theory, I should see an increase in expression given an increase in the total score between the paired samples.

This is the code that I have generated to run the analysis. I am using the section "9.4.1 Paired Samples" of the limma user guide as a reference.

ID <- factor(mSetSqFlt$Collaborator.ID)

Visit <- factor(mSetSqFlt$TimePoint, levels = c("Baseline","Followup"))

Score <- factor(mSetSqFlt$Total.Score)


design <- model.matrix(~ID + Visit + Visit:Score)

fit <- lmFit(mVals, design)

fit2 <- eBayes(fit)

topTable(fit2, coef = "VisitFollowup")

I get outputs using this method and would just like some insights into if I am on the right track.

Any insights would be very much appreciated as I am new to limma.

You are including too many terms in the linear model. TimePoint and Total.Score give almost the same information. You can include either one of them in the model but not both.

The model you have defined in your question is recommended in the User's Guide for time-course experiments in independent groups, rather than for paired experiments. Defining an interaction is not at all relevant in the paired comparison context.