# R studio
## 回歸上課筆記

### 載入資料
getwd()
setwd("C:\\Users\\user\\Downloads")
#複製下載的位置
getwd()
mf1 <- read.table("myfile1.txt",header=T)
#下載>組合管理>資料與搜尋>檢視>把隱藏已知檔案之附檔名勾勾取消
mf1a <- read.table("file:///C:/Users/user/Downloads/myfile1.txt"
,header=T)
#把file1的檔案右鍵按複製
mf2 <- read.table("myfile2.txt"
,skip=3 #跳過前x筆資料
,sep="," #用逗點分隔
,header=T) #header=參數
mf3 <- read.table("wages_pp.txt"
,header = T
,sep=",")
pizza <- read.csv("pizza.csv",header=T)
alv <- read.csv("https://data.montgomerycountymd.gov/api/views/4tja-rkhg/rows.csv?accessType=DOWNLOAD")
write.csv(alv,"Myalv.csv",row.names = F)
#把資料變成自己的
write.table(pizza,"Mypizza.txt",row.names = F)

#### 小考1
setwd('C:\\Users\\user\\Downloads')
getwd()
Canc. <- read.table('Cancer.csv',header=T,sep=",")

plot(PSA~Cv,data=Canc.)

set.seed(234)
n=nrow(Canc.)
Sindex=sample(n,round(n*0.9))
Train=Canc.[Sindex,]
Test=Canc.[-Sindex,]
mod1=lm(PSA~Cv,data=Train)
summary(mod1)

confint(mod1,level=0.95)
range(Train$Cv)
#有95%的信心，當cv每增加一單位，PSA會增加2.37%~4.16%，B0因x沒包含0，所以無意義

### 畫圖/信賴區間
setwd('C:\\Users\\user\\Downloads\\新資料夾\\HardnessExample')
Hardness=read.table('hardness.txt',header=T)

boxplot(Hardness$Temp,col='blue')
boxplot(Hardness$Hard,col='blue')
stem(Temp)
stem(Hard)
plot(Hardness$Temp,Hardness$Hard,col='2')
cor(Hardness$Temp,Hardness$Hard)
plot(Hard~Temp,data=Hardness,col='2')

x <- Hardness$Temp
y <- Hardness$Hard
xbar <- mean(x)
ybar <- mean(y)

sxx <- sum((x-xbar)^2)
sxy <- sum((x-xbar)*(y-ybar))
b1 <- sxy/sxx
b0 <- ybar-b1*xbar

mod1 = lm(Hard~Temp,data=Hardness)
summary(mod1)
#residual=殘差=yi-yi(hat)

range(Hardness$Temp)
#溫度不包含0，所以B0沒有解釋意義
attributes(mod1)

ei <- Hardness$Hard-mod1$fitted.values
s2 <- sum(ei^2)/(14-2)
s <- sqrt(s2)

anova(mod1)
confint(mod1,level=0.95)
confint(mod1,level=0.99)
dnew <- data.frame(Temp = c(35,55,65))
#溫度對於彈簧硬度的平均值的信賴區間
predict(mod1,dnew,interval="confidence",level=0.95) #prediction and corresponding CI for mean response#預測平均值的信賴區間
predict(mod1,dnew,interval="prediction",level=0.95)
#prediction and corresponding CI for new obs
#單一預測值的信賴區間
#(36)fit=彈簧平均硬度 有95%的信心彈簧硬度會落在Temp35=(48.68,50.96);Temp55=(23.35,25.64)
#(37)當Temp35=我有95%信心「預測」彈簧硬度會落在Temp35=(48.68,50.96);Temp55=(23.35,25.64)
library(visreg)
visreg(mod1)

### 線性回歸
library(car)
#library(lmtest)

#Reading the data as a data frame in R
data(Prestige)
View(Prestige)

