library(MASS) a = 10 beta_0 = 2 sd = .5 y = rnorm(12, mean = a, sd = sd) # the response x = rnorm(12, mean = a - beta_0, sd = .2) #case 1 ############################################################# ################# the second formulation x1_1<-c(rep(1,12)) x1_2<-x X1<-cbind(x1_1,x1_2) colnames(X1) = c("beta0", "beta1") beta1<-solve(t(X1)%*%X1)%*%t(X1)%*%y sd_beta1<-(sd*sd)*solve(t(X1)%*%X1) #the calculations in R: NB this is default in R fit.lm=lm(y~x) model.matrix(fit.lm) summary(fit.lm) #case 2 x.avg<-mean(x) x2_1<-c(rep(1,12)) x2_2<-x-x.avg X2<-cbind(x2_1,x2_2) colnames(X2) = c("a", "b") beta2<-solve(t(X2)%*%X2)%*%t(X2)%*%y sd_beta2<-(sd*sd)*solve(t(X2)%*%X2) #the calculations in R: NB this is default in R z<-x-x.avg fit.lm=lm(y~(z)) model.matrix(fit.lm) summary(fit.lm) ## S.t, s.t. X1=X2 x S.t S.t = ginv(X2)%*%X1 beta2.from.beta1 <- S.t%*%beta1 print(cbind(beta1,beta2,beta2.from.beta1)) # thus gamma.from.beta == gamma !!!