x1 <- rnorm(50)
x2 <- rt(50, 2)
# add boxplots to scatterplots
par(fig = c(0, 0.9, 0, 0.9))
plot(x1, x2)
par(fig = c(0, 0.9, 0.4, .9), new = T)
boxplot(x1, horizontal = TRUE, axes = F)
par(fig = c(0.7, 1, 0, 0.9), new = T)
boxplot(x2, axes = F)

par(fig = c(0, .5, 0, 0.5), new = T)
hist(x1, axes = F,main="",col=2)

###############################################
## histogram
x<-rnorm(50,10)
hist(x)
?hist
a<-hist(x)
a
hist(x,breaks=c(7,9,11,13))
a<-hist(x,nclass=3)
a

?persp
##############################################
#### chapter 8 ANOVA Experimental Design  ####
##############################################
### Design Of Experiments ## DOE 

### see example1 page 179 

#Example1: PlantGrowth  
data()
data(PlantGrowth)
#?PlantGrowth

PlantGrowth

PlantGrowth$group

is.factor(PlantGrowth$group)

boxplot(PlantGrowth$weight)


boxplot(weight ~ group, data = PlantGrowth)

boxplot(weight ~ group, data = PlantGrowth, main = "PlantGrowth data", ylab = "Dried weight of plants", col = "lightgray", notch = FALSE)

M<-lm(weight ~ group, data = PlantGrowth)
M
summary(M)
anova(M)

aggregate(weight ~ group,data = PlantGrowth,FUN=mean)

plot.design(weight ~ group, data = PlantGrowth)


?aov

M1<-aov(weight ~ group, data = PlantGrowth)
M1
summary(M1)
anova(M1)

model.tables(M1)


library(ggplot2)

qplot(group, weight, data = PlantGrowth,geom="point")

qplot(group, weight, data = PlantGrowth,geom="boxplot")

qplot(group, weight, data = PlantGrowth,geom="line")

qplot(weight, data = PlantGrowth,geom="density",color= group,fill= group,alpha=.3)

qplot(weight, data = PlantGrowth,geom="histogram")

#########################################
### paired comparison ### 
#########################################

## Tukey Honest Significant Differences
?TukeyHSD

TukeyHSD(M1)

## package multcomp
# install.packages("multcomp",dep=T)
library(multcomp)

?glht

Tuky<-glht(M1, linfct = mcp(group = "Tukey"))
summary(Tuky)

### alternatively, define contrast matrix directly
L <- rbind("trt1 - ctrl" = c(-1, 1, 0), 
           "trt2 - ctrl" = c(-1, 0, 1), 
           "trt1 + trt2 - 2*ctrl" = c(-2, 1, 1),
           "trt1+ctrl-2*trt2"=c(1,1,-2))
L

Mt<-glht(M1, linfct = mcp(group = L))
Mt
summary(Mt)
####################################

?pairwise.t.test

pairwise.t.test(PlantGrowth$weight , PlantGrowth$group, p.adj = "bonf")

## package agricolae
#install.packages("agricolae",dep=T)
library(agricolae)
?LSD.test
out<-LSD.test(M1,"group")

names(out)
out$statistics
out$means
out$groups
plot(out)

?SNK.test
out<-SNK.test(M1,"group")
plot(out)
out$snk
out$groups

### check residuals

#Test of  normality
R<-resid(M1)

#plot

par(mfrow=c(2,1),mar=c(2,3,2,1))
qqnorm(R,pch=16,cex=.5)
qqline(R)
hist(R)

shapiro.test(R)
ks.test(R,"pnorm")
#####################################
# Test of Homogeneity of Variances

#plots
plot(R,fitted(M1),ylim=c(4.5,6))

?var.test

bartlett.test(R~PlantGrowth$group)
# nonparametric test
?fligner.test
fligner.test(R~PlantGrowth$group)

## uncorrelated residuals
#plots

ts.plot(R)
points(R,pch=16,cex=.5,col=4)
abline(h=0,lty=2,col=2)
####
library(lmtest)
?dwtest 

dwtest(M1)

### Kruskal-Wallis Rank Sum Test
 ?kruskal.test

kruskal.test(weight ~ group, data = PlantGrowth)


###################################
#Example2: ToothGrowth
###################################
ToothGrowth

ToothGrowth$dose

is.factor(ToothGrowth$dose)
table(ToothGrowth$dose)

#Convert the variable dose from a numeric to a factor variable
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
head(ToothGrowth)

is.factor(ToothGrowth$dose)

is.factor(ToothGrowth$supp)


library(ggplot2)
# Basic box plot
ggplot(ToothGrowth, aes(x = dose, y = len))+ geom_boxplot()

ggplot(ToothGrowth, aes(x = dose, y =len, col=dose))+ geom_boxplot()

ggplot(ToothGrowth, aes(x = dose, y = len,col=dose))+ geom_point()


??geom
# change defult theme
?theme_set()
?theme_bw

theme_set(theme_bw())
ggplot(ToothGrowth, aes(x = dose, y = len,col=dose))+ geom_boxplot()


theme_set(theme_classic())
ggplot(ToothGrowth, aes(x = dose, y = len,col=dose))+ geom_boxplot()



