#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) 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)) inverse.response.plot(m1,key=TRUE) #R output on page 179 library(alr3) 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) library(alr3) 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 on page 202 logDArea <- log(DArea) logCCost <- log(CCost) logDwgs <- log(Dwgs) logLength <- log(Length) logSpans <- log(Spans) 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") m1 <- lm(log(Time)~log(DArea)+log(CCost)+log(Dwgs)+log(Length)+log(Spans)) 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)