7.1.1 Coplots

coplot(L2PhonemicAwareness~Grade1L1Reading|NonVerbalReasoning, panel= 
function(x,y,...) panel.smooth(x,y,span=.8,...),data=lafrance) 

coplot(L2PhonemicAwareness~Grade1L1Reading|NonVerbalReasoning*NamingSpeed, panel= function(x,y,...) 
panel.smooth(x,y,span=.8,...),data=lafrance) 

7.1.2 3D Graphs

scatter3d(lafrance$KinderL2Reading, lafrance$Grade1L1Reading,    
lafrance$L2PhonemicAwareness, fit=c("linear","smooth"), bg="white",
axis.scales=TRUE, grid=TRUE, ellipsoid=FALSE, xlab="KinderL2Reading",   
ylab="Grade1L1Reading", zlab="L2PhonemicAwareness") 

7.1.3 Tree Models 

library(tree) 
model=tree(lafrance5) 
plot(model) 
text(model) 

model=tree(Grade1L1ReadingPerformance~., lafrance.tree) 

7.3 Doing the Same Type of Regression as SPSS

model=lm(Grade1L1ReadingPerformance~NonverbalReasoning+KinderL2Reading 
Performance+NamingSpeed+WorkingMemory+PhonologicalAwarenessInL2,
data=lafrance5) 

summary(model)

confint(model)

model.standard=lm(scale(Grade1L1ReadingPerformance)~scale(Nonverbal 
Reasoning)+scale(KinderL2Reading Performance)+ scale(NamingSpeed)+
scale(WorkingMemory)+scale(PhonologicalAwarenessInL2), data=lafrance5) 

library(relaimpo) 
calc.relimp(model) 

RegModel.1 <- lm(Grade1L1ReadingPerformance~KinderL2ReadingPerformance+NamingSpeed+
NonverbalReasoning+PhonologicalAwarenessInL2+WorkingMemory,
data=lafrance5) 

7.5 Finding the Best Fit 

model = y~A + B + C + A:B + B:C + A:B:C  
  
model = y~A*B*C = A + B + C + A:B + B:C + A:B:C  
 
model = y~(A + B + C)^2 = A + B + C + A:B + A:C  
 
model = y~A*B*C – A:B:C

model = y~A + I(A^2) + B + I(B^2) 

7.5.1 First Steps to Finding the Minimal Adequate Model in R

library(mice) 
imp<-mice(lafrance5) 
complete(imp) 
lafrance<-complete(imp)  

lafrance$PAL2<-lafrance$PhonologicalAwarenessInL2 
lafrance$KL2WR<-lafrance$KinderL2ReadingPerformance 
lafrance$NS<-lafrance$NamingSpeed 
lafrance$G1L1WR<-lafrance$Grade1L1ReadingPerformance
 

model=lm(G1L1WR ~ ., data=lafrance, na.action=na.exclude) 

summary(model1)

model2=update(model1,~.- PAL2:KL2WR:NS, data=lafrance) 
summary(model2) 

anova(model1,model2)

AIC(model1,model2)

model3=update(model2,~.- PAL2:NS, data=lafrance)
summary(model3) 

anova(model2,model3)

model4=update(model3,~.- KL2WR:NS, data=lafrance)
summary(model4) 

anova(model3,model4)

model5=update(model4,~.- PAL2:KL2WR, data=lafrance) 
summary(model5) 

anova(model4,model5)

model6=update(model5,~.-NS, data=lafrance)
summary(model6)

anova(model5,model6)

model7=update(model6,~.-KL2WR, data=lafrance)
summary(model7) 

model8=lm(G1L1WR~1,data=lafrance,na.action=na.exclude) 
anova(model7,model8) 

library(bootStepAIC)
boot.stepAIC(model1,data=lafrance) 

7.5.2 Reporting the Results of Regression in R

deviance(model1) 
[1] 3.484850 

7.6 Further Steps to Finding the Best Fit: Overparameterization and Polynomial Regression 

lafrance$NR<-lafrance$NonverbalReasoning 
lafrance$WM<-lafrance$WorkingMemory  

model1=lm(G1L1WR~NR+I(NR^2)+WM+I(WM^2)+NS+I(NS^2)+PAL2+
I(PAL2^2)+KL2WR+I(KL2WR^2),data=lafrance)  

summary(model1)

library(bootStepAIC)
boot.stepAIC(model1, data=lafrance)  

model2=lm(G1L1WR~ NS + I(NS^2) + PAL2 + I(PAL2^2), data=lafrance)
summary(model2) 

getfullmodel=lm(G1L1WR~NR*WM*NS*PAL2*KL2WR,data=lafrance) 
summary(getfullmodel)
#Output not printed but gives me what I need for following command
interactions=c("NR:WM", "NR:NS", "WM:NS", "NR:PAL2", "WM:PAL2", "NS:PAL2",
"NR:KL2WR", "WM:KL2WR", "NS:KL2WR", "PAL2:KL2WR") 

model3= lm(G1L1WR~NR+WM+ NS+ PAL2+ KL2WR+
PAL2:KL2WR + NR:NS + NS:PAL2 + NR:WM+NS:KL2WR, data=lafrance)
model4= lm(G1L1WR~NR+WM+ NS+ PAL2+ KL2WR+ WM:KL2WR + NR:KL2WR+
WM:PAL2+WM:NS+NR:PAL2, data=lafrance)

PAL2:KL2WR
NR:WM
WM:KL2W 
NR:KL2WR 
NR:PAL2