theme_set(theme_void())
ggplot(ToothGrowth, aes(x = dose, y = len,col=dose))+ geom_boxplot()


theme_set(theme_light())
ggplot(ToothGrowth, aes(x = supp, y = len,col=supp)) + 
	geom_boxplot() +geom_point(col="black",size=.8,alpha=.2)


ggplot(ToothGrowth, aes(x = dose, y = len, fill=supp)) + 
	geom_boxplot() 
	#+geom_dotplot(binaxis="y",stackdir='center',position=position_dodge(1))

######################  Generate Factor Levels #############
### gl
?gl

###one way anova: 5 treatment in 5 replications
gl(5,5,labels=c("15%","20%","25%","30","35"))


### two factors: A in 2 levels and A in 3 levels, all in 2 replications
A<-gl(2,1,12,labels=c("a1","a2"))
A

B<-gl(3,2,2*3*2,labels=c("b1","b2","b3"))
B

data.frame(A,B)
####### gl
species<-gl(3,1,4*3*2, labels=c("Species 1","Species 2","Species 3") )
#   Male  Male ... Male;  Female  Female ... Female
gender <- gl(2,3*4,2*3*4, labels=c("Male", "Female") )

#############
A<-rep(1:2,6)
A<-factor(A,labels=c("a1","a2"))
A

B<-rep(rep(1:3,each=2),2)
B<-factor(B,labels=c("b1","b2","b3"))
B

data.frame(A=A,B=B)


##############################################
## computing means by the two factors
aggregate(len ~ dose,data = ToothGrowth,FUN=mean)

aggregate(len ~ supp,data = ToothGrowth,FUN=mean)

aggregate(len ~ dose+supp,data = ToothGrowth,FUN=mean)

plot.design(len ~ dose+supp,data = ToothGrowth)

M1<-aov(len ~ dose*supp,data = ToothGrowth)

anova(M1)

model.tables(M1)

## Tukey Honest Significant Differences

TukeyHSD(M1)

library(agricolae)

out<-LSD.test(M1,"dose")
out

out<-duncan.test(M1,"dose")
out

out<-LSD.test(M1,c("dose","supp"))
out
plot(out)

?interaction.plot

interaction.plot(ToothGrowth$dose,ToothGrowth$supp,ToothGrowth$len)

interaction.plot(ToothGrowth$supp,ToothGrowth$dose,ToothGrowth$len)


### check residuals
#?? exercise


####################################
###### Central Limit Theorem #######
####################################

# x~binomial(5,p)
n=100
p=0.1


x<-rbinom(n,5,p)
xbar<-mean(x)

rep=1000
xbar<-0
for(i in 1:rep) {
	x<-rbinom(n,5,p)
	xbar[i]<-mean(x)
}

hist(xbar,col="gray80")
barplot(table(x))

layout(matrix(c(1,1,2,3),2))
hist(xbar,col="orange")
barplot(table(x))
qqnorm(xbar,pch=16,cex=.5)
qqline(xbar)

#####
?distributions
# x~cauchy(0,1)
n=100

x<-rcauchy(n)
xbar<-mean(x)

rep=1000
xbar<-0
for(i in 1:rep) {
	x<-rcauchy(n)
	xbar[i]<-mean(x)
}

hist(xbar,col="gray80")
hist(x)

layout(matrix(c(1,1,2,3),2))
hist(xbar,col="orange")
hist(x,col="gray80")
qqnorm(xbar,pch=16,cex=.5)
qqline(xbar)

## without for!! 

##by using matrix


#################
## by user function
 library(ggplot2)
clt.binom<-function(n=5,p=0.5){
	rep=1000
	xbar<-0
	
for(i in 1:rep) {
	x<-rbinom(n,5,p)
	xbar[i]<-mean(x)
}
layout(matrix(c(1,1,2,3),2))
hist(xbar,col="gray85")
barplot(table(x),col="green")
qqnorm(xbar,pch=16,col=4,cex=.5)
qqline(xbar,col="gray70")
}

clt.binom()

clt.binom(n=50)   
## paste
clt.binom(n=100,p=.99)
clt.binom(n=100,p=.01)

## density function with different samples

clt.binom.d<-function(n=c(5,30,50),p=0.5){
	rep=1000
	L<-length(n)
	xbar<-matrix(NA,ncol=L,nrow=rep)
for(i in 1:L) {
	x<-rbinom(n[i]*rep,5,p)
	xbar[,i]<-apply(matrix(x,n[i],rep),2,mean)
}

gg<-list()
hh<-list()
for(i in 1:L){
gg[[i]]<-ggplot(data=data.frame(x=xbar[,i]),aes(x=x)) + geom_density()+xlab(paste("n=",n[i]))
hh[[i]]<-ggplot(data=data.frame(x=xbar[,i]),aes(x=x)) + geom_histogram(bins=10,fill="gray80",col="red")+xlab(paste("n=",n[i]))
}

plot(density(xbar[,1]),ylim=c(0,3))
for(i in 2:L) lines(density(xbar[,i]),lty=i,col=i)
return(list(xbar=xbar,gg=gg,hh=hh))
}