contrasts(Prestige$type) #編碼
m1 <- lm(prestige~.,data=Prestige) #.=剩下的其他變數
summary(m1)
round(coef(m1),2)
#迴歸係數只取到小數第二位。在其他自變數不變的情況下，職業類別是prof他的Y(社會經濟地位)會比職業類別是bc的高10.77個單位
#在其他自變數不變的情況下，職業類別是wc他的Y(社會經濟地位)會比職業類別是bc的高0.29個單位
round(confint(m1,level=0.95),1)
#係數的信賴區間，四捨五入到小數第一位。在其他自變數不變的情況下，我有百分之95信心，教育程度每增加一分，社會經濟地位會平均增加2.6到5.2分
#在其他自變數不變的情況下，我有百分之95信心，職業類別是prof社會經濟地位的平均值會比職業類別藍領為高1.5到20.1分
summary(m1)#教育程度不顯著的影響一個人的社會經濟地位(因為p-value0.66大於0.5=不拒絕)
#m1=full model
m0 <- lm(prestige~.-type,data=Prestige)#m0=reduced model
summary(m0)

anova(m0,m1)
#p-value0.0018<0.5所以拒絕B21B22同時為0的虛無假設，所以type整體是有顯著的影響
Prestige$type <- relevel(Prestige$type,ref="prof")
contrasts(Prestige$type)

m1a <- lm(prestige~.,data=Prestige)
summary(m1a)
#在其他自變數不變的情況下，職業類別是藍領社會經濟地位的平均值會比職業類別prof為低10.77分
#在其他自變數不變的情況下，職業類別是白領社會經濟地位的平均值會比職業類別prof為低10.48分
#Multiple R-squared: 0.841=模型所解釋的變異占總變異0.841
#< 2.2e-16拒絕b0以外至少由一個B顯著不為0

PrestigeC <- na.omit(Prestige)#有遺失值
m1a <- lm(prestige~.,data=PrestigeC)
summary(m1a)

m1b <- lm(prestige~1,data=PrestigeC)
summary(m1b)

anova(m1b,m1a)

MSR <- 23841/6
MSE <- 4505.9/91
F1 <- MSR/MSE
F1

### 各式圖及模型的解釋

library(car)
library(lmtest)
library(nortest)
setwd("C:\\Users\\user\\Downloads")
IOWAtest=read.table("iowatest.txt",header=T,sep=",")
IOWAtest=IOWAtest[,-1]
stem(IOWAtest$Poverty)
stem(IOWAtest$Test)
barplot(table(IOWAtest$City))
contrasts(IOWAtest$City)
#n=dim(NGasData)[1]
#Sindex=sample(n,round(n*0.7))
#Train=NGasData[Sindex,]
#Test=NGasData[-Sindex,]

set.seed(123)
ind1 <- sample(nrow(IOWAtest),110)
TrainD=IOWAtest[ind1,]
TestD=IOWAtest[-ind1,]

pairs(Test~.,data=TrainD)
#散步圖矩陣
cor(cbind(TrainD$Poverty,TrainD$Test))
#有很強的負相關:-0.834
summary(TrainD$City)
mod1=lm(Test~Poverty+City,data=TrainD)
#.=test以外的所有變數都要當自變數#Full Model
mod1=lm(Test~.,data=TrainD)
summary(mod1)
#在其他自變數不變的情況之下，來自CityDavenport的學生考試平均值會比參考類別Cedar Rapids的學生低8.1951分
contrasts(TrainD$City)
#參考類別:Cedar Rapids(因0 0 0 0 0)
TrainD$City=relevel(TrainD$City,ref="Iowa City")
#在其他自變數不變的情況之下，來自CityDavenport的學生考試平均值會比參考類別Iowa City的學生低8.1951分
#有**(小)拒絕H0跟City有關的dummy varable參數皆為0，City全部留下

mod2<-lm(Test~.-City,data=TrainD)
#Reduced Model
summary(mod2)
anova(mod1,mod2)
#以上複習
e=residuals(mod1)
es=rstandard(mod1)
yhat=fitted.values(mod1)

