GWAS with phenotype adjusted for covariate (age, sex, PC, batch, centre) using PLINK or external linear model and take residual?

Hi everyone, I am conducting GWAS analysis using PLINK but I wonder should I use PLINK to adjust the phenotype for fix effect (age, sex, PC, batch, centre) or should I use a separate regression (such as lm(pheno ~ covar) in R) and get the residuals for the GWAS?

Also, how about for Binary trait (case/control)? the residuals from logistic is not interpreted the same way as the residuals from the linear model.

And how about interaction effect (GxE) for binary cases, I know how to use the –parameters option to select the variables of interest but I wonder how to conduct GWAS adjusted for fix effect properly.

Using PLINK considering covariate / fix effects take quite a long time, I just have been able to run with the residuals getting from R.

Your help is really appreciated!