#Simple Regression

sqrt(9)
sin(2)
2/5
c(1,2,4)
mdat <- matrix(c(1,2,3, 11,12,13), nrow = 2, ncol = 3, 
               dimnames = list(c("row1", "row2"),
                               c("C.1", "C.2", "C.3")),byrow=TRUE)
mdat
y=c(70,65,90,95,110,115,120,140,155,150)
x=c(80,100,120,140,160,180,200,220,240,260)
mean(x)
sum(x^2)
var(x)
cov(x,y)
cor(x,y)

model=lm(y~x)
summary(model)
anova(model)


USairpollution=matrix(c(46,11,24,47,11,31,110,23,65,26,9,17,17,35,56,10,28,14,14,13,30,10,10,16,29,18,9,31,14,69,10,61,94,26,28,12,29,56,29,8,36,47.6,56.8,61.5,55.0,47.1,55.2,50.6,54.0,49.7,51.5,66.2,51.9,49.0,49.9,49.1,68.9,52.3,68.4,54.5,61.0,55.6,61.6,75.5,45.7,43.5,59.4,68.3,59.3,51.5,54.6,70.3,50.4,50.0,57.8,51.0,56.7,51.1,55.9,57.3,56.6,54.0,44,46,368,652,391,35,3344,462,1007,266,641,454,104,1064,412,721,361,136,381,91,291,337,207,569,699,275,204,96,181,1692,213,347,343,197,137,453,379,775,434,125,80,116,244,497,905,463,71,3369,453,751,540,844,515,201,1513,158,1233,746,529,507,132,593,624,335,717,744,448,361,308,347,1950,582,520,179,299,176,716,531,622,757,277,80,8.8,8.9,9.1,9.6,12.4,6.5,10.4,7.1,10.9,8.6,10.9,9.0,11.2,10.1,9.0,10.8,9.7,8.8,10.0,8.2,8.3,9.2,9.0,11.8,10.6,7.9,8.4,10.6,10.9,9.6,6.0,9.4,10.6,7.6,8.7,8.7,9.4,9.5,9.3,12.7,9.0,33.36,7.77,48.34,41.31,36.11,40.75,34.44,39.04,34.99,37.01,35.94,12.95,30.85,30.96,43.37,48.19,38.74,54.47,37.00,48.52,43.11,49.10,59.80,29.07,25.94,46.00,56.77,44.68,30.18,39.93,7.05,36.22,42.75,42.59,15.17,20.66,38.79,35.89,38.89,30.58,40.25,135,58,115,111,166,148,122,132,155,134,78,86,103,129,127,103,121,116,99,100,123,105,128,123,137,119,113,116,98,115,36,147,125,115,89,67,164,105,111,82,114),41)
cities = c("Albany", "Albuquerque","Atlanta","Baltimore","Buffalo","Charleston","Chicago", "Cincinnati","Cleveland","Columbus","Dallas","Denver","DesMoines","Detroit","Hartford","Houston","Indianapolis","Jacksonville","Kansas City","Little Rock","Louisville","Memphis","Miami", "Milwaukee","Minneapolis","Nashville","New Orleans","Norfolk","Omaha","Philadelphia", "Phoenix","Pittsburgh","Providence","Richmond","Salt Lake City","San Francisco","Seattle","St. Louis","Washington","Wichita","Wilmington")
variables = c("SO2","temp","manu","popul","wind","precip","predays")
colnames(USairpollution) = variables
rownames(USairpollution)= cities
d<-USairpollution
d[,1]
d[1,1]
y<-d[,1]

x1<-d[,2]
x2<-d[,3]
x3<-d[,4]
data1=cbind(y,x1)

m1=lm(y~x1)
summary(m1)


plot(x1,y)
abline(a=-10,b=1)
abline(lm(y~x1),col="red")

yhat<-fitted(m1)
predict(m1)
confint(m1,level=0.90)

#regression diagnostics
lm.influence(m1)
par(mfrow=c(2,2))
plot(m1)
lev = hat(model.matrix(m1))
plot(lev)
data1[lev >0.2]
#hat : a vector containing the diagonal of the hat matrix.
#coefficients: a matrix whose i-th row contains the change in the estimated 
#coefficients which results when the i-th case is dropped from the regression. Note that aliased coefficients are not included in the matrix.
#sigma:a vector whose i-th element contains the estimate of the residual 
#standard deviation obtained when the i-th case is dropped from the regression (The approximations needed for GLMs can result in this being NaN.)
#wt.res: a vector of weighted residuals.