plot(yhat,es,col='2',pch=16)
#plot(yhat,es,col='2')
abline(h=0) #e-yhat plot
residualPlot(mod1,type="rstandard",quadratic=F)
residualPlots(mod1,type="rstandard",quadratic=F)#從最後一個圖開始看(跟剛剛的圖一樣)，假設Fitted Value有問題往上看其他圖看是哪個變數出問題
#如果中心點沒有躺在線上就代表出了問題

resettest(mod1,power=seq(0.5,2,by=0.5),type = "regressor") # Ramsey's reset goodness of fit test
#除了模型我再加一個模型，後從0.5到2每次增加一單位次方(這裡是0.5)，如果成功代表元模型不好(這裡沒有比較好因為0.18>0.05)
#不一定要會
raintest(mod1)#Rainbow goodness of fit test
bptest(mod1) #Breusch and Pagan test 檢定變異數是否一致

qqPlot(mod1)
#常態性(呈現一直線)，不要跳出虛線太多點就接受他是常態
lillie.test(es)
#KS test for normality 假設服從常態結果拒絕
shapiro.test(es)
#Shapiro-Wilk Normality Test 假設服從常態然後不拒絕

plot(es,type = "l",col='2')
dwtest(mod1)
#Durbin-Watson test(alternative hypothesis: true autocorrelation is greater than 0=右尾)
#接受前後不相關的假設，接受資料為獨立
acf(es,ci=0.99)
#虛線是convidence interval
#detach(Train)

confint(mod1,level=0.95)
#我有95%信心，poverty每增加一個單位考試成績就下降0.56-0.42個單位
pv1=predict(mod1,TestD,interval="confidence")
plot(TestD$Test,pv1[,1])
pv2=predict(mod1,TestD,interval="prediction")

### Anova表

library(car)
library(lmtest)
setwd("C:\\Users\\user\\Downloads")
GasData=read.table("gasconsumption.txt",header=T)
m0 <- lm(MPG~GPM+WT+DIS,data=GasData)
summary(m0)
#n=dim(GasData)[1]
#Sindex=sample(n,round(n*0.7))
#Train=GasData[Sindex,]
#Test=GasData[-Sindex,]

attach(GasData)
stem(MPG)
stem(WT)
stem(DIS)
stem(NC)

pairs(MPG~GPM+WT+DIS,data=GasData) #散布圖矩陣
scatterplotMatrix(GasData)
scatterplot(MPG~WT,, id.method="identify")

M1=lm(MPG~WT+DIS+NC)
#複回歸模型，x=WT,DIS,NC
summary(M1)#MPGhat=59.61-14.43WT+0.07DIS-1.09NC。14.43=在其他自變數不變的情況之下，車重每增加一個單位，MPG平均會降低14.43個單位
anova(M1)#除非所有的自變數包含0，截距項才有意義。在所有的自變數為0的情況之下，MPG的平均估計值為59.61；否則只要任一自變數不包含0就沒意義
#DIS會顯著的影響MPG(anova裡的pr)，NC不顯著影響
#(T-TEST:單一變數是否為0的虛無假設SUMMARY)
#(所有變數是否全部為0的虛無假設ANOVA)

y <- GasData$MPG
ybar <- mean(GasData$MPG)
SST <- sum((y-ybar)^2)
SSE <- sum(M1$residuals^2)
SSR <- SST-SSE

M0 <- lm(MPG~1,data=GasData)
summary(M0)
anova(M0,M1)#SST=1586.09(df=37),SSE=193.04,sum of square=SSR

F1 = (1393/3)/(193.04/34)#F

confint(M1,level=0.95)
#有百分之95的信心，在其他的自變數不變的情況之下，車重每增加一個單位，油耗的平均值會降低18到10個單位之間
confint(M1,level=0.99)

