For theoretical develpments, refer
In R, linear models can be directly implemented using the base function lm( ) . The set of arguments / inputs include a formula connecting response variable (numeric) and set of predictors; that is of the form $ y \sim x_1 + x_2 + \cdots x_k$. Intercept or constant term will be directly added to the model matrix.
Further details about lm( ) can be found in the official reference pages
From the theoretical reference mentioned, it can be observed that the OLS solution is basically a “closed form” solution and all the computations are arithmetic and matrix operations. Hence, implementing in a programming environment may not be very difficult task; here, is one such reproducible coding steps
#An artificial data set
x1=c(37,45,38,42,31)
x2=c(4,0,5,2,4)
y1=c(71,66,75,70,65) #Response variable
#-------------------------------------------------
#Necessary Inputs
k=2 #No of regressors
n=length(x1) #No of data points
p = k+1
alp = 0.05 #Alpha for (1-Alpha)100% CI
#-------------------------------------------------
#Given Predictors in matrix form
xd=matrix(c(x1,x2),n,k)
#pairs(xd)
xo=matrix(1,n,1)
x=cbind(xo,xd) #Modal matrix
m=t(x)%*%x #X'X
d=det(m) #Determinant of m to chk singularity
y=as.matrix(y1,n,1)
#-------------------------------------------------
#OLS Estimations - Det(model matrix) = 0 or not
if (d!=0)
{
print("Matrix X'X is NOT singular")}else
print("Matrix X'X is singular")
#-------------------------------------------------
#OLS Estimations
#mi is inverse of X'X
#mi=solve(m)
#-------------------------------------------------
#OLS Estimations - Regression beta coeffts
betao = mi%*%t(x)%*%y
#-------------------------------------------------
#Calculations based on OLS
#Residuals
#Sum of Squares_Residue
#MSReg = estimate of SigmaSq
#Covariance of reg coeffts (Betas)
SSRes=as.vector(t(y)%*%y-t(betao)%*%t(x)%*%x%*%betao)
MSRes=SSRes/(n-k-1)
covao=MSRes*mi
cbind(SSRes,MSRes)
#-------------------------------------------------
#Regression
#Sum of Squares_Regression
SSR=t(betao)%*%t(x)%*%y-(((sum(y))^2)/n)
MSR=SSR/k #Mean Square_Regression
#Total
SST=SSRes+SSR
#-------------------------------------------------
#Overall Adequacy
F_O = MSR/MSRes
pval_F_O = 1-pf(F_O,k,n-p)
Rsquare_O = SSR/(SSR+SSRes)
AdjRsquare_O = 1-((SSRes/(n-k-1))/((SSR+SSRes)/(n-1)))
#-------------------------------------------------
#Standard error for individual coeffts
dia_elt=0 #To pick diagonal elements of Covariance matrix - covao
for (i in 1:p)
{
dia_elt[i]=sqrt(covao[i,i])
}
#t-statistic, P values
calc_t=betao/dia_elt
pval_O=2*pt(-abs(calc_t),df=n-k-1)
#CI FOR BETA ESTIMATES
tval=abs(qt(alp/2,n-p))
LL=0
UL=0
for(nn in 1:p)
{
LL[nn]=as.matrix(betao[nn]-(tval*dia_elt[nn]),p,1)
UL[nn]=as.matrix(betao[nn]+(tval*dia_elt[nn]),p,1)
}
#-------------------------------------------------
#Scaling_Unit Length to find standardized regression estimates
dif=matrix(0,n,k)
ssu=0
w=matrix(0,n,k)
for(h in 1:k)
{
for(q in 1:n)
{
dif[q,h]=(xd[q,h]-mean(xd[,h]))^2
ssu[h]=sqrt(sum(dif[,h]))
w[q,h]=(xd[q,h]-mean(xd[,h]))/ssu[h]
}
}
y_w=0
for(q in 1:n)
{
y_w[q]=as.matrix((y[q]-mean(y))/sqrt(SST),n,1)
}
#-------------------------------------------------
#VIF COMPUTATION
mcm=solve(t(w)%*%w) #Multicollinearlty matrix
dia_elt_1=0 #To pick diagonal elements of Covariance matrix - covao
for (i in 1:k)
{
dia_elt_1[i]=mcm[i,i]
}
dia_elt_w=as.matrix(dia_elt_1,1,k)
Const=NA
VIF=as.matrix(rbind(Const,dia_elt_w),p,1)
st_beta_OLS=solve(t(w)%*%w)%*%t(w)%*%y_w
st_beta_o=rbind(Const,st_beta_OLS)
#-------------------------------------------------
#RESULTS_OLE
SE_O=as.matrix(dia_elt,p,1) #SE
ols_ans=cbind(betao,SE_O,calc_t,pval_O,LL,UL,st_beta_o,VIF)
colnames(ols_ans)=c("Beta","SE","t-stat","p-val","LL","UL","St.b's","VIF")
rownames(ols_ans)=c("Intercept","X1","X2")
ols_overall=cbind(Rsquare_O,AdjRsquare_O,MSRes,F_O,pval_F_O)
colnames(ols_overall)=c("R-2","AdjR-2","MSE","F-Stat", "p-val")
op=list(round(ols_overall,4),round(ols_ans,4))
names(op)=c("OverallSummary","Estimates")
opRounding off digits – Intermediate Steps
If an attempt is made to round the digits (whatever the number of digits) there will be serious consequences in the final answers. In some cases it may even provide an implausible values so that subsequent computation cannot be proceeded.
Linear model can be used to illustrate this idea, by rounding the inverse of model matrix. This has a serious impact in the SS – Sum of Squares calculations, which are always non-negative.
R code for that (derived from the above) is,
#Intention of this example is to show how rounding of digits in the #intermittent stage of computation leads to implausible values of sum of #squares
#-------------------------------------------------
x1=c(37,45,38,42,31)
x2=c(4,0,5,2,4)
y1=c(71,66,75,70,65) #Response variable
#-------------------------------------------------
#Necessary Inputs
k=2 #No of regressors
n=length(x1) #No of data points
p = k+1
alp = 0.05 #Alpha for (1-Alpha)100% CI
#-------------------------------------------------
#Given Predictors in matrix form
xd=matrix(c(x1,x2),n,k)
#pairs(xd)
xo=matrix(1,n,1)
x=cbind(xo,xd) #Modal matrix
m=t(x)%*%x #X'X
d=det(m) #Determinant of m to chk singularity
y=as.matrix(y1,n,1)
#-------------------------------------------------
#OLS Estimations - Det(model matrix) = 0 or not
if (d!=0)
{
print("Matrix X'X is NOT singular")}else
print("Matrix X'X is singular")
#-------------------------------------------------
#OLS Estimations - Inverse of modal matrix
#mi is inverse of X'X - WITH ROUNDING - MAY CREATE TROUBLE
mi=round(solve(m),2)
#-------------------------------------------------
#OLS Estimations - Regression beta coeffts
betao = mi%*%t(x)%*%y
#-------------------------------------------------
#Calculations based on OLS
#Residuals
#Sum of Squares_Residue
#MSReg = estimate of SigmaSq
#Covariance of reg coeffts (Betas)
SSRes=as.vector(t(y)%*%y-t(betao)%*%t(x)%*%x%*%betao)
MSRes=SSRes/(n-k-1)
covao=MSRes*mi
cbind(SSRes,MSRes)
#------------------------------------------------If we run the above part we get,
> cbind(SSRes,MSRes)
SSRes MSRes
[1,] -2109919 -1054959This is just because of the intermediate rounding the numbers. That has an adverse effect and subsequent steps cannot be implemented