r<-resid(lm(y~x1))
sr=rstandard(m1)
dffits(m1)
cooks.distance(m1)
dfbeta(m1)
covratio(m1)


plot(x1,yhat)
lines(x1,yhat)
qqnorm(r)
qqline(m1$res)
sr=rstandard(m1)
hist(sr)
plot(yhat,r)

c<-predict(lm(y~x1),int="c");c

yhat<-fitted(lm(y~x1))
cor.test(x1,y,alternative = "two.sided",method = "pearson",conf.level = 0.98)
anova(lm(y~x1), test="F")

pred <- lm(y~x1)
pc <- predict(pred, interval=c("confidence"))

pp<-predict(pred,newdata=data.frame(x=x1), interval=c("prediction"))
plot(x1,yhat)
plot(x1,y,ylim=range(y,yhat, pp,pc, na.rm=T))
matlines(x1,yhat,lty=c(3,1,2), col="green")
matlines(x1,pc[,c("lwr","upr")],lty=c(1,1,2), col="black")
matlines(x1,pp[,c("lwr","upr")], col="brown")

par(mfrow=c(2,2), mex=.6)
plot(pred)
abline(v=4/41)
abline(h=2)
abline(h=-2)


par(mfrow=c(2,2), mex=0.6)
plot(rstandard(pred))
plot(rstudent(pred))
plot(dffits(pred),type="l")
matplot(dfbetas(pred),type="l", col="black")
lines(sqrt(cooks.distance(pred)), lwd=2)






####Simulation###########################
library(MASS)
require(graphics)


mu=c(1.2,1,2,3,2)
Sigma=matrix(c(.81,.4,.3,-.2,-.1,.4,2.25,.4,.3,-.2,.3,.4,1,.4,.3
,-.2,.3,.4,.64,.4,-.1,-.2,.3,.4,.49)
,5,5)

data=mvrnorm(200, mu, Sigma)
x1<-data[,1]
x2<-data[,2]
x3<-data[,3]
x4<-data[,4]
x5<-data[,5]
e=runif(200, min=0, max=2.25)
y<-5+3*x1+2*(x2)+x3-2*x4-3*x5+e
sd<-cbind(y,x1,x2,x3,x4,x5)
x=cbind(rep(1,200),x1,x2,x3,x4,x5)







cor(sd)
#multiple Regression

pairs(sd, gap=0, cex.labels=0.9)
#gap and cex.labels control the visual appearance by
#removing the space between subplots and decreasing the font size.

R1<-lm(y~x1+x2+x3+x4+x5)
summary(R1)
yhat1<-fitted(R1)
sqrt(t(y-yhat1)%*%(y-yhat1)/194)
#B1=B2=...=B5=0
anova(R1)

#Test of B0=B4=B5=0;F0=(SSE(RM)-SSE(FM))/(q*MSE(FM))
R2<-lm(y~-1+x1+x2+x3)
summary(R2)
anova(R2,R1)

#Test of B2=B5,B0=2,B4=-2
z2=x2+x5
c1=matrix(1,200,1)
w=y-2*c1+2*x4
summary(lm(w~-1+x1+z2+x3))
R3<-lm(w~-1+x1+z2+x3)
anova(R3)
F0=((2270.7-82.9)/3)/0.4;F0


vcov(lm(y~x1+x2+x3+x4+x5))
coef(lm(y~x1+x2+x3+x4+x5))
alias(lm(y~x1+x2+x3+x4+x5))

############Hale Tamrin Dr Shahkar###############
x1<-c(-1,1,-1,1,0,0,0)
x2<-c(-1,-1,1,1,0,1,2)
x3<-x1*x1
#x3<-c(1,1,1,1,0,0,0)

A<-matrix(c(1,-1,-1,1,1,1,-1,1,1,-1,1,1,1,1,1,1,1,0,0,0,1,0,1,0,1,0,2,0),4,7)
x1<-A[2,]
x2<-A[3,]
x3<-A[4,]
y<- matrix(c(1,4,8,9,3,8,9),7,1,byrow=TRUE)