aa<-clt.binom.d()

cowplot::plot_grid(plotlist = aa$gg)
cowplot::plot_grid(plotlist = aa$hh)


aa<-clt.binom.d(n=c(5,20,50,80,100,200))
cowplot::plot_grid(plotlist = aa$hh)

aa<-clt.binom.d(n=c(5,20,50,80,100,200),p=.05)
cowplot::plot_grid(plotlist = aa$gg)
cowplot::plot_grid(plotlist = aa$hh)

########################
W<-c(60,61,80,81,76,64,57,58,44,65)
hist(W)
length(W)
n<-3
x<-sample(W,n,replace=T)
xbar<-mean(x)
xbar

mean(W)


######################
gg1 <- ggplot()
gg1 <- gg1 + geom_density(aes(x=xbar[,1]),fill=1,alpha=.1)+xlab("xbar")

gg1
############
## chisq, t, cauchy
?dcauchy



###################################################
######### density plots for the some distributions
###################################################
## beta
set.seed(1)
n=1000
v1 <- rbeta(n, 0.5, 0.1)
v2 <- rbeta(n, 0.25, 0.3)
data1<-data.frame(x=c(v1, v2),g=c(rep("a", n), rep("b", n)))
head(data1)

gg <- ggplot(data1)
gg <- gg + geom_density(aes(x=x, y=..scaled.., fill=g), alpha=1/2)
gg

####################################
### normal and t
set.seed(1)
n=5000
v1 <- rnorm(n)
v2 <- rt(n, 1)
data1<-data.frame(x=c(v1, v2),g=c(rep("Normal", n), rep("t", n)))
head(data1)

gg <- ggplot(data1)
gg <- gg + geom_density(aes(x=x, fill=g), alpha=.1)+xlim(-10,10)
gg

####
set.seed(301)
n=100000
v1 <- rnorm(n)
gg <- ggplot()
gg <- gg + geom_density(aes(x=v1),fill="red",alpha=.1)+xlim(-10,10)
gg

v2 <- rt(n, 1)
gg <- gg + geom_density(aes(x=v2),lty=2,col=2)
gg

v3 <- rt(n, 5)
gg<-gg+ geom_density(aes(x=v3),lty=3,col=3)
gg

v4 <- rt(n, 30)
gg<-gg+ geom_density(aes(x=v4),lty=4,col=4)
gg

############################################################
   ##### different plots by ggplot for distributions #####
############################################################
p1 <- ggplot(data = data.frame(x = c(-4, 4)),aes(x)) + stat_function(fun = dnorm, n = 101, args = list(mean = 0, sd = 1)) + scale_y_continuous(breaks = NULL)+ stat_function(fun = dnorm, args = list(mean = 1, sd = 1),lty=2,col=3) 
p1    
###############

###############################
#### Discrete Analysis #####
###############################
?binom.test
#########
n<-15
x<-10
binom.test(x,n,alt="g")

#### 

# H0:p1=p2=p3 ?
n1<-50
n2<-60
n3<-45
x1<-28
x2<-40
x3<-20

prop.test(c(x1,x2,x3),c(n1,n2,n3))

#H0: p1=.6, p2=0.7 , p3=.7
p01=0.6
p02=0.7
p03=0.7
prop.test(c(x1,x2,x3),c(n1,n2,n3),p=c(p01,p02,p03))


?table
?as.table
?xtabs
#important: goodness-of-fit tests for continuous data : ks.test

#Pearson's Chi-squared Test for Count Data
#chi-squared contingency table tests and goodness-of-fit tests
# chisq.test
?chisq.test

### Example 5 page 229
## From Agresti(2007) p.39
tab <- as.table(rbind(c(279, 73, 225), c(165, 47, 191)))
tab
dimnames(tab) <- list(gender = c("F", "M"),
                    party = c("Democrat","Independent", "Republican"))
tab

M<-chisq.test(tab)
M
M$obs
M$expected
M$residuals

## glm ##

?glm
#### 5th competition 1383 page 62 ###

setwd("/Users/macbookpro/Downloads/competition Questions /5")

SSC1<-read.csv("5-SSC1.csv")

head(SSC1)
## data summary
summary(SSC1$age)
table(SSC1$sector)
table(SSC1$y)
dim(SSC1)

is.numeric(SSC1$age)
is.factor(SSC1$y)
is.factor(SSC1$sector)

SSC1$y<-as.factor(SSC1$y)
SSC1$sector<-as.factor(SSC1$sector)

is.factor(SSC1$y)
is.factor(SSC1$sector)

library(ggplot2)
ggplot(SSC1,aes(x=y,fill=sector))+geom_bar()
ggplot(SSC1,aes(x=y,fill=sector))+geom_bar(position = "dodge2")
# a)  model? y is discrete  
## glm ##

?glm

# The gaussian family accepts the links (as names) identity, log and inverse; 
# the binomial family the links logit, probit, cauchit, (corresponding to logistic, normal and Cauchy CDFs respectively) log and cloglog (complementary log-log); 
#the Gamma family the links inverse, identity and log; 
#the poisson family the links log, identity, and sqrt; 
#and the inverse.gaussian family the links 1/mu^2, inverse, identity and log.


