Endale B Altaye
Bowling Green State University
(419) 378 - 9131
ebaltay@bgsu.edu
Bowling Green State University
(419) 378 - 9131
ebaltay@bgsu.edu
Multivariate multiple linear regression: several y’s and several x’s. Multivariate refers to the dependent variables and multiple pertains to the independent variables.
Steps:
* Overall significance test. Testing the existence of at least one predictor variable to explain at least one of the
response variable.
* Variable selection using Backward elimination
* Model improvement
+ Interaction included
+ R-square and S used for assessment
* Multivariate Normality assessed
* Residual assessment
libs <- c("car","xlsx","psych","dplyr","ggplot2","knitr","gridExtra")
lapply(libs,require, character.only= TRUE)
library(repr)
options(repr.plot.width=7,repr.plot.height=4)
multdata=read.xlsx("C:/Users/selamina/...data.xlsx",1)
#message("H0: B = O vs H1: B != O, (We want to know whether at least one Bjk != 0,
#to explain at least one of the response variable.")
multdataModel="cbind(y1,y2,y3)~x3+x4+x5+x6+x7+x8+x9+x10+x11+x12+x13+x14+x15+x16"
multdataFullModel<- lm(multdataModel, data=multdata)
c=diag(15)[2:15,] # first column is for the const term
(OverallFit <- linearHypothesis(model = multdataFullModel, hypothesis.matrix = c)) #Spits H, E and Sign for over all
#message("From all the four multivariate test statistics, we can conclude as there exist at least one potential predictor variable
#to explain at least one of the response variable.")
As we can see in all four multivariate test statistics, we can conclude that there exist at least one potential predictor variable to explain at least one of the response variable at 5 percent level of significance.
I chose to implement the Backward elimination, but could be tried Stepwise or forward selection as well. The concept for backward elimination is to compute $X_j$ given $X_{-j}$ , using Sum of Square of Type III; then drop the one with most insignificant one. Note that $X_{-j}$ means all the predictor variables except $j^{th}$.
#message("The concept is to compute Xj-given-all_X, using Sum of Square of Type III;
#then drop the one with most insignificant one or the least contributer.")
Manova(multdataFullModel,test.statistic="Wilks")
#message("The P-value per variable has been assessed and the one with largest p-value is dropped, which is X10 with p-value=0.89213")
The p-value for each variable has been assessed and the one with largest p-value is dropped, which is X10 with p-value=0.89213. Then X10 is dropped to move to next step.
reduced1=update(multdataFullModel,.~.-x10)
Manova(reduced1,test.statistic="Wilks")
Steps deleted here
reduced6=update(reduced5,.~.-x4)
Manova(reduced6,test.statistic="Wilks")
reduced7=update(reduced6,.~.-x11)
Manova(reduced7,test.statistic="Wilks")
#summary(Manova(reduced7))
Since the biggest p-value is 0.26889, its acceptable not to drop any variable afterwards, just to be on the safe side.
#1. R-Square For the final model selected using BackWard elimination
message("#1.R-Square")
Intercept=c(rep(1,35))
matX=multdata %>% select(x5 , x6, x9 , x12, x13, x14 , x16)
matX$x16=ifelse(matX$x16=="Sec",1,0) #Sec=0, else 1
matX=as.matrix(cbind(Intercept,matX))
yY=multdata %>% select(y1,y2,y3)
meanY=as.matrix(apply(yY,2,mean))
yY=as.matrix(yY)
bhat=solve(t(matX)%*%matX)%*%t(matX)%*%yY # Pf_H. behaves strange
rSquare=solve(t(yY)%*%yY-35*meanY%*%t(meanY))%*%(t(bhat)%*%t(matX)%*%yY-35*meanY%*%t(meanY))
tr(rSquare)/length(meanY)
#rSquare
#2.S
message("S-Pooled")
Sp=t(yY)%*%yY-t(bhat)%*%t(matX)%*%yY
Sp
message("trace of S-Pooled")
(tr(Sp))
The backward elimination algorithm returns a model of seven predictor variables. According to this model, 47.43% of the variability in the response variables could be explained by the variables in the model.
Picking the most significant once in the selected model above and try to improve model by including interaction among these variables which is between $x_5$,$x_{14}$ and $x_{16}$.
modelInteraction=lm(cbind(y1,y2,y3)~x5+x6+x9+x12+x13+x14+x16,data=multdata)
#Manova(modelInteraction,statistics="Wilks")
#message("Picked the most significant once and try to improve model by including interaction among these x5,x14 ,x16")
modelInteraction1=lm(cbind(y1,y2,y3)~x5+x6+x9+x12+x13+x14+x16+x5*x14+x5*x16+x14*x16,data=multdata)
Manova(modelInteraction1,test.statistics="Wilks")
anova(modelInteraction, modelInteraction1, test = "Wilks") #
At 20% level of significance, it looks significant to include interaction to the model. However, the model could be improved by droping the most insignificant interactions.
modelInteraction2=lm(cbind(y1,y2,y3)~x5+x6+x9+x12+x13+x14+x16+x5*x16+x14*x16,data=multdata)
#Manova(modelInteraction2,statistics="Wilks")
message("## R-Square for updated model")
Intercept=c(rep(1,35))
matX=multdata %>% select(x5 , x6, x9 ,x12, x13, x14 , x16) %>% mutate(x16=ifelse(x16=="Sec",1,0),x5x16=x5*x16,x14x16=x14*x16)
matX=as.matrix(cbind(Intercept,matX))
yY=multdata %>% select(y1,y2,y3)
meanY=as.matrix(apply(yY,2,mean))
yY=as.matrix(yY)
bhat=solve(t(matX)%*%matX)%*%t(matX)%*%yY # Pf_H. behaves strange
rSquare=solve(t(yY)%*%yY-35*meanY%*%t(meanY))%*% (t(bhat)%*%t(matX)%*%yY-35*meanY%*%t(meanY))
tr(rSquare)/length(meanY)
#2. S
message("## 2. S")
Sp=t(yY)%*%yY-t(bhat)%*%t(matX)%*%yY
message("S-Pooled")
Sp
message("Trace of S-Pooled")
(tr(Sp))
After including interaction, we manage R-square to go from 47.43% to 54.62% and S to go from 76.2537 to 61.17 which is good improvement to our preselected model.
Manova(modelInteraction2)
message("Final model Coefficients:")
kable(modelInteraction2$coefficients, digits=4)
The predictor variables $x_5$, $x_9$,$x_{13}$ and interaction $x_5:x_{16}$ and $x_{14}:x_{16}$ have a positive linear relationship with $Y_1$ which means that increasing one of these variables will increase the value of $y_1$ and vies versa. Whereas, the other variables have a negative linear relationship with $y_1$. Note also that $x_5$ and $x_5:x_{16}$ have a relationship close to zero towards $y_1$. It could be interpreted in a similar fashion for the other response variables.
normmultdata=data.frame(multdata %>% select( y1,y2,y3,x5 , x6, x9 ,x12, x13, x14 , x16))
normmultdata$x16=ifelse(normmultdata$x16=="Sec",1,0) #Sec=0, else 1
normmultdata=as.matrix(normmultdata)
Di_square=c(rep(0,35))
meanNmultdata=apply(normmultdata,2,mean)
CovNmultdata=cov(normmultdata)
for(i in 1:35)
{
Di_square[i]=t(normmultdata[i,]-meanNmultdata)%*%solve(CovNmultdata)%*%(normmultdata[i,]-meanNmultdata)
}
Di_square_sort=sort(Di_square)
p=ncol(normmultdata)
n=nrow(normmultdata)
alpha1=(p-2)/(2*p)
beta1=(n-p-3)/(2*(n-p-1))
Vi=(c(1:35)-alpha1)/(35-alpha1-beta1+1)
Normdata=cbind.data.frame(Ui=Di_square_sort,Vi=Vi)
ggplot(Normdata, aes(x=Normdata[,1],y=Normdata[,2]))+geom_point(stat="identity")+labs(x="U(i)",y="Vi (Quantiles of Beta)")
maxDi=round(max(Di_square_sort),2)
tabvalue=paste("at Alpha=0.05 and p=10,n=35 we get D-value=15.75") #Table A.6.
multdataResid=data.frame(modelInteraction2$residuals)
p1=ggplot(multdataResid, aes(x=c(1:35),y=y1)) + geom_point()+labs(x=" ",y="y1 Residual")
p2=ggplot(multdataResid, aes(x=c(1:35),y=y2)) + geom_point()+labs(x=" ",y="y2 Residual")
p3=ggplot(multdataResid, aes(x=c(1:35),y=y3)) + geom_point()+labs(x=" ",y="y3 Residual")
grid.arrange(p1,p2,p3,ncol=2)
According to the Normality plot and the Residuals plot, we can conclude that the selected model is a very legit one.
[Alvin_C._Rencher] Methods of Multivariate Analysis, 2nd edition.