S<-cbind(y,x1,x2,x3)
S

FM<-lm(y~x1+x2+x3)
FM
summary(FM)
yhat<-fitted(FM)
sqrt(t(y-yhat)%*%(y-yhat)/3)

## Azmon B1=-B2-2B3 & B0=-B2-3B3 ##

u<-x1+x2-1
u2<--2*x1+x3-3
RM<-lm(y~-1+I(-x1+x2-1)+I(-2*x1+x3-3))
summary(RM)
anova(RM)
anova(RM,FM)


#multiple Regression with matrix algebra

C=solve(t(x)%*%x)
ybar=mean(y)
bhat=C%*%t(x)%*%y
yhat=x%*%bhat
sse=t(y-yhat)%*%(y-yhat)
mse=sse/194

vhat=matrix(0,6,6)
for(i in 1:6) 
for(j in 1:6)
vhat[i,j]=mse*C[i,j]

SSR=t(yhat)%*%yhat-200*ybar^2;SSR






###CHAPTER 3###################################################################

setwd("H:/Ch 3")

responsetransformation <- read.table("responsetransformation.txt",header=TRUE)
attach(responsetransformation)

#Figure 3.25 on page 84
plot(x,y)

#Figure 3.26 on page 85
plot(x,y)
m1 <- lm(y~x)
summary(m1)
par(mfrow=c(1,2))
StanRes1 <- rstandard(m1)
absrtsr1 <- sqrt(abs(StanRes1))
plot(x,StanRes1,ylab="Standardized Residuals")
plot(x,absrtsr1,ylab="Square Root(|Standardized Residuals|)")

