I am analyzing a RNA-seq data having some batch effect using DESeq2.

To account for the effect, I add it as a covariate to the regression model.

Next, to determine, how well the effect is corrected, I plot
two PCA plots

1. plotMDS(assays(dds)[["counts"]

] - raw count

2. plotMDS(assays(dds)[["mu"]] - fitted values.

My question is: Is it the correct way to check how well the regression
model worked on accounting for the batch effect?

No the "mu" here is not the right matrix to look at, as it contains the batch effect terms. In ?DESeq, you can rearrange the formula to get:

mu_ij = s_j 2^(beta_intercept + beta_batch + ...)

It's not easy to observe how the model corrects for batch effect, as it happens inside the model, in the balance of the coefficients from batch and the other variables. The model accounts for mean shift (in log common scale counts) which can be explained by batch.

You might look at the difference between the counts and fitted values on the log common scale to see if they associate with batch:

log(counts(dds, normalized=TRUE) + 1) - log(t( t(assays(dds)[["mu"]]) / sizeFactors(dds) ) + 1)