#在R執行SAS的glm lsmeans
#使用的資料為R內建的dataset- airquality
#可以先看一下該資料的描述(紐約某一年的空氣品質指標數據)
?airquality
head(airquality)
#輸出為csv檔供SAS使用
airquality$Month<- as.factor(airquality$Month)
write.csv(airquality, file = "airquality.csv")
- 10月 06 週六 201819:15
在R中執行 SAS的glm lsmeans
#在R執行SAS的glm lsmeans
#使用的資料為R內建的dataset- airquality
#可以先看一下該資料的描述(紐約某一年的空氣品質指標數據)
?airquality
head(airquality)
#輸出為csv檔供SAS使用
airquality$Month<- as.factor(airquality$Month)
write.csv(airquality, file = "airquality.csv")
- 9月 06 週四 201808:36
資料處理技巧(2)_R語言
#資料處理技巧(2)
install.packages(c("magrittr","tidyr","dplyr"))
install.packages("tidyverse") #上面三個包都在其中
library(tidyverse)
cars %>% summary() #把cars叫進summary()。等同summary(cars)
#應用運算符號%>%
birth<- 1995
age<- Sys.Date() %>%
format(format= "%Y") %>%
as.numeric() %>%
`-` (birth)
#調整輸入位置
cars_lm<- lm(formula = dist~ speed, data= cars) #傳統方式
cars_lm<- cars %>%
lm(formula = dist~ speed, data = .) #以. 指定資料輸入的位置
#gather() 將多個數值變數堆積在同一個數值變數中(value)，再用一個類別變數(key)紀錄數值變數的來源。
team_name<- c("Bull", "Warrior")
wins<- c(72,73)
losses<- c(10,9)
team<- data.frame(team_name,wins,losses)
team
gather(team, key = variable, value = values, wins, losses)
install.packages(c("magrittr","tidyr","dplyr"))
install.packages("tidyverse") #上面三個包都在其中
library(tidyverse)
cars %>% summary() #把cars叫進summary()。等同summary(cars)
#應用運算符號%>%
birth<- 1995
age<- Sys.Date() %>%
format(format= "%Y") %>%
as.numeric() %>%
`-` (birth)
#調整輸入位置
cars_lm<- lm(formula = dist~ speed, data= cars) #傳統方式
cars_lm<- cars %>%
lm(formula = dist~ speed, data = .) #以. 指定資料輸入的位置
#gather() 將多個數值變數堆積在同一個數值變數中(value)，再用一個類別變數(key)紀錄數值變數的來源。
team_name<- c("Bull", "Warrior")
wins<- c(72,73)
losses<- c(10,9)
team<- data.frame(team_name,wins,losses)
team
gather(team, key = variable, value = values, wins, losses)
- 8月 30 週四 201818:20
基因富集分析 (gene set enrichment set analysis)_GO/ GSEM/ KEGG_R語言
#基因富集分析 (gene set enrichment set analysis)
#數據庫下載: 與自己的差異基因進行搜尋及比對
#AnnotationHub是生物數據庫的中轉站，方便搜尋目標數據，另一個相似套件為biomaRt
#參考網址:https://www.jianshu.com/p/ae94178918bc
source("https://bioconductor.org/biocLite.R")
install.packages("BiocManager")
BiocManager::install(c("AnnotationHub","BiocGenerics","parallel"))
library(AnnotationHub)
ah <- AnnotationHub()
ah
- 8月 28 週二 201817:00
資料探索的檢定_主成分分析/多變量變異數分析/集群分析_R語言
#資料探索的檢定
#資料型態: 幾乎任何資料都可以拿來做資料探索，實務上，當變數數目跟觀測數都夠多時，資料探索的用處比較大。資料探索的目的為處理觀測值之間的交互關係，並將期間模式凸顯給實驗者。
#主成分分析和因素分析 (Principle component analysis (PCA) and factor analysis): 都是藉由加權變數來將每個個體的差異最大化。許多方面，主成分分析跟相關係數和回歸分析很類似。只要資料中每個個體都有兩個以上的觀測值就可使用。主成分分析假第資料是連續的且遵從常態分佈。若執行此分析的目的是為了建立新假設就可以忽略此假定。
#視覺化: 藉由「每個個體指有兩個觀測值」的特例將主成分分析的運作視覺化。想像將兩個相關的變數製成散佈圖，其資料形成橢圓形雲狀物。PCA會計算通過此橢圓形長軸的直線，並將此直線當成第一個主要軸線(PCA1)。通過資料雲，且和第一個軸線垂直的直線就是PCA2。分析過程是電腦以多維空間的形式進行，每個變數都代表一個維度，每條通過資料雲的直線都來自對於個變數做適當的加權。
- 8月 27 週一 201816:10
兩個變數之間有因果關係嗎_回歸_R語言
#兩個變數之間有有因果關係嗎
#許多情況下，一組觀測值的取值很明顯是取決於另一組觀測值。本次的案例中每個「個體」都有兩個觀測值。一個觀測值是「原因變數」、「x變數」、「預測變數」、「自變數」，此變數的取值為實驗者的設定或選擇；另一個觀測值是「效果變數」、「y變數」、「應變數」，此變數的取值非實驗者設定。有一系列方式可以判斷原因和效果之間的關係形式和強度，每個方法對變數及其關係的假定各有不同，這裡考慮五種檢定: 線性回歸、Kendall最佳配適線、羅吉斯回歸、第二型模式回歸和多項式回歸。
- 8月 26 週日 201817:00
兩個變數之間有相關或相關性嗎?_R語言
#檢定關係的檢定
#兩個變數之間有相關或相關性嗎?
#將資料中每次的觀測都賦予定性的值(如殼蛋的重量分級為"LL","L","M","MS","S")
#將觀測結果分類到各類別，即可運用檢定判斷觀測結果的類別是否獨立，可用的檢定有"卡方相關性檢定"、"phi係數"和"Cramér係數"
#卡方相關性檢定 (Chi-square test of association): 觀測結果可被分配兩個變數各自的組別中。本檢定並未對資料的形式限制，是一種無母數檢定，但卻沒有對應的母數分析。
#本檢定統計量為各類別組合的預期觀察數和實際觀測數的差值平方後，再全部相加。將此數值查卡方表，自由度為欄位數減一乘以列數減一。
#須注意 (1)預期觀察數不可小於1，應將各類別合併避免此問題 (2)預期觀測數小於5的比率應超過20%，若沒有提升到5以上的辦法，可改用Fisher精確性檢定。
#卡方檢定一定要用在頻率 (觀測次數)，不能用百分比或是轉換過的資料。統計結果會輸出統計量，自由度和P值的相關資訊。
#Example: 調查雞隻打架頻率是否和品種有關。在種雞場中隨機挑選了50隻雞作為定性調查。打架頻率被分為"低(L)"和"高(H)"，雞隻被分為四個品種(類別): "絲羽烏骨雞(Silkies)"、"竹崎(Shek-Ki)"、"花東(Hua-Tung)"、"金門(Quemoy)"。
chicken<- matrix(c(10,8,2,7,2,6,10,5),nrow=2,byrow=T, dimnames = list(c(1:2),c("Silkies","Shek-Ki","Hua-Tung","Quemoy")))
#虛無假設為兩種打架頻率的雞隻類別分配比例相等，意旨每個品種和某一打架頻率共同出現的機率相等。該檢定無法判定哪些類別組合比預期出現多還是少，不過可以直接檢視原始資料跟期望值大概就可得知。
chisq.test(chicken) #虛無假設被拒絕，回頭看原始資料可以發現Hua-Tung(花東)較常和打鬥頻率高有關。
#兩個變數之間有相關或相關性嗎?
#將資料中每次的觀測都賦予定性的值(如殼蛋的重量分級為"LL","L","M","MS","S")
#將觀測結果分類到各類別，即可運用檢定判斷觀測結果的類別是否獨立，可用的檢定有"卡方相關性檢定"、"phi係數"和"Cramér係數"
#卡方相關性檢定 (Chi-square test of association): 觀測結果可被分配兩個變數各自的組別中。本檢定並未對資料的形式限制，是一種無母數檢定，但卻沒有對應的母數分析。
#本檢定統計量為各類別組合的預期觀察數和實際觀測數的差值平方後，再全部相加。將此數值查卡方表，自由度為欄位數減一乘以列數減一。
#須注意 (1)預期觀察數不可小於1，應將各類別合併避免此問題 (2)預期觀測數小於5的比率應超過20%，若沒有提升到5以上的辦法，可改用Fisher精確性檢定。
#卡方檢定一定要用在頻率 (觀測次數)，不能用百分比或是轉換過的資料。統計結果會輸出統計量，自由度和P值的相關資訊。
#Example: 調查雞隻打架頻率是否和品種有關。在種雞場中隨機挑選了50隻雞作為定性調查。打架頻率被分為"低(L)"和"高(H)"，雞隻被分為四個品種(類別): "絲羽烏骨雞(Silkies)"、"竹崎(Shek-Ki)"、"花東(Hua-Tung)"、"金門(Quemoy)"。
chicken<- matrix(c(10,8,2,7,2,6,10,5),nrow=2,byrow=T, dimnames = list(c(1:2),c("Silkies","Shek-Ki","Hua-Tung","Quemoy")))
#虛無假設為兩種打架頻率的雞隻類別分配比例相等，意旨每個品種和某一打架頻率共同出現的機率相等。該檢定無法判定哪些類別組合比預期出現多還是少，不過可以直接檢視原始資料跟期望值大概就可得知。
chisq.test(chicken) #虛無假設被拒絕，回頭看原始資料可以發現Hua-Tung(花東)較常和打鬥頻率高有關。
- 8月 24 週五 201818:03
資料有大於兩種互相獨立的分組方法_R語言
#資料有大於兩種互相獨立的分組方法
#多因子檢定 (Multifatorial testing): 需分析因子數越多，需要分析的交互作用項也會快速增加。
#三因子變異數分析 (無重複試驗) [Three-way ANOVA (without replication)]
#Exzample: 白羅曼鵝的性別、受日照長短、試驗區域對其採食量影響
intake<- c(78.1,76.3,
69.5,73.2,
82.4,83.0,
72.3,70.0)
sex<- rep(c("F","F","M","M"),times= 2)
light<- c("L","L","L","L","S","S","S","S")
region<-c("A","B","A","B",
"A","B","A","B")
data.frame(intake=intake, sex=sex, light=light,region=region)
summary(aov(intake~sex*light*region))
summary(aov(intake~sex*light*region-sex:light:region)) #為無重複試驗時須要扣除同時有三因子之交互作用 (-sex:light:region)
#多因子檢定 (Multifatorial testing): 需分析因子數越多，需要分析的交互作用項也會快速增加。
#三因子變異數分析 (無重複試驗) [Three-way ANOVA (without replication)]
#Exzample: 白羅曼鵝的性別、受日照長短、試驗區域對其採食量影響
intake<- c(78.1,76.3,
69.5,73.2,
82.4,83.0,
72.3,70.0)
sex<- rep(c("F","F","M","M"),times= 2)
light<- c("L","L","L","L","S","S","S","S")
region<-c("A","B","A","B",
"A","B","A","B")
data.frame(intake=intake, sex=sex, light=light,region=region)
summary(aov(intake~sex*light*region))
summary(aov(intake~sex*light*region-sex:light:region)) #為無重複試驗時須要扣除同時有三因子之交互作用 (-sex:light:region)
- 8月 23 週四 201817:00
以R語言解答動物科學統計習題- 習題七 屬性資料的分析
#屬性資料的分析
#擲兩個骰子，兩個骰子的點數一樣的機率為多少?兩個骰子總共點數不足5的機率為多少
dbinom(2, size = 2, prob = 1/6)*6
first<- (1:6) ;second<- (1:6)
results<- c()
i<-1
while (length(results) < 100000) {
results[i] <- (sample(first , 1)+ sample(second,1)) <5
i<-i+1
}
sum(results)/100000
#擲兩個骰子，兩個骰子的點數一樣的機率為多少?兩個骰子總共點數不足5的機率為多少
dbinom(2, size = 2, prob = 1/6)*6
first<- (1:6) ;second<- (1:6)
results<- c()
i<-1
while (length(results) < 100000) {
results[i] <- (sample(first , 1)+ sample(second,1)) <5
i<-i+1
}
sum(results)/100000
- 8月 21 週二 201817:00
R語言 資料處理技巧(1)
#資料處理技巧
#解構資料框
#解構為單一變數值: 取出內建資料iris 第一筆觀測值得第五個變數 (Species)的值
iris[1,5]
iris[1,"Species"]
#解構資料框
#解構為單一變數值: 取出內建資料iris 第一筆觀測值得第五個變數 (Species)的值
iris[1,5]
iris[1,"Species"]
- 8月 20 週一 201817:00
R語言 檢定多於兩組(多組)的觀察值之間是否有差異
#檢定多於兩組(多組)的觀察值之間是否有差異
#重複測量(Repeated measures): Friedman檢定、重複測量變異數分析
#Friedman檢定 (重複測量): 二因子變異數分析對應的無母數檢定。對資料分布沒有任何假設，只要求資料是可排序的。只有在每個因子值組合僅有一個觀察值時適用。
#Example: 選擇四天分別取出六個池塘的水樣計算藍藻數(個/mm3)，每次取樣只取一個樣本，就可使用該檢定。
pondcells<- matrix(c(130,125,350,375,225,235,
115,120,375,200,250,200,
145,170,235,275,225,155,
200,230,140,325,275,215),
nrow=4,byrow = T,
dimnames = list(1:4,c("A","B","C","D","E","F"))) #此檢定要將資料整理成矩陣形式
friedman.test(pondcells) #自由度為池塘數減一
#重複測量變異數分析 (Repeated-measures ANOVA): 此種試驗設計下，每個因子值組合都只有一個才可以執行二因子變異數分析。在重複測量的設計下，影響因子需互相獨立或是得大量降低自由度。
pond<-c("A","A","A","A",
"B","B","B","B",
"C","C","C","C",
"D","D","D","D",
"E","E","E","E",
"F","F","F","F")
cells<-c(130,115,145,200,
125,120,170,230,
350,375,235,140,
375,200,275,325,
225,250,225,275,
235,200,155,215)
day<- c(rep(1:4,6))
summary(aov(cells~pond*day + Error(cells/(pond*day))))
#重複測量(Repeated measures): Friedman檢定、重複測量變異數分析
#Friedman檢定 (重複測量): 二因子變異數分析對應的無母數檢定。對資料分布沒有任何假設，只要求資料是可排序的。只有在每個因子值組合僅有一個觀察值時適用。
#Example: 選擇四天分別取出六個池塘的水樣計算藍藻數(個/mm3)，每次取樣只取一個樣本，就可使用該檢定。
pondcells<- matrix(c(130,125,350,375,225,235,
115,120,375,200,250,200,
145,170,235,275,225,155,
200,230,140,325,275,215),
nrow=4,byrow = T,
dimnames = list(1:4,c("A","B","C","D","E","F"))) #此檢定要將資料整理成矩陣形式
friedman.test(pondcells) #自由度為池塘數減一
#重複測量變異數分析 (Repeated-measures ANOVA): 此種試驗設計下，每個因子值組合都只有一個才可以執行二因子變異數分析。在重複測量的設計下，影響因子需互相獨立或是得大量降低自由度。
pond<-c("A","A","A","A",
"B","B","B","B",
"C","C","C","C",
"D","D","D","D",
"E","E","E","E",
"F","F","F","F")
cells<-c(130,115,145,200,
125,120,170,230,
350,375,235,140,
375,200,275,325,
225,250,225,275,
235,200,155,215)
day<- c(rep(1:4,6))
summary(aov(cells~pond*day + Error(cells/(pond*day))))
生物資訊 (4)