detach(GasData)
WT=c(2.5,3.5,4.5)
DIS=c(250,350,450)
NC=c(4,5,6)
New=data.frame(cbind(WT,DIS,NC))
predict(M1,New,interval="confidence") #prediction and corresponding CI for mean response 平均值的信賴區間
predict(M1,New,interval="prediction") #prediction and corresponding CI for new obs 單一值的信賴區間

### 矩陣/資料處理

x <- 1:8
dim(x) <- c(2,4)
x #經過dim之後變成矩陣
x <- matrix(1:8,2,4,byrow=F) #把1 2 列合併的結果 #先把行填滿 再填列
x <- matrix(1:8,2,4,byrow=T)# T:把列先填滿
v1<-c(1,2,3,4)
v2 <- c(5,6,7,8)
m1 <- rbind(V1,V2) #用列填滿
m2 <-cbind(v1,v2) #用行填滿
m1<-matrix(rnorm(1000),10,100)
##1000筆資料
m2<-matric(rnorm(2000,2,1),20,100)
##20列 100個變數
m <-rbind(m1,m2)
v1[3]
m[30,1] #取30列第1格個元素
m[,1] #要一整列 就空格
m[30,] #取第30列 一整行
mm[seq(1,30,2),seq(2,30,2)] #seq1~30 每隔2取數
mm1<- m[c(1,7,8),c(9,4,8,7)]#1 7 8 列的 9 4 8 7 行
mm2 <-m[m[,1]>=2,1] #m大於等於2的值
mm3 <-m[-30,] #不要第30列 前面加負號
mm4 <-m[,-seq(2,100,2)] #2~30行 每隔2去掉
gender <- rep(c("male","female"),c(200,300))
income <- runif(500,min=22000,max=55000) #500筆薪水資料 最低22000 最高55000
mydate <- cbind(gender,income) #不同的資料型態全部變成字串
mydate1 <- cbind(index=1:500,gender,income)
mean(mydata[,2])
mydata2 <- data.frame(gender,income)
summary(mydata2)
mean(mydata2[,2])
mean(mydata2$income) #和36 結果一樣
mydata2$Index <- 1:500 #對資料做編號的指標
data()
data(USArrests)
View(USArrests)
is.data.frame(USArrests)
USArrests["Florida",] #取佛羅里達的資料
row.names(USArrests) #取出來是向量
plot(Murder~Assault,data=USArrests)
text(USArrests$Assault,USArrests$Murder,labels = row.names(USArrests),pos=1)
colnames(USArrests)
USArrests[,c("Murder","Assault")]
SN<-row.names(USArrests) #取MURDER大於10的洲名
SN[USArrests$Murder>=10]

#從USArrests 取murder>=10
SN[USArrests$Assault>=mean(USArrests$Assault)]
USArrests[,2]#取第2列所有資料
USArrests[3,]#取第3行所有資料
newd <- USArrests[,-3]
newd <- USArrests[,!colnames(USArrests)%in%c("Rape")] #削掉rape #colnames 是變數名稱 !=not
test<-matrix(rnorm(21),7,3) #21筆資料7列3行
testd<-data.frame(test)
colnames(testd)<-c("name1","name2","names3")
testd
m1<- matrix(1:8,4,2) #1~8 4列2行的矩陣
v1 <-c("Marry","Jone")
L1 <- list(mt1=m1,ve1=v1)
L1
L1[[1]] #list要取元素 要2個中括號
L1[[2]]
t1<- L1[1]
t2 <-L1[[1]]
e1<-lm(Murder~UrbanPop,data=USArrests)
e1[[1]] #取係數
e1[[2]] #取殘差
L1$mt1
e1$coefficients

## Data mining

### 支持度，信心度

setwd("D:\\DM")
ti<-read.csv("titanic-train.csv")
str(ti)

library(infotheo)

ti$Age<-discretize(ti$Age, "equalwidth", 3)
ti$Age<-as.factor(ti$Age$X)