FM<-glm(y~age+sector,data=SSC1,family=binomial(link="logit"))
FM
summary(FM)

# odds Ratio

coef(FM)

OR<-exp(coef(FM)[3])
OR

?anova.glm
anova(FM,test="Chisq")
anova(FM,test="LRT")

# check residuals

# reduce model
RM1<-glm(y~age,data=SSC1,family=binomial(link="logit"))
RM1
summary(RM1)


RM2<-glm(y~sector,data=SSC1,family=binomial(link="logit"))
RM2
summary(RM2)

#
anova(RM1,FM,test="LRT")

anova(RM2,FM,test="LRT")

##


M2<-glm(y~age+sector,data=SSC1,family=binomial(link="probit"))
M2
summary(M2)

M3<-glm(y~age+sector,data=SSC1,family=binomial(link="cloglog"))
M3
summary(M3)

# MSE 

MSE.FM<-mean(resid(FM)^2)

MSE.M2<-mean(resid(M2)^2)

MSE.M3<-mean(resid(M3)^2)
#AIC 
AIC(FM)
AIC(M2)
AIC(M3)
###
#b
M2<-glm(y~age+sector,data=SSC1,family=binomial(link="probit"))
M2
summary(M2)
confint(M2)
### c)

#predict

?predict.glm
New<-data.frame(age=50)

RM3<-glm(y~age,data=SSC1,family=binomial(link="probit"))
RM3

predict(RM3,New,type="link")

predict(RM3,New,type="response")

predict(RM3,New,type="terms")

######
#d)
M2<-glm(y~age+sector,data=SSC1,family=binomial(link="probit"))
M2

New<-data.frame(age=50,sector="1")
New
predict(M2,New,type="response")

New<-data.frame(age=50,sector="2")
New
predict(M2,New,type="response")


## e)


# full model and reduce model
anova(M1)
anova(M2)

?anova.glm

anova(RM1,FM,test="LRT")

anova(RM2,FM,test="LRT")
# 
#######################################
#### 7th competition 1385 page 95 ###
#######################################

setwd("/Users/macbookpro/Downloads/competition Questions /7")

data7<-read.csv("Prob2.csv")

data7
head(data7)

is.factor(data7$SEX)

data7$SEX<-as.factor(data7$SEX)
data7$EDUCATION <-as.factor(data7$EDUCATION)
data7$EXRCS<-as.factor(data7$EXRCS)
data7$SMOKING<-as.factor(data7$SMOKING)
data7$HTN<-as.factor(data7$HTN)

BMI<-data7$WEIGHT/(data7$HEIGHT/100)^2


max(BMI)
BMI1<-cut(BMI,breaks=c(0,25,29.99,40))
BMI1

data7$BMI<-factor(BMI1,labels=c("Normal","littleFat","Fat"))
head(data7)


#### a)

table(data7$BMI, data7$HTN)

chisq.test(table(data7$BMI, data7$HTN))

table(data7$EDUCATION, data7$HTN)

chisq.test(table(data7$EDUCATION, data7$HTN))

chisq.test(table(data7$SEX, data7$HTN))

# xtabs to create 
?xtabs
xtabs(~BMI+HTN,data=data7)

xtabs(~BMI+HTN+SEX,data=data7)


chisq.test(xtabs(~BMI+HTN,data=data7))

## Cochran-Mantel-Haenszel test  # 3-dimensional array
mantelhaen.test(xtabs(~BMI+HTN+SEX,data=data7))

mantelhaen.test(xtabs(~EDUCATION+HTN+BMI,data=data7))

mantelhaen.test(xtabs(~EDUCATION+BMI+HTN,data=data7))


###### b)

#1)
#glm family =binomial



data7$SEX<-factor(data7$SEX,labels=c("Male","Fmale"))
data7$SEX

data7$EXRCS<-factor(data7$EXRCS,labels=c("NO","Yes"))
data7$EXRCS

data7$SMOKING<-factor(data7$SMOKING,,labels=c("NO","Yes"))

M<-glm(HTN~SEX+SMOKING+BMI+EXRCS,data=data7,family=binomial(link="logit"))
M

summary(M)

#define levels and reference levels 

?relevel


data7$SEX<-relevel(data7$SEX,ref="Fmale")
data7$SEX

data7$EXRCS<-relevel(data7$EXRCS,ref="Yes")
data7$EXRCS

M<-glm(HTN~SEX+SMOKING+BMI+EXRCS,data=data7,family=binomial(link="logit"))

summary(M)


#2) 
# Odds Ratio 

OR<-exp(coef(M))
OR

1/0.08684245

## confidence interval for  OR
exp(confint(M))



### 3

M<-glm(HTN~SEX+SMOKING+BMI+EXRCS+AGE,data=data7, family=binomial(link="logit"))

summary(M)

# Odds Ratio 

OR<-exp(coef(M))
OR

## confidence interval for  OR
exp(confint(M))

#### 4)


xnew<-data.frame(SEX="Male", SMOKING="Yes", BMI="Fat", EXRCS="Yes", AGE=59)
xnew

predict(M,xnew,type="response")




