#兩樣本問題t檢驗、方差齊次性檢驗、分布檢驗
#相關與回歸分析
#例一某種元件的壽命X（小時），服從正態分布，N（mu,sigma2），其中mu,sigma2均未知，16隻元件的壽命如下：問是否有理由認為元件的平均壽命大於225小時。
#分析：按題意，需檢驗
#H0： μ ≤ 225 H1: μ > 225
#此問題屬於單邊檢驗問題，可以使用R語言t.test
#t.test(x,y=NULL,
alternative=c(「two.sided」,「less」,「greater」),mu=0,paired=FALSE,var.equal=FALSE,conf.level=0.95)#其中x,y構成兩向量，（如果只提供x,則作單個正態總體的均值檢驗，如果提供x,y則作兩個總體的均值檢驗)，
#alternative表示備擇假設，two.sided(預設)，雙邊檢驗(H1:μ≠H0),less表示單邊檢驗(H1:μ<μ0),greater表示單邊檢驗(H1:μ>μ0)，
#mu表示原假設μ0，paired表示是否配對樣本T檢驗，conf.level置信水平，即1-α，通常是0.95，var.equal是邏輯變量，var.equal=TRUE表示兩樣品方差相同，var.equal=FALSE（預設）表示兩樣本方差不同。
X<-c(159, 280, 101, 212, 224, 379, 179, 264,222, 362, 168, 250, 149, 260, 485, 170)
t.test(X, alternative = 「greater」, mu = 225)
#可見P值為0.257 > 0.05 ，不能拒絕原假設，接受H0，即平均壽命不大於225小時。
#例二 在平爐上進行的一項試驗以確定改變操作方法的建議是否會增加鋼的得率，試驗時在同一個平爐上進行的，每煉一爐剛時除操作方法外，其它條件都儘可能做到相同，先用標準方法煉一爐，然後用新方法煉一爐，以後交替進行，各煉了10爐，其得率分別為
#標準方法 78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3
#新方法 79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1
#設這兩個樣本相互獨立，且分別來自正態總體N（μ1, σ2)和N（μ2, σ2)，其中μ1，μ2和σ2未知。問新的操作能否提高得率？（取α=0.05）
#分析：根據題意，需要假設
#H0： μ1 ≥ μ2 H1: μ1 < μ2
#這裡假定σ12=σ22=σ2，因此選擇t.test，var.equal=TRUE
#R代碼：
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
t.test(X,Y,var.equal = TRUE,alternative = 「less」)
#可見P值<0.05，拒絕原假設，即新的操作能夠提高得率。
#因為數據是成對出現的，所以採用成對數據t檢驗比上述的雙樣本均值檢驗更準確。
#所謂成對t檢驗就是Zi=Xi-Yi，再對Z進行單樣本均值檢驗（與0比較）
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
t.test(X-Y, alternative = 「less」)
t.test(X,Y,paired = TRUE,alternative = 「less」)
#例二方差相同嗎？檢驗方差是否相同
#H0： σ12 = σ22 H1: σ12 ≠ σ22
#方差檢驗可以用var.test
#var.test(x, y, ratio = 1,
alternative = c(「two.sided」, 「less」, 「greater」),conf.level = 0.95, …)#x,y是來自兩樣本數據構成的向量，ratio是方差比的原假設，預設值為1.alternative是備擇假設，two.sided表示雙邊檢驗,greater表示單邊檢驗
#R代碼：
X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
var.test(X,Y)
#可見P值為0.559>0.05，不拒絕原假設，認為兩者方差相同
#二項分布的參數檢驗
#有一批蔬菜種子的平均發芽率p0=0.85，現隨即抽取500粒，用種衣劑進行浸種處理，結果有445粒發芽。試檢驗種衣劑對種子發芽率有無效果。
#根據題意，所檢驗的問題為
#H0：p=p0=0.85， H1：p≠p0
#可以用R語言的binom.test
#binom.test(x, n, p = 0.5,
alternative = c(「two.sided」, 「less」, 「greater」),conf.level = 0.95)#其中x是成功的次數；或是一個由成功數和失敗數組成的二維向量。n是試驗總數，當x是二維向量時，此值無效。P是原假設的概率。
#R語言代碼：
binom.test(445,500,p=0.85)
#可知P值0.01207<0.05，拒絕原假設，說明種衣劑對種子的發芽率有顯著效果。
#擬合優度檢驗
#例：某體育網站為了確定球迷對5支具備奪取世界盃冠軍球隊的支持情況，隨即抽取了1000名球迷進行問卷調查：每個人選一支最支持的隊伍。這5支球隊分別記為A、B、C、D、E。下表是根據樣本資料整理的5支球隊支持人數的頻數分布。
#試根據這些數據判斷球迷對五支球隊的支持數量有無明顯差異？
#分析：
#如果球迷對5支球隊的支持無顯著差異，那麼，就可以認為喜好這5支球隊的人呈均勻分布，即支持每隊的人數各佔20%。據此假設
#H0：支持5支球隊的人數服從分布均勻
#可以使用Pearson χ2擬合優度檢驗，R語言中調用chisq.test(X)
#chisq.test(x, y = NULL, correct = TRUE,
p = rep(1/length(x), length(x)), rescale.p = FALSE,simulate.p.value = FALSE, B = 2000)#其中x是由觀測數據構成的向量或者矩陣，y是數據向量（當x為矩陣時，y無效)。correct是邏輯變量，標明是否用於連續修正，TRUE(預設值)表示修正，FALSE表示不修正。p是原假設落在小區間的理論概率，預設值表示均勻分布，rescale.p是邏輯變量，選擇FALSE（預設值）時，要求輸入的p滿足和等於1；選擇TRUE時，並不要求這一點，程序將重新計算p值。simulate.p.value邏輯變量（預設值為FALSE），當為TRUE，將用仿真的方法計算p值，此時,B表示仿真的此值。
#R語言代碼：
X<-c(210, 312, 170, 85, 223)
chisq.test(X)
#p值非常接近於0，拒絕原假設
#泊松分布檢驗
#為研究電話總機在某段時間內接到的呼叫次數是否服從Poisson分布，現收集了42個數據，如下表所示，通過對數據的分析，
#問能否確認在某段時間內接到的呼叫次數服從Poisson分布（α = 0.1）？
#接到呼叫次數 0 1 2 3 4 5 6
#頻數 7 10 12 8 3 2 0
X<-0:6; Y<-c(7, 10, 12, 8, 3, 2, 0)
#計算理論分布，其中mean(rep(X,Y))為樣本均值
q<-ppois(X, mean(rep(X,Y))); n<-length(Y) #生成<=x的累積概率密度
p=rep(0,n)
p[1]<-q[1]; p[n]<-1-q[n-1]
for (i in 2:(n-1))
p[i]<-q[i]-q[i-1]
#作檢驗
chisq.test(Y,p=p)
#提示結果可能不準確，因為皮爾森卡方擬合由度檢驗要求分組後每組的頻數至少要大於等於5，而後三組中出現的頻率分別為3，2，0，均小於5，
#解決問題的方法是將後三組合成一組，此時的頻數為5，滿足要求，重寫R語言代碼
#輸入數據
X<-0:6; Y<-c(7, 10, 12, 8, 3, 2, 0)
#計算理論分布，其中mean(rep(X,Y))為樣本均值
q<-ppois(X, mean(rep(X,Y))); n<-length(Y)
p<-rep(0,n)
p[1]<-q[1]; p[n]<-1-q[n-1]
for (i in 2:(n-1))
p[i]<-q[i]-q[i-1]
#重新分組
Z<-c(7, 10, 12, 8, 5)
#重新計算理論分布
n<-length(Z); p<-p[1:n-1]; p[n]<-1-q[n-1]
#作檢驗
chisq.test(Z,p=p)
#p值大於0.1，不拒絕原假設，表明服從泊松分布
#正態分布的檢驗
#已知15名學生體重如下，問是否服從正態分布
#75.0, 64.0, 47.4, 66.9, 62.2, 62.2, 58.7, 63.5,
#66.6, 64.0, 57.0, 69.0, 56.9, 50.0, 72.0
w <- c(75.0, 64.0, 47.4, 66.9, 62.2, 62.2, 58.7, 63.5,
66.6, 64.0, 57.0, 69.0, 56.9, 50.0, 72.0)
shapiro.test(w)
#P值>0.05，接受原假設，認為來自正態分布總體。
#也可以用ks.test檢驗正態性
#統計學裡, Kolmogorov–Smirnov 檢驗(亦稱：K–S 檢驗)是用來檢驗數據是否符合某種分布的一種非參數檢驗，通過比較一個頻率分布f(x)與理論分布g(x)或者兩個觀測值分布來判斷是否符合檢驗假設。其原假設H0:兩個數據分布一致或者數據符合理論分布。拒絕域構造為：D=max| f(x)- g(x)|，當實際觀測值D>D(n,α)則拒絕H0，否則則接受H0假設。由於KS檢驗不需要知道數據的分布情況，在小樣本的統計分析中效果比較好。（大樣本數據下，使用t-檢驗；小樣本數據，使用t-檢驗會出現較大的偏差）
#ks.test(x, y, …,
alternative = c(「two.sided」, 「less」, 「greater」),exact = NULL)#x:觀測值向量
#y:第二觀測值向量或者累計分布函數，如pnorm(正態分布函數，一般做正態檢測的時候直接輸入pnorm)，只對連續CDF有效
#alternative = c(「two.sided」, 「less」, 「greater」):雙側檢驗還是單側檢驗
#exact：默認為NULL，也可以是其他邏輯值，表明是否需要計算精確的P值
#結果解釋：D：D值越小，越接近0，表示樣本數據越接近正態分布(簡單來說，D越小越好),p：p-value小於顯著性水平α(0.05)，則拒絕H0（p越大越好）
x1=c(87,77,92,68,80,78,84,77,81,80,80,77,92,86,76,80,81,75,77,72,81,90,84,86,80,68,77,87,76,77,78,92,75,80,78)
x1<-scale(x1)#標準化
ks.test(x1,pnorm)
#注意事項：
#在做單樣本K-S檢驗或者正態檢驗時，有時會有錯誤提示「Kolmogorov - Smirnov檢驗裡不應該有連結」，
#這是因為K-S檢驗只對連續CDF有效，而連續CDF中出現相同值的概率為0，也就是說數據中倘出現相同值，則連續分布的假設不成立，因此R會報錯。這也提醒我們，在做正態性檢驗之前，要先對數據進行描述性分析，對數據整體要先有個大致的認識，這也才後續才能選擇正確的檢驗方法。
#檢驗兩個未知的分布是否為同一分布
X<-c(0.61, 0.29, 0.06, 0.59, -1.73, -0.74, 0.51, -0.56,1.64, 0.05, -0.06, 0.64, -0.82, 0.37, 1.77,
2.36, 1.31, 1.05, -0.32, -0.40, 1.06, -2.47,0.39, 1.09, -1.28)
Y<-c(2.20, 1.66, 1.38, 0.20, 0.36, 0.00,0.96, 1.56, 0.44, 1.50, -0.30, 0.66, 2.31, 3.29, -0.27, -0.37, 0.38, 0.70,0.52, -0.71)
ks.test(X,Y)
#p=0.5286>0.05,不能拒絕原假設，兩個服從相同的分布
#注意：如果樣本太少，可能出現錯誤結論
a1=c(-1020,-1040);a2=c(1020,1040)
t.test(a1,a2)
#一般情況，樣本越小，往往容易得到不能拒絕原假設的結論
#pearson相關係數
a=c(1,3,5,7,9);b=c(1,4,6,9,10)
cor(a,b)
cor.test(a,b) #檢驗相關係數的顯著性
cor(iris[1:4]) #相關係數,參數填數據集，則計算相關係數矩陣
a1=rnorm(5);b1=rnorm(5);cor(a1,b1);cor.test(a1,b1) #自己模擬生成兩個變量
#spearman相關係數,亦即秩相關係數
#spearman和kendall都是等級相關係數，亦即其值與兩個相關變量的具體值無關，而僅僅與其值之間的大小關係有關。
#spearman相關係數，亦即秩相關係數，根據隨機變量的等級而不是其原始值衡量相關性的一種方法。
m=c(1,2,4,3);n=c(100,101,102,103)
m1=c(30,31,35,34);n1=c(85,87,90,93)
cor(m,n);cor(m1,n1)
cor(m,n,method = 「spearman」);cor(m1,n1,method = 「spearman」)
cor.test(m,n,method = 「spearman」);cor.test(m1,n1,method = 「spearman」)
#spearman相關係數的計算可以由計算pearson係數的方法,只需要把原隨機變量中的原始數據替換成其在隨機變量中的等級順序即可:
acf #自相關和協方差函數
acf(airmiles,type=『correlation』,lag.max=10) #自相關
pacf(airmiles,lag.max=10) #偏自相關
pairs(~Sepal.Length+Sepal.Width+Petal.Length+Petal.Width,data=iris,
main=「Simple Scatterplot Matrix」) #散點圖矩陣
install.packages(「scatterplot3d」) #3D散點圖
library(scatterplot3d)
scatterplot3d(irisSepal.Length,irisSepal.Length, irisSepal.Length,irisPetal.Length, iris$Petal.Width)
install.packages(「corrgram」) #有興趣的同學自己練習
library(corrgram)
#1、設置排序處理
corrgram(mtcars,order=TRUE)
#2、設置上下三角面板形狀
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie)
#3、只顯示下三角部分
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=NULL)
#4、調整面板顏色
corrgram(mtcars,order=TRUE,lower.panel=panel.shade,upper.panel=panel.pie,
col.regions=colorRampPalette(c(「darkgoldenrod4」,「burlywood1」,「white」, 「darkkhaki」,「darkgreen」)))
install.packages(「corrplot」)
library(corrplot)
#1、使用不同的method繪製相關矩陣圖
methods<-c(「circle」,「square」,「ellipse」,「pie」,「shade」,「color」)
par(mfrow=c(2,3))
t0=mapply(function(x){corrplot(cor(mtcars), method=x,order=「AOE」)},methods)
par(mfrow=c(1,1))
#2、設置method=color繪製熱力矩陣圖
corrplot(cor(mtcars), method=「color」, order = 「AOE」,tl.col=「black」,tl.srt=45,
addCoef.col=「black」,col=colorRampPalette(c("#7F0000",「red」,"#FF7F00",
「yellow」,「white」, 「cyan」, 「#007FFF」, 「blue」,"#00007F"))(20))
#3、繪製上下三角及不同色彩的相關矩陣圖
library(RColorBrewer)
par(mfrow=c(2,2))
corrplot(cor(mtcars),type=「lower」)
corrplot(cor(mtcars),type=「lower」,order=「hclust」,
col=brewer.pal(n=8,name=「RdYlBu」))
corrplot(cor(mtcars),type=「upper」,order=「AOE」,
col=c(「black」,「white」),bg=「lightblue」)
corrplot(cor(mtcars),type=「upper」,order=「FPC」,
col=brewer.pal(n=8, name=「PuOr」))
par(mfrow=c(1,1))
d<-sqrt(1-cor(mtcars)^2)
hc<-hclust(as.dist(d))
plot(hc)
rect.hclust(hc,k=3)
install.packages(「pvclust」)
library(pvclust)
cluster.bootstrap <- pvclust(mtcars, nboot=1000, method.dist=「correlation」)
plot(cluster.bootstrap)
pvrect(cluster.bootstrap) #自己練習部分結束
#回歸分析
#回歸分析的形式：fitted.model=lm(formula,data=data.frame)
#yx;yx+1 表示y對x的一元線性回歸，默認包含常數項
y1=rnorm(30);x1=1:30
y2=3*x2-5;
x2=rnorm(30);k2=lm(y2~x2);summary(k2)
k1=lm(y1~x1)
summary(k1)
#log(y)~x1+x2 表示y的對數對自變量的回歸，默認帶常數項
#yploy(x,2),y1+x+I(x^2) 對x的一元二次多項式回歸
k3=lm(y1~1+x1+I(x1^2));k3;summary(k3)
#y~A， 表示單因素方差分析，A是分類變量
#y~A+x 單因素協方差分析，A分類變量，x連續型變量
install.packages(「readxl」) #讀取Excel數據的包
library(readxl)
y=c(16.3,16.8,19.2,18, 19.5,20.9,21.1,20.9,20.3,22,23.7,34.2,36.2)
x1=c(1.1, 1.4, 1.7, 1.7, 1.8, 1.8, 1.9, 2, 2.3 ,2.4,2.5,2.6,2.7)
x2=c(1.1, 1.5, 1.8, 1.7, 1.9 ,1.8, 1.8 ,2.1, 2.4, 2.5,2.7,2.9,2.94)
A=c(0, 1 ,1, 0, 0 ,1, 1, 1, 0, 0,2,2,2)
A=as.factor(A) #注意R語言中分類變量的處理
B=c(「wuhan」,「changsha」,「wuhan」,「hefei」,「changsha」,「wuhan」,「hefei」,「changsha」,「wuhan」,「hefei」,「wuhan」,「changsha」,「hefei」)
m1=lm(y~x1+x2+A)
m1
m2=aov(y~A) #單因素方差分析
m2
summary(m2)
m3=lm(y~x1+I(x1**2)) #多項式要用I括起來表示整體
m4=lm(y~x1+B)
m4
m5=lm(y~x1+B+x1:B);m5 #交互項用：
m6=lm(y~.,data=) #數據集中除因變量的所有其它變量作為自變量
m7=lm(y~x1x2B);m7 #三變量單獨、兩兩交互、三者交互的組合
m8=lm(y~(x1+x2+B)**2);m8 #單獨、兩兩交互
#summary(詳細結果)coefficients(模型參數)confint(置信區間)
#fitted(擬合值)residuals(殘差值)AIC(赤池信息)predict(預測值)
m9=lm(y~x1+x2)
summary(m9)
coefficients(m9)
confint(m9)
AIC(m9)
fitted(m9)
residuals(m9)
predict(m9)
formula(m3)
plot(x1,y) #散點圖
boxplot(x1~B) #箱線圖
plot(m3) #默認畫4個圖，殘差對擬合值圖，殘差正態概率圖，標準化殘差絕對值平方根對擬合值圖，cook距離圖
k1=par(mfrow=c(1,1)) #圖的展示方式，2行2列
plot(m3)
k2=lm(y~x1)
plot(y,x1)
abline(lm(y~x1)) #散點圖加擬合線
#回歸模型診斷
install.packages(「car」) #很多回歸診斷功能
library(car)
library(carData)
durbinWatsonTest(m3) #自相關檢驗
crPlots(m3) #看線性趨勢
ncvTest(m3) #原假設是同方差
vif(m3) #方差膨脹因子
outlierTest(m3) #原假設無離群點