model5= lm(G1L1WR~NR+WM+ NS+ PAL2+ KL2WR+
PAL2:KL2WR+ NR:WM +WM:KL2WR +NR:KL2WR +NR:PAL2, data=lafrance) 

boot.stepAIC(model5, data=lafrance) 

model6=update(model5,~.-NS-PAL2:KL2WR-NR:WM, data=lafrance) 

interactions=c("NR:WM:NS", "NR:WM: PAL2", "NR:NS: PAL2", "WM:NS: PAL2", 
"NR:WM: KL2WR", "NR:NS: KL2WR", "WM:NS: KL2WR", "NR: PAL2: KL2WR", 
"WM: PAL2: KL2WR", "NS: PAL2: KL2WR") 
model7= lm(G1L1WR~NR+WM+ NS+ PAL2+ KL2WR+ 
NS:PAL2:KL2WR+NR:WM:PAL2+NR:NS:KL2WR+NR:WM:NS+ 
WM:PAL2:KL2WR, data=lafrance) 
model8= lm(G1L1WR~NR+WM+ NS+ PAL2+ KL2WR+ 
NR:PAL2:KL2WR+WM:NS:KL2WR+NR:WM:KL2WR+NR:NS:PAL2+ 
WM:NS:PAL2, data=lafrance) 

model9= lm(G1L1WR~NR+WM+ NS+ PAL2+ KL2WR+ 
NR:WM:NS:PAL2 + NR:WM:NS:KL2WR + NR:WM:PAL2:KL2WR + 
NR:NS:PAL2:KL2WR + WM:NS:PAL2:KL2WR + NR:WM:NS:PAL2:KL2WR, 
data=lafrance)

model10=lm(G1L1WR~ NR + WM + NS + I(NS^2) + PAL2 + I(PAL2^2) + KL2WR+ 
WM:KL2WR +NR:KL2WR +NR:PAL2+ #two-way interactions 
NR:WM:NS +WM:PAL2:KL2WR +  #three-way interactions NR:WM:NS:KL2WR + WM:NS:PAL2:KL2WR, #four-way interactions 
data=lafrance)  

model11=lm(G1L1WR~ NR + WM + NS+ I(NS^2)+PAL2 +KL2WR+ 
WM:KL2WR + NR:PAL2+NR:WM:NS, data=lafrance) 

model12=lm(G1L1WR~PAL2, data=lafrance) 

>anova(model10,model11) 

7.7 Examining Regression Assumptions

model12=lm(G1L1WR~PAL2,data=lafrance) plot(model12,cex=1.5) 

model13=update(model12, subset=(lafrance !=22)) 

plot(model12,which=4)

vif(model)

model.vif=lm(G1L1WR~NR+KL2WR+NS+WM+PAL2,data=lafrance)
vif(model.vif) 

library(MASS)
plot(row.names(lafrance),stdres(model12),cex=1.5) 


plot(model1)

plot(model1,which=4)

plot(model1,which=6)
 
vif(model1)

library(MASS) 
plot(row.names(lafrance),stdres(model1))

7.9.1 Visualizing Robust Reggression

library(robustbase) 
covPlot(cbind(Lafrance3$G1L1WR,Lafrance3$PAL2),which="tolEllipsePlot", 
classic=T, cor=T) 

library(robust) 
lafrance1.fit=fit.models(  
list(Robust="covRob",  
Classical="ccov"),  
data=Lafrance3)

Lafrance3<-na.omit(Lafrance3) 

plot(lafrance1.fit)

covPlot(cbind(lafrance$G1L1WR,lafrance$PAL2),which="tolEllipsePlot", 
classic=T) 

library(robust) 
your.fit=fit.models(  
list(Robust="covRob",  
Classical="ccov"),  
data=yourdata) 
plot(your.fit)

7.9.2 Robust  Methods

library(robust) lafrance3.fit=fit.models(   
list(Robust="lmRob", LS="lm"), 
formula=G1L1WR~ PAL2*KL2WR*NS,   
data=Lafrance3) 
summary(lafrance3.fit)  

lafrance3.fit=fit.models(  
list(Robust="lmRob", LS="lm"), 
formula=G1L1WR~PAL2*KL2WR*NS, 
data=lafrance3)

m1.robust=lmRob(G1L1WR~PAL2*KL2WR*NS,data=Lafrance3) summary(m1.robust) 
m2.robust=lmRob(G1L1WR~PAL2+KL2WR+NS+PAL2:KL2WR+PAL2:NS+KL2WR: 
NS,data=Lafrance3) #this deletes the 3-way interaction term PAL2:KL2WR:NS 
summary(m2.robust) 

m3.robust=lmRob(G1L1WR~PAL2+KL2WR+NS,data=Lafrance3) 
step.lmRob(m3.robust) 

m4.robust=lmRob(G1L1WR~PAL2+NS,data=Lafrance3)
summary(m4.robust) 

plot.lmRob(m1.robust) 

plot(lafrance3.fit)

my.control=lmRob.control(initial.alg="Auto", final.alg="MM", efficiency=.95) 

m1.robust=lmRob(G1L1WR~PAL2*KL2WR*NS,data=Lafrance3) 
summary(m1.robust) 
m2.robust= lmRob(G1L1WR~PAL2+KL2WR+NS- PAL2:KL2WR:NS, 
data=Lafrance3)  
anova(m1.robust,m2.robust) 
plot(m1.robust) 

#for diagnostic plots