##### glm

M<-glm(HTN~BMI+EDUCATION+SEX,data=data7, family=binomial(link="logit"))

M
summary(M)

anova(M,test="Chisq")


#####################################
##### Time Series Analysis by R #####
#####################################

######## Data definition
# example 1:
xt<-c(8,9,12,15,25,30,32,33,20,25,18,10,7,10,11,14,26,28,33,32,22,26,17,9)
xt
length(xt)
yt<-ts(xt)
yt
plot(yt)
yt<-ts(xt,start=c(2010,1),end=c(2012,12),freq=12)
yt
plot(yt)

# adding Errors to data!

#Example 2:
t<-1:36
zt<-yt*t
zt<-ts(zt,start=c(2010,1),end=c(2012,12),freq=12)
plot(zt)

par(mfrow=c(2,1),mar=c(2,2,2,2))
plot(yt)
plot(zt)
############
# Classical Seasonal Decomposition by Moving Averages

?decompose

examp1.add<-decompose(yt,type="additive")
plot(examp1.add)

examp1.mult<-decompose(yt,type="multiplicative")
plot(examp1.mult)

# camparing two decompose: additive and  multiplicative


## additive
names(examp1.add)

examp1.add$seasonal
examp1.add$trend
examp1.add$random

examp1.add$figure  # seasonal estimates


fitted1<-examp1.add$seasonal+ examp1.add$trend

fitted1

plot(yt)
lines(fitted1,lty=2,col=2)

# MSE additive model Yt=St+Tt+Et
resid<-examp1.add$random
MSE1.add<-mean(resid^2,na.rm=T)
MSE1.add

### multiplicative
fitted2<-examp1.mult$seasonal*examp1.mult$trend

fitted2

lines(fitted2,lty=3,col=4)
# MSE multiplicative model Yt=St*Tt*Et
resid<-examp1.mult$random
MSE1.mult<-mean(resid^2,na.rm=T)
MSE1.mult

######

predict(examp1.add)
forecast
##### 
examp2.add<-decompose(zt,type="additive")
plot(examp2.add)

examp2.mult<-decompose(zt,type="multiplicative")
plot(examp2.mult)

resid<-examp2.add$random
MSE2.add<-mean(resid^2,na.rm=T)
MSE2.add

resid<-examp2.mult$random
MSE2.mult<-mean(resid^2,na.rm=T)
MSE2.mult

## Seasonal Decomposition of Time Series by Loess

?stl

exam1<-stl(yt, s.window="period")
exam1
plot(exam1)
library(forecast)
plot(predict(exam1),36)
plot(predict(exam1))

exam2<-stl(zt, s.window="period")
exam2
plot(exam2)
plot(predict(exam2))
plot(predict(exam2),36)

# MSE

names(exam1)
resid<-exam1$time.series[,"remainder"]
MSE1.stl<-mean(resid^2,na.rm=T)
MSE1.stl

resid<-exam2$time.series[,"remainder"]
MSE2.stl<-mean(resid^2,na.rm=T)
MSE2.stl

####

library(forecast)
forecast(exam1,h=2*12)

plot(forecast(exam1,h=3*12))

plot(forecast(exam2,h=3*12))


monthplot(yt)
monthplot(zt)


forecast::seasonplot(yt) # package forecast
forecast::seasonplot(zt) # package forecast


############################
#### AirPassengers data in R
############################
data(AirPassengers)
AP<-AirPassengers
AP
plot(AP)
is.ts(AP)
?stl
fit1 <- stl(yt, s.window="period")
fit1
plot(fit1)
fit2 <- stl(AP, s.window="period")
fit2
plot(fit2)
library(forecast)
?ets

########################
#HoltWinters
?HoltWinters

#exam1:yt
m1<-HoltWinters(yt,seasonal = "additive")
m1
plot(m1)
m2<-HoltWinters(yt,seasonal = "multiplicative")
m2
plot(m2)
par(mfrow=c(1,2))
plot(m1,main="m1")
plot(m2,main="m2")
names(m1)
m1$SSE
m2$SSE
m1$fitted
sum((yt-(m1$fitted)[,1])^2)
sum((yt-(m2$fitted)[,1])^2)

##exam2:zt
m1<-HoltWinters(zt,seasonal = "additive")
m1
plot(m1)
m2<-HoltWinters(zt,seasonal = "multiplicative")
m2
plot(m2)
par(mfrow=c(1,2))
plot(m1,main="m1")
plot(m2,main="m2")
names(m1)
m1$SSE
m2$SSE


################
m3<-HoltWinters(AP,seasonal = "additive")
m3
m4<-HoltWinters(AP,seasonal = "multiplicative")
m4
par(mfrow=c(1,2))
plot(m3,main="m3")
plot(m4,main="m4")
m3$SSE
m4$SSE

names(m4)
yhat<-m4$fitted[,"xhat"]
resid<-AP-yhat
MSE.AP.m4<-mean(resid^2)
MSE.AP.m4
sum(resid^2) # SSE

yhat<-m3$fitted[,"xhat"]
resid<-AP-yhat
MSE.AP.m3<-mean(resid^2)
MSE.AP.m3