ti$Ticket<-discretize(ti$Ticket, "equalfreq", 3)
ti$Ticket<-as.factor(ti$Ticket$X)

ti$Pclass<-as.factor(ti$Pclass)
ti$SibSp<-as.factor(ti$SibSp)
ti$Parch<-as.factor(ti$Parch)
ti$Fare<-as.factor(ti$Fare)
str(ti)

require(arules)
rule<-apriori(ti, parameter = list(supp=0.1,conf=0.7),
appearance=list(rhs=c("Survived=Yes","Survived=No")))
sort.rule<-sort(rule, by="support")
subset.matrix<-as.matrix(is.subset(x=sort.rule, y=sort.rule))
subset.matrix[lower.tri(subset.matrix, diag=T)]<-NA
redundant<-colSums(subset.matrix, na.rm=T)>=1
sort.rule<-sort.rule[!redundant]
sort.rule<-as(sort.rule, "data.frame")
write.csv(sort.rule, "ti-rules.csv")

### 正確率
setwd("C:\\Users\\lydia\\Downloads\\DM")
titrain<-read.csv("#5-titanic-train.csv")
titest<-read.csv("#5-titanic-test.csv")

set.seed(123)
n <- nrow(titrain)
n <- nrow(titest)
sindex <- sample(n,round(n*0.7))

library(RWeka)

ctree<- J48(Survived~.,data = titrain, control = Weka_control(M=9,C=0.25))
print(ctree)

library(partykit)
rparty.tree<-as.party(ctree)
plot(rparty.tree)

titest_predicted=predict(ctree, titest, type="class")
titest_predicted
titest$predict=titest_predicted

cm<-table(titest$Survived,titest$predict,dnn=c("實際","預測"))
cm

accuracy<-sum(diag(cm))/sum(cm)
accuracy

### 訓練資料集/測試資料集/正確率
setwd("C:\\Users\\lydia\\Downloads\\DM")
titrain<-read.csv("#5-titanic-train.csv")
titest<-read.csv("#5-titanic-test.csv")
str(titest)
str(titrain)

titest$Pclass<-factor(titest$Pclass)
titest$SibSp<-factor(titest$SibSp)
titest$Parch<-factor(titest$Parch)
titrain$Pclass<-factor(titrain$Pclass)
titrain$SibSp<-factor(titrain$SibSp)
titrain$Parch<-factor(titrain$Parch)

set.seed(123)
n <- nrow(titrain)
n <- nrow(titest)
sindex <- sample(n,round(n*0.7))

library(rpart)

ti_clf<- rpart(Survived~.,data = titrain, method="class")

library(rattle)
fancyRpartPlot(ti_clf)

library(rpart.plot)
prp(ti_clf, faclen = 0, fallen.leaves = TRUE, shadow.col = "gray")

library(partykit)
rparty.tree<-as.party(ti_clf)
plot(rparty.tree)

titest_predicted=predict(ti_clf, titest, type="class")
titest_predicted

cm<-table(titest$Survived,titest_predicted,dnn=c("實際","預測"))
cm

accuracy<-sum(diag(cm))/sum(cm)
accuracy

printcp(ti_clf)
plotcp(ti_clf)
prunetree_ti<-prune(ti_clf,
cp = ti_clf$cptable[which.min(ti_clf$cptable[,"xerror"]),"CP"])

prunetree_pred<- predict(prunetree_ti, newdata = titest, type="class")
table(real=titest$Survived, predict=prunetree_pred)

prunetree_cm <- table(real=titest$Survived, predict = prunetree_pred)
prunetree_cm
sum(diag(prunetree_cm))/sum(prunetree_cm)

library(caret)
library(e1071)
train_control <- trainControl(method="cv", number = 10)
train_control.model<-train(Survived~., data = titrain, method="rpart", trcontrol=train_control)
train_control.model