Experimental design before and after with control group I have a project which test effectiveness of two factors (two levels, receive and not receive the factor). There are three stage consecutively through time, and the factors are imposed in the second stage. Let's call the first stage baseline period, the second stage short term test, the third stage long term test.    
Question is how to test the effectiveness of the two factors? What kind of statistical method would you recommend?   
I'm trying to use multiple linear regression with interaction term, while the before and after with control group through three stages in time series confused me on how to interpret the results?   
 A: Multivariate regression sounds like a good fit for your problem in principle. However, you should be careful as to how you compute your standard errors, as you have panel data. Two common solutions for modeling panel data are mixed models and generalized estimating equations. With a mixed model you must specify a model for dependence across time and failing to specify a correct dependence model may lead to incorrect point estimates/standard errors. 
Consequently, as long as there is a good sample size available I prefer to use generalized estimating equations which produce consistent estimates and consistent standard errors. Here is an example for how to use generalized estimating equations. 
n <- 500
treatment <- rbinom(n, 1, 0.5)
datlist <- list()
library(mvtnorm)
sigma <- matrix(c(1, 0.75, 0.5, 0.75, 1, 0.75, 0.5, 0.75, 1), ncol = 3) * 0.3
for(i in 1:length(treatment)) {
  datlist[[i]] <- data.frame(treatment = treatment[i], timepoint = 0:2, ptid = i,
                             rintercept = as.numeric(rmvnorm(3, sigma = sigma)))
}
simdata <- do.call("rbind", datlist)

X <- model.matrix(~ treatment * factor(timepoint), data = simdata)
colnames(X)
beta <- c(0, 0, 0, 0, 1, 0.1)
mu <- as.numeric(X %*% beta)
y <- rnorm(nrow(simdata), mu, 1) + simdata$rintercept

library(geeM)
simdata$timepoint <- factor(simdata$timepoint)
geefit <- geem(y ~ treatment * timepoint, data = simdata,
               id = ptid, waves = as.integer(simdata$timepoint))
summary(geefit)

It is important to note that in generalized linear models, GEEs and mixed models estimate different parameters. However, in linear regression, the parameters estimated by the two methods are the the same and the difference is in how efficient/robust the methods are. Mixed effects models are more efficient if the variance model is specified correctly and the GEEs are less efficient but more robust to model misspecification.  