sum(resid^2)  # SSE

################
library(forecast)
?forecast
f1<-forecast(m1,h=2*12)
plot(f1)
f2<-forecast(m4,h=5*12)
plot(f2)
####################
##ARIMA##
layout(matrix(c(1,2,1,3),2))
plot(yt)
plot(m1)
plot(m2)

layout(matrix(c(1,2,3,2),2))
plot(m1)
plot(yt)
plot(m2)
########################################
############ ARIMA Models ###############
?arima
par(mfrow=c(1,2))
acf(yt,lag.max=10,type='cor',main='ACF')
acf(yt,lag.max=10,type='par',main='PACF')
## fitting arima model in R with library(forecast)
#model1=ar(1)
fit1=arima(yt, order=c(1,0,0))
fit1
fit2=arima(yt, order=c(2,0,0))
fit2
fit3=arima(yt, order=c(3,0,0))
fit3
summary(fit3)
fit4=arima(yt, order=c(3,0,1))
summary(fit4)
fit5=arima(yt, order=c(3,0,2))
summary(fit5)
fit6=arima(yt, order=c(3,0,3))
summary(fit6)

fit11<-arima(yt, order=c(3,0,2),seasonal = list(order = c(1, 0, 0)))
fit11

fit7<-auto.arima(yt)
summary(fit7)

yt1<-diff(yt,lag=12)
plot(yt1)
plot(yt)

fit8=arima(yt1, order=c(0,0,1))
summary(fit8)

fit9=arima(yt, order=c(0,0,1),seasonal=c(0,1,0))
summary(fit9)

fit10=arima(yt, order=c(0,0,1),seasonal=c(0,1,3))
summary(fit10)
f1<-forecast(fit7,h=2*12)
f2<-forecast(fit10,h=2*12)
layout(matrix(c(1,2),2))
plot(f1)
plot(f2)
############
layout(matrix(c(1,2,1,3),2))
plot(AP)
acf(AP,type='cor',main='ACF')
acf(AP,type='par',main='PACF')
AP1<-diff(AP,1)
plot(AP1)
AP2<-diff(AP1,12)
plot(AP2)
layout(matrix(c(1,2,1,3),2))
plot(AP2)
acf(AP2,type='cor',main='ACF')
acf(AP2,type='par',main='PACF')
auto.arima(AP)
auto.arima(AP2)
## fitting arima model in R with library(forecast)
#model1=ar(1)
fit1=arima(AP, order=c(2,1,1),seasonal=c(0,1,0))
summary(fit1)

fit2=arima(AP, order=c(2,0,0))
fit2
fit3=arima(AP, order=c(3,0,0))
fit3
summary(fit3)

f1<-forecast(fit1,h=4*12)
plot(f1)


#########################################################################################################
###          Regression         ###
##################################
#####Mesal 1 - page 143
x<-c(42,35,50,43,48,62,31,36,44,39)
y<-c(12,8,14,9,11,16,7,9,12,10)
plot(x,y,pch=16)
cor(x,y)
f<-cor.test(x,y)
f$estimate
cor.test(x,y,method = "spearman")
rx<-rank(x)
ry<-rank(y)
cor(rx,ry)
cor.test(x,y,method = "kendall")

library(car)
scatterplot(y~x)


#####Mesal 2 - page 148
m<-c(31,37,28,29,38,42,46,33,45,44,47,42)
f<-c(385,400,395,365,475,440,490,420,560,525,480,510)
lm(m~f)
k<-lm(m~f)
k$coefficients
k$df.residual
class(k)
summary(k)
anova(k)
k1<-lm(m~(-1)+f)
k1$coefficients
library(car)

#####Mesal 3 - page 150
y<-c(49,28,19,21,26,17,30,33,13,9,13,19,27,6,13,11,21,21,31,21,10,21,19,26,22,34)
x1<-c(18,14,11,9,10,8,11,13,7,2,6,9,12,4,7,7,7,5,11,10,5,8,9,11,11,13)
x2<-c(21,12,8,9,11,7,13,14,5,4,6,8,11,2,5,4,9,9,13,9,4,8,7,11,9,14)
n<-data.frame(y,x1,x2)
w<-lm(y~x1+x2,data = n)
summary(w)
coef(w)
anova(w)
#####Mesal 4 - page 154
r<-w$residuals
yh<-w$fitted.values
plot(n$x1,r)
plot(n$x2,r)
plot(yh,r)
library(lmtest)
dwtest(y~x1+x2,data = n)
bptest(y~x1+x2,data = n)
vif(w)

# How to do hypothesis in Beta?
##### H0: beta1=beta2

FM<-lm(y~x1+x2,data = n)

# under H0
RM<-lm(y~I(x1+x2),data=n)

anova(RM,FM)

### H0: beta1=beta2/2 ?
## H0: beta2=2*beta1 
# under H0
RM<-lm(y~I(x1+2*x2),data=n)
anova(RM,FM)

## H0: beta2=10*beta1 
# under H0
RM<-lm(y~I(x1+10*x2),data=n)
anova(RM,FM)

### H0: beta1=.5 
RM<-lm(y~offset(0.5*x1)+x2,data=n)
anova(RM,FM)