#Figure 3.27 on page 86
par(mfrow=c(3,2))
plot(density(y,bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="y")
rug(y)
boxplot(y,ylab="Y")
qqnorm(y, ylab = "Y")
qqline(y, lty = 2, col=2)
sj <- bw.SJ(x,lower = 0.05, upper = 100)
plot(density(x,bw=sj,kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="x")
rug(x)
boxplot(x,ylab="x")
qqnorm(x, ylab = "x")
qqline(x, lty = 2, col=2)

#Figure 3.28 on page 87
#install.packages("alr3")
install.packages("car")
library(car)


#You will be asked to 
#--- Please select a CRAN mirror for use in this session ---
#library(alr3)
#inverse.response.plot(m1,key=TRUE)

par(mfrow=c(1,1))

k=invResPlot(m1)
invResPlot(m1,key=TRUE)
# Click on the section of the plot that you wish to put the figure legend


#Figure 3.29 on page 88
k=inverseResponsePlot(m1,lam=seq(-3,3,by=0.01))
plot(k[[1]],k[[2]],pch=".",ylab=expression(RSS(lambda)),xlab=expression(lambda))
MM<-k[sort.list(k[,2]),]
MM[1,]

#Figure 3.29 on page 88 without using library alr3
lambda<-seq(-1,1,0.001)
n<-length(lambda)
sse=0
A<-matrix(0,ncol=length(y),nrow=n)
A
for(i in 1:n)
{
if(lambda[i]==0) A[i,]<-log(y)
else A[i,]<-(y^lambda[i]-1)/lambda[i]
m<-lm(y~x)
z<-fitted(m)
k<-lm(z~A[i,])
resid(k)
sse[i]<-sum(resid(k)^2)
}
sse
plot(lambda,sse,ylab='sse(lambda)',xlab='lambda',pch=".")
F<-cbind(lambda,sse)
H<-F[sort.list(F[,2]),]
H[1,]


#Figure 3.30 on page 92
library(MASS)
par(mfrow=c(1,2))
boxcox(m1,lambda=seq(0.28,0.39,length=100))
boxcox(m1,lambda=seq(0.332,0.333,length=100))

m=boxcox(m1,lambda=seq(0.28,0.39,length=100))
mm=cbind(m$x,m$y)
NN<-mm[sort.list(mm[,2]),]
NN[1,]
NN[100,]

#Regression output & Figure 3.31 on page 93
ty <- y^(1/3)
par(mfrow=c(2,2))
sj <- bw.SJ(ty,lower = 0.05, upper = 100)
plot(density(ty,bw=sj,kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab=expression(Y^(1/3)))
rug(ty)
boxplot(ty,ylab=expression(Y^(1/3)))
qqnorm(ty, ylab = expression(Y^(1/3)))
qqline(ty, lty = 2, col=2)
m2 <- lm(ty~x)
plot(x,ty,ylab=expression(Y^(1/3)))
abline(m2)
summary(m2)

detach(responsetransformation)



install.packages("alr4")
library(alr4)
data(salarygov)
attach(salarygov)

#Figure 3.32 on page 96
m1 <- lm(MaxSalary~Score)
par(mfrow=c(2,2))
plot(Score,MaxSalary)
abline(m1,lty=2,col=2)
StanRes1 <- rstandard(m1)
absrtsr1 <- sqrt(abs(StanRes1))
plot(Score,StanRes1,ylab="Standardized Residuals")
plot(Score,absrtsr1,ylab="Square Root(|Standardized Residuals|)")
abline(lsfit(Score,absrtsr1),lty=2,col=2)

#Output from R on page 96

summary(tranxy <-powerTransform(cbind(MaxSalary,Score)))

#Figure 3.33 on page 97
par(mfrow=c(3,2))
plot(density(MaxSalary,bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="MaxSalary")
rug(MaxSalary)
boxplot(MaxSalary,ylab="MaxSalary")
qqnorm(MaxSalary, ylab = "MaxSalary")
qqline(MaxSalary, lty = 2, col=2)
plot(density(Score,bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="Score")
rug(Score)
boxplot(Score,ylab="Score")
qqnorm(Score, ylab = "Score")
qqline(Score, lty = 2, col=2)

#Figure 3.34 on page 97
par(mfrow=c(1,1))
plot(sqrt(Score),log(MaxSalary),xlab=expression(sqrt(Score)))
abline(lsfit(sqrt(Score),log(MaxSalary)),lty=2,col=2)

#Figure 3.35 on page 98
par(mfrow=c(3,2))
plot(density(log(MaxSalary),bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab="log(MaxSalary)")
rug(log(MaxSalary))
boxplot(log(MaxSalary),ylab="log(MaxSalary)")
qqnorm(log(MaxSalary), ylab = "log(MaxSalary)")
qqline(log(MaxSalary), lty = 2, col=2)
sj <- bw.SJ(sqrt(Score),lower = 0.05, upper = 100)
plot(density(sqrt(Score),bw=sj,kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab=expression(sqrt(Score)))
rug(sqrt(Score))
boxplot(sqrt(Score),ylab=expression(sqrt(Score)))
qqnorm(sqrt(Score), ylab=expression(sqrt(Score)))
qqline(sqrt(Score), lty = 2, col=2)

#Figure 3.36 on page 99
m2 <- lm(log(MaxSalary)~sqrt(Score))
par(mfrow=c(1,2))
StanRes2 <- rstandard(m2)
absrtsr2 <- sqrt(abs(StanRes2))
plot(sqrt(Score),StanRes2,ylab="Standardized Residuals",xlab=expression(sqrt(Score)))
plot(sqrt(Score),absrtsr2,ylab="Square Root(|Standardized Residuals|)",xlab=expression(sqrt(Score)))
abline(lsfit(sqrt(Score),absrtsr2),lty=2,col=2)

#R output on page 99

summary(tranx <- powerTransform(Score))

#Figure 3.37 on page 100
m3 <- lm(MaxSalary~sqrt(Score))
par(mfrow=c(1,1))
inverseResponsePlot(m3,key=TRUE)
#Click on the plot where you want to put the legend

#Figure 3.38 on page 101
par(mfrow=c(2,2))
plot(density(MaxSalary^-0.25,bw="SJ",kern="gaussian"),type="l",
main="Gaussian kernel density estimate",xlab=expression(MaxSalary^-0.25))
rug(MaxSalary^-0.25)
boxplot(MaxSalary^-0.25,ylab=expression(MaxSalary^-0.25))
qqnorm(MaxSalary^-0.25,ylab=expression(MaxSalary^-0.25))
qqline(MaxSalary^-0.25,lty=2, col=2)

#Figure 3.39 on page 102
par(mfrow=c(2,2))
plot(sqrt(Score),MaxSalary^-0.25,xlab=expression(sqrt(Score)),ylab=expression(MaxSalary^-0.25))
abline(lsfit(sqrt(Score),MaxSalary^-0.25),lty=2,col=2)
m3 <- lm(MaxSalary^-0.25~sqrt(Score))
StanRes3 <- rstandard(m3)
absrtsr3 <- sqrt(abs(StanRes3))
plot(sqrt(Score),StanRes3,ylab="Standardized Residuals",xlab=expression(sqrt(Score)))
plot(sqrt(Score),absrtsr3,ylab="Square Root(|Standardized Residuals|)",xlab=expression(sqrt(Score)))
abline(lsfit(sqrt(Score),absrtsr3),lty=2,col=2)
summary(m2);summary(m3)
detach(salarygov)


#CHAPTER 4###################################################################


cleaningwtd <- read.table("cleaningwtd.txt",header=TRUE)
attach(cleaningwtd)

#Regression output on page 117
wm1 <- lm(Rooms~Crews,weights=1/(round(StdDev,4))^2)
summary(wm1)
predict(wm1,newdata=data.frame(Crews=c(4,16)),interval="prediction",level=0.95)
predict(wm1,newdata=data.frame(Crews=c(4,16)),interval="confidence",level=0.95)



#Regression output on page 120
ynew <- Rooms/StdDev
x1new <- 1/StdDev
x2new <- Crews/StdDev
wm2 <- lm(ynew~x1new + x2new - 1)
summary(wm2)
predict(wm2,newdata=data.frame(x1new=c(1/4.966555,1/12.000463),x2new=c(4/4.966555,16/12.000463)),interval="prediction",level=0.95)

detach(cleaningwtd)

CHAPTER 6###################################################################
setwd("C:/Users/Mahdi/Desktop/Sheather/Data")
magazines <- read.csv("magazines.csv", header=TRUE)
attach(magazines)

#R output on page 177
library(car)
summary(powerTransform(cbind(AdPages,SubRevenue,NewsRevenue)~1))


#Figure 6.21 on page 178
pairs(AdRevenue~AdPages+SubRevenue+NewsRevenue)

#Figure 6.22 on page 179
tAdPages<- log(AdPages)
tSubRevenue <- log(SubRevenue)
tNewsRevenue <- log(NewsRevenue)
m1 <- lm(AdRevenue~log(AdPages)+log(SubRevenue)+log(NewsRevenue))
par(mfrow=c(1,1))
inverseResponsePlot(m1,key=TRUE)

#R output on page 179
summary(tranxy <-powerTransform(cbind(AdRevenue,AdPages,SubRevenue,NewsRevenue)~1))

#Figure 6.23 on page 180
pairs(log(AdRevenue)~log(AdPages)+log(SubRevenue)+log(NewsRevenue))

#Figure 6.24 on page 181
m3 <- lm((AdRevenue)~(AdPages)+(SubRevenue)+(NewsRevenue))
m2 <- lm(log(AdRevenue)~log(AdPages)+log(SubRevenue)+log(NewsRevenue))
m1 <- lm((AdRevenue^0.23)~log(AdPages)+log(SubRevenue)+log(NewsRevenue))
summary(m1);summary(m2);summary(m3)
par(mfrow=c(2,2))
StanRes2 <- rstandard(m2)
plot(log(AdPages),StanRes2,ylab="Standardized Residuals")
plot(log(SubRevenue),StanRes2,ylab="Standardized Residuals")
plot(log(NewsRevenue),StanRes2,ylab="Standardized Residuals")
plot(m2$fitted.values,StanRes2,ylab="Standardized Residuals",xlab="Fitted Values")

#Figure 6.25 on page 181
par(mfrow=c(1,1))
plot(m2$fitted.values,log(AdRevenue),xlab="Fitted Values")
abline(lsfit(m2$fitted.values,log(AdRevenue)))

#Figure 6.26 on page 182
par(mfrow=c(2,2))
plot(m2)
abline(v=2*4/204,lty=2)
abline(h=c(2,-2),lty=3)


#Figure 6.27 on page 183
library(car)
par(mfrow=c(2,2))
avPlot(m2,variable=log(AdPages),ask=FALSE,identify.points=FALSE)
avPlot(m2,variable=log(SubRevenue),ask=FALSE,identify.points=FALSE)
avPlot(m2,variable=log(NewsRevenue),ask=FALSE,identify.points=FALSE)

#Regression output on page 183
summary(m2)

detach(magazines)



bridge <- read.table("bridge.txt", header=TRUE)
attach(bridge)

summary(tranxy <- powerTransform(cbind(Time,DArea,CCost,Dwgs,Length,Spans)~1))
testTransform(tranxy, c(-0.25, 0,-0.25, -0.25,-0.25,-0.5))

# fit linear model with transformed response:
coef(tranxy, round=TRUE)
summary(lm(bcPower(Time,tranxy$lam)~bcPower(DArea,tranxy$lam)+bcPower(CCost,tranxy$lam)+bcPower(Dwgs,tranxy$lam)+bcPower(Length,tranxy$lam)+bcPower(Spans,tranxy$lam)))

# save the rounded transformed values, plot them with a separate
# color for each highway type
transformedY <- bcPower(with(bridge,cbind(Time,DArea,CCost,Dwgs,Length,Spans)),
                coef(tranxy))
pairs(transformedY, col=as.numeric(bridge$Time)) 

#Figure 6.40 page 198
pairs((Time)~(DArea)+(CCost)+(Dwgs)+(Length)+(Spans),data=bridge)
pairs(log(Time)~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans),data=bridge)

#Figure 6.41 page 199
m1 <- lm(log(Time)~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans))
StanRes1 <- rstandard(m1)
par(mfrow=c(2,3))
plot(log(DArea),StanRes1, ylab="Standardized Residuals")
plot(log(CCost),StanRes1, ylab="Standardized Residuals")
plot(log(Dwgs),StanRes1, ylab="Standardized Residuals")
plot(log(Length),StanRes1, ylab="Standardized Residuals")
plot(log(Spans),StanRes1, ylab="Standardized Residuals")
plot(m1$fitted.values,StanRes1, ylab="Standardized Residuals",xlab="Fitted values")

#Figure 6.42 page 199
par(mfrow=c(1,1))
plot(m1$fitted.values,log(Time),xlab="Fitted Values")
abline(lsfit(m1$fitted.values,log(Time)))

#Figure 6.43 page 200
par(mfrow=c(2,2))
plot(m1)
abline(v=2*6/45,lty=2)
abline(h=c(2,-2),lty=3)

#Regression output on page 200
summary(m1)

#Figure 6.44 page 201
library(alr3)
mmps(m1,layout=c(2,3))

#R output for detecting the mulicollinearity
logDArea <- log(DArea)
logCCost <- log(CCost)
logDwgs <- log(Dwgs)
logLength <- log(Length)
logSpans <- log(Spans)
X <- cbind(rep(1,45),logDArea,logCCost,logDwgs,logLength,logSpans)
c <- cor(X)
round(c,3)

X <- cbind(logDArea,logCCost,logDwgs,logLength,logSpans)
c <- cor(X)
round(c,3)

#Figure 6.45 on page 202
library(car)
par(mfrow=c(2,3))
avPlot(m1,variable=log(DArea),ask=FALSE,identify.points=FALSE)
avPlot(m1,variable=log(CCost),ask=FALSE,identify.points=FALSE)
avPlot(m1,variable=log(Dwgs),ask=FALSE,identify.points=FALSE)
avPlot(m1,variable=log(Length),ask=FALSE,identify.points=FALSE)
avPlot(m1,variable=log(Spans),ask=FALSE,identify.points=FALSE)

#R output on page 203
library(car)
vif(m1)

detach(bridge)

CHAPTER 7###################################################################
setwd("C:/Users/Mahdi/Desktop/Sheather/Data")
#install.packages("olsrr")
library(olsrr)

m1 <- lm(log(Time)~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans))
ols_step_all_possible(m1)
logDArea <- log(DArea)
logCCost <- log(CCost)
logDwgs <- log(Dwgs)
logLength <- log(Length)
logSpans <- log(Spans)
X <- cbind(logDArea,logCCost,logDwgs,logLength,logSpans)
#install.packages("leaps")
library(leaps)
b <- regsubsets(as.matrix(X),log(Time))
rs <- summary(b)
par(mfrow=c(1,2))
plot(1:5,rs$adjr2,xlab="Subset Size",ylab="Adjusted R-squared")
library(car)
subsets(b,statistic=c("adjr2"))

#Table 7.1 on page 235
#Calculate adjusted R-squared
rs$adjr2
om1 <- lm(log(Time)~log(Dwgs))
om2 <- lm(log(Time)~log(Dwgs)+log(Spans))
om3 <- lm(log(Time)~log(Dwgs)+log(Spans)+log(CCost))
om4 <- lm(log(Time)~log(Dwgs)+log(Spans)+log(CCost)+log(DArea))
om5 <- m1
#Subset size=1
n <- length(om1$residuals)
npar <- length(om1$coefficients) +1
#Calculate AIC
extractAIC(om1,k=2)
#Calculate AICc
extractAIC(om1,k=2)+2*npar*(npar+1)/(n-npar-1)
#Calculate BIC
extractAIC(om1,k=log(n))
#Subset size=2
npar <- length(om2$coefficients) +1
#Calculate AIC
extractAIC(om2,k=2)
#Calculate AICc
extractAIC(om2,k=2)+2*npar*(npar+1)/(n-npar-1)
#Calculate BIC
extractAIC(om2,k=log(n))
#Subset size=3
npar <- length(om3$coefficients) +1
#Calculate AIC
extractAIC(om3,k=2)
#Calculate AICc
extractAIC(om3,k=2)+2*npar*(npar+1)/(n-npar-1)
#Calculate BIC
extractAIC(om3,k=log(n))
#Subset size=4
npar <- length(om4$coefficients) +1
#Calculate AIC
extractAIC(om4,k=2)
#Calculate AICc
extractAIC(om4,k=2)+2*npar*(npar+1)/(n-npar-1)
#Calculate BIC
extractAIC(om4,k=log(n))
#Subset size=5
npar <- length(om5$coefficients) +1
#Calculate AIC
extractAIC(om5,k=2)
#Calculate AICc
extractAIC(om5,k=2)+2*npar*(npar+1)/(n-npar-1)
#Calculate BIC
extractAIC(om5,k=log(n))

#Regression output on pages 235 and 236
summary(om2)
summary(om3)

#Output from R on page 237
backAIC <- step(m1,direction="backward", data=bridge)
backBIC <- step(m1,direction="backward", data=bridge, k=log(n))

#Output from R on page 238
mint <- lm(log(Time)~1,data=bridge)
forwardAIC <- step(mint,scope=list(lower=~1, 
upper=~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans)),
direction="forward", data=bridge)
forwardBIC <- step(mint,scope=list(lower=~1, 
upper=~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans)),
direction="forward", data=bridge,k=log(n))

StepwiseReg <- step(mint,scope=list(lower=~1, 
upper=~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans)),
direction="both", data=bridge)

detach(bridge)


PROJECT#######################################

library(sn)
library(MASS)
n=200;
mu=c(1.2,1,2,3,2,5,2)
Sigma=cor(longley)
alpha <- c(2,6,1,4,8,2,3)
X <- rmsn(n,  mu, Sigma, alpha)
X1=X[,1];X2=X[,2];X3=X[,3];X4=X[,4];X5=X[,5];X6=X[,6];X7=X[,7];
Z1=rnorm(n,5,1)
Z2=rnorm(n,4,1)
X<-cbind(rep(1,n),Z1,Z2,X1,X2,X3,X4,X5,X6,X7)
b=matrix(c(1.5,1,3,0,2,3,-1,4,2,-3),10,1)
e=mvrnorm(1, rep(0,n), diag(n))

Y=(1.5+log(Z1)+3*sqrt(Z2)+2*log(X[,2])+3*log(X[,3])-log(X[,4])+4*sqrt(X[,5])+2*(X[,6])^2
-3*X[,7]+e)^(1)
data=cbind(Y,X);data





mm=lm(Y~X[,2]+X[,3]+X[,4]+X[,5]+X[,6]+X[,7])
summary(mm)