### H0: beta1=.5 
RM<-lm(y~offset(0.5*x1)+x2,data=n)
anova(RM,FM)


library(car)

?linearHypothesis

linearHypothesis(FM, c("(Intercept) = 0", "x1= 1"))

linearHypothesis(FM, c( "x1= 0.5"))
linearHypothesis(FM, c( "x1=x2"))

linearHypothesis(FM, c("(Intercept) = 0", "x2= 10*x1"))

#####Mesal 5 - page 155
k<-lm(m~f)
plot(k)
w<-lm(y~x1+x2,data = n)
plot(w)
r2<-residuals(k)
qqplot(r2)
qqline(r2)
r3<-residuals(w)
qqplot(r3)
qqline(r3)
#####Mesal 6 - page 159
y<-c(30,34,33,26,41,4,5,20,31,38,43,47,45,45,11,10,30,29,23,16,37,50,36,54,44)
x1<-c(84,84,79,81,84,74,73,75,84,86,88,90,88,88,81,79,84,84,84,77,87,89,89,93,93)
x2<-c(85,86,83,83,88,77,78,84,89,91,91,94,94,92,87,83,87,87,88,83,92,92,94,92,93)
x3<-c(398,345,388,406,379,478,462,464,430,406,393,385,405,392,448,436,392,392,398,431,379,393,394,386,389)
A<-data.frame(y,x1,x2,x3)
y<-A[,1]
x<-A[,-1]
#install.packages("leaps")
library(leaps)
leaps(x,y,method="adjr2")
leaps(x,y,method="Cp")$Cp
leaps(x,y,method="r2")$r2

#####Mesal 7 - page 163



##########################
### Example real data regression
##########################

# DATA
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
#USairpollution

head(USairpollution)

SO2=USairpollution[,1]

temp=USairpollution[,2]
manu=USairpollution[,3]
popul=USairpollution[,4]
wind =USairpollution[,5]
precip =USairpollution[,6]
predays=USairpollution[,7]




library(car)
scatterplotMatrix(USairpollution)

cov(USairpollution)
cor(USairpollution)

model=lm(SO2~.,data=data.frame(USairpollution))
coef(model)
confint(model)
predict(model)
predict(model,int="c")
yhat=fitted(model)
yhat
ei=resid(model)
ei

summary(model)
anova(model)
vif(model)
################ Classic assumption with ei
######## 1) normality
#### a
qqnorm(ei);qqline(ei)

#### b
shapiro.test(ei)

######## 2) Homoscedasticity
#### a
par(mfrow=c(2,3))
plot(temp,ei);abline(h=0)

plot(manu,ei);abline(h=0)

plot(popul,ei);abline(h=0)

plot(wind,ei);abline(h=0)

plot(precip,ei);abline(h=0)

plot(predays,ei);abline(h=0)

#### b
library(lmtest)
bptest(model)

######## 3) Autocorrelation
#### a
ts.plot(ei);abline(h=0, lty=2,col=2)

#### b
library(lmtest)
dwtest(model,data=data.frame(USairpollution))


################ Classic assumption with ei
######## 2) Homoscedasticity
#### a
par(mfrow=c(2,3))
plot(temp,ei);abline(h=0)

plot(manu,ei);abline(h=0)

plot(popul,ei);abline(h=0)

plot(wind,ei);abline(h=0)

plot(precip,ei);abline(h=0)

plot(predays,ei);abline(h=0)



################ Classic assumption with ri
######## 2) Homoscedasticity
#### a
ri=rstandard(model)
par(mfrow=c(2,3))
plot(temp,ri);abline(h=0)

plot(manu,ri);abline(h=0)

plot(popul,ri);abline(h=0)

plot(wind,ri);abline(h=0)

plot(precip,ri);abline(h=0)

plot(predays,ri);abline(h=0)


################ Classic assumption with si
######## 2) Homoscedasticity
#### a
si=sqrt(abs(ri))
par(mfrow=c(2,3))
plot(temp,si);abline(h=0)

plot(manu,si);abline(h=0)

plot(popul,si);abline(h=0)

plot(wind,si);abline(h=0)

plot(precip,si);abline(h=0)

plot(predays,si);abline(h=0)

##
plot(ei,yhat)

###### leverage
### a
lm.influence(model)
lev = hat(model.matrix(model))
plot(lev,ylim=c(-0.5,1),pch=16,col=4,cex=.5)
text(lev,abbreviate(row.names(USairpollution)),cex=0.6,col="gray70")
n=41;p=7 ;astane=2*p/n
abline(h=astane,col="red",lty=4,lwd=.8)

a=lev>(2*p/n)
a
b=lev[lev>(2*p/n)]
b

### b ==> Outlier with rstsndard
ri<-rstandard(model)
plot(1:41,ri,ylim=c(-4,4),pch=16,col=4);abline(h=c(-2,2),col="red",lty=2)
identify(1:41,ri, labels=names(ri),cex=.6,col="gray50")

##### ==> Outlier with rstudent
plot(1:41,rstudent(model),ylim=c(-4.5,4.7),pch=16);abline(h=c(-2,2),col="red",lty=2)
identify(1:41,rstudent(model),labels=names(ri),cex=.6,col="gray50")

library(car)
outlierTest(model)


###### 
par(mfrow=c(2,2))
plot(model)
p=7
abline(v=(2*(p+1))/n,lty=2,col="blue",lwd=2)
abline(h=c(-2,2),col="green",lwd=2)

#### statiticts for outlier
dffits(model)
d=dffits(model)[dffits(model)>(2*sqrt(p/n))]
d
cooks.distance(model)
e=cooks.distance(model)[cooks.distance(model)>(4/n-2)]
e
 

###################################################
## removing outlier
###################################################
USairpollution.new=USairpollution[-c(32,33),]
dim(USairpollution.new)
library(car)
scatterplotMatrix(USairpollution.new)

cov(USairpollution.new)
cor(USairpollution.new)

model.new=lm(SO2~.,data=data.frame(USairpollution.new))
coef(model.new)
confint(model.new)
predict(model.new)
predict(model.new,int="c")
yhat.new=fitted(model.new)
yhat.new
ei.new=resid(model.new)
ei.new

summary(model.new)
anova(model.new)


#######################################################
##### Choose a model by AIC in a Stepwise Algorithm
######################################################

####### Goodness of fit ==> AIC & BIC
n=41
backAIC=step(model,direction="backward")
backBIC=step(model,direction="backward", k=log(n))       


null<-lm(SO2~1) # null model
FM<-lm(SO2~temp+manu+popul+wind+precip+predays) # full model

# backward method
step(FM,scope=list(lower=null),direction="backward")

# forward methods
step(null,scope=list(upper=FM),direction="forward")


# stepwise method
step(null,scope=list(lower=null,upper=FM),direction="both" )

##  Mallows' CP
#scale=(summary(full)$sigma)^2 

# stepwise method by Cp Mallows
MSE<-(summary(FM)$sigma)^2
MSE
sum(resid(FM)^2)/(n-p) # computing  MSE by residuals

step(null,scope=list(lower=null,upper=FM),scale=MSE,direction="both" )


BIC=step(model,direction="both",k=log(n) )



####### Goodness of fit ==> Cp & R2adj

library(leaps)
data1=cbind(SO2,temp,manu,popul,wind,precip,predays)
A=data1[,-1]
B=data1[,1]
Cpresult=leaps(A,B,method="Cp")
Cpresult

length(Cpresult$size)

R2adjresult=leaps(A,B,method="adjr2")
R2adjresult

length(R2adjresult$size)
max(R2adjresult$adjr2)


library(leaps)
b=regsubsets(A,B)
rs=summary(b)
par(mfrow=c(1,2))
plot(1:6,rs$adjr2,xlab="Subset Size",ylab="Adjusted R-squared")
library(car)
subsets(b,statistic="adjr2")


##################################################
############################## Roykard tabdil
############### halat 1 ==> just response
########## method inverse response plot
library(alr3)
inverseResponsePlot(model)
SO2new<-SO2^1.264841
model1<-lm(SO2new~temp+manu+popul+wind+precip+predays)
scatterplotMatrix(cbind(SO2new,temp,manu,popul,wind,precip,predays))

############### halat 2 ==> pridictors
library(car)
library(MASS)
boxCox(model,plotit=T)
MF<-powerTransform(cbind(temp,manu,popul,wind,precip,predays))
summary(MF)

tempnew=temp^-1.2352
manunew=manu^0.2036
populnew=popul^0.2030
windnew=wind^0.5571 
precipnew=precip^1.6024
predaysnew=predays^0.2548

model2<-lm(SO2~tempnew+manunew+populnew+windnew+precipnew+predaysnew)
scatterplotMatrix(cbind(SO2,tempnew,manunew,populnew,windnew,precipnew,predaysnew))


############ halat 3 ==> all varible
######## roykard 1
library(car)
library(MASS)
boxCox(model,plotit = TRUE)
B<-powerTransform(cbind(temp,manu,popul,wind,precip,predays))
summary(B)

tempnew=temp^-1.2352
manunew=manu^0.2036
populnew=popul^0.2030
windnew=wind^0.5571 
precipnew=precip^1.6024
predaysnew=predays^0.2548

library(alr3)
inverseResponsePlot(model)

ynew<-Y^1.264841

model3<-lm(SO2new~tempnew+manunew+populnew+windnew+precipnew+predaysnew)
scatterplotMatrix(cbind(SO2new,tempnew,manunew,populnew,windnew,precipnew,predaysnew))


######## roykard 2
library(MASS)
boxcox(model,plotit = TRUE)
D<-powerTransform(cbind(SO2,temp,manu,popul,wind,precip,predays))
summary(D)

SO2new=SO2^-0.4241
tempnew=temp^-0.7865
manunew=manu^0.2764
populnew=popul^0.2321
windnew=wind^1.0969
precipnew=precip^1.6428
predaysnew=predays^0.1791

model3<-lm(SO2new~tempnew+manunew+populnew+windnew+precipnew+predaysnew)
scatterplotMatrix(cbind(SO2new,tempnew,manunew,populnew,windnew,precipnew,predaysnew))