R INTRODUCTION
Biểu đồ trong R
So sánh 2 nhóm
2 biến phân nhóm (phân loại 2 giá trị)
So sánh hai nhóm: Ratio risk
Biểu đồ tương quan đa biến
Tương quan 2 biến
plot(fnbmd~age,pch=16,xlab="Age",ylab="FNBMD", xlim=c(0,1.2), col="blue")
Tương quan giữa nhiều biến
require psych
pairs.panels(dat)
Biên tập dữ liệu
data=data.frame(id,x,y)
Xem xét cấu trúc dữ liệu fix(dat), names(dat), head(dat). na.strings="."
Tạo ra biến mới
data$sum=data$x1+data$x2 (sum không có từ trước)
Biến điều kiện
data$sex[gender=="male"]<-1 (sex không có trước, gender có male, female)
data$group[id<4]<-"A", data$group[id>=4]<-"B"
Hoán chuyển số- text
as.numeric(), as.character(), as.vector(), as.matric(), as.data.frame()
Trích xuất dữ liệu
data2 = subset(data,id<=8&x<30)
Sắp xếp dữ liệu
data2 = data[order(x), ]
Gộp dữ liệu
dat=merge(dat1,dat2,by="id",all.x=T,all.y=T)
Chuyển dòng- cột
require(reshape2)
data1=melt(data,id=c("id","sex","group"),measure.vars=c("day1","day2","day3"))
data2=cast(data1,id+sex+group ~ variable)
Phân tích mô tả
Biến phân nhóm
Biến liên tục
require(psych)
describe.by(data, gender, range=F) (phân tích theo gender)
require(gmodels)
CrossTable(gender, anyfx, digits=3, chisq=T, fisher=T)
Kiểm định cơ bản
Kiểm định t test (2 nhóm)
Giả định cho t test
Dữ liệu tuân theo phân phối chuẩn
Hai nhóm độc lập, cỡ mẫu lớn > 30
Thu thập từ lấy mẫu ngẫu nhiên
require(psych)
describe.by(pcs,homeless)
t = t.test(pcs~homeless), (hoặc t.test(x1,x2), print (t)
Beeswarm
require(beeswarm), beeswarm(pcs~homeless, data=ex, col=16, pch=16)
boxplot(pcs~homeless, add=T, data=ex) (vẽ boxplot chồng trên beewarm)
Ghép data
x=c(smokers,nonsmokers)
g1=rep("Smokers",10),g2=rep("Nonsmoker",10),group=c(g1,g2)
Kiểm định phân phối chuẩn
Phương sai của 2 nhóm tương đương nhau
Shapiro-Wilk test
Histogram
hist(data, breaks=20,col="blue", border="white" prob=T), lines(density(na.omit(fnbmd)),col="red",lwd=3)
Qqnorm
plot(density(x)), qqnorm(x), qqline(x,col=2)
shapiro.test(data)
Anderson-Darling test
require(nortest), ad.test(data)
So sánh hai nhóm phương pháp phi tham số (khi không theo phân phối chuẩn)
dat = as.table(matrix(c(101,964,141,921),nrow=2,byrow=T)
epi.2by2(dat,method="cohort.count",conf.level=0.95)
Wilcoxon Rank Sum Test: so sánh rank
wilcox.test(jitter(g1),jitter(g2)) (hoặc wilcox.test(x~group))
Phương pháp hoán vị permutation
sample(x,5,replace=F) lấy 5 phần tử từ x không hoàn lại
require(coin), group=c(rep(1,21),rep(2,23)), oneway_test(x~as.factor(group),distribution=approximate(B=1000))
So sánh 2 nhóm: biến nhị phân z test
require(epiR)
tab=table(gender,anyfx)
epi.2by2(tab,method="cohort.count",conf.level=0..95)
Tương tự ratio risk
Phương pháp Chi squared Wald: prop.test(tab) (cho biết thêm khoảng tin cậy 95%)
Phương pháp hoán vị
chisq_test(tab,data=os,distribution=approximate(B=1000))
scatterplot(fnbmd~age|gender,pch=c(1,16), col=c("red","blue"))
Boxplot
boxplot(wt~sex, main="Weight by gender",xlab="Weight", ylab="Kg", notch=T, col=c("red","blue"), names=c("Female", "Male"), horizontal=T)
require(sciplot), bargraph.CI(anyfx,fnbmd,col="red", names=c("No facture","Facture"))
bargraph.CI(sex,fnbmd,group=anfx,col=c("blue","red"), names=c("Women","Men"))
lineplot.CI(age,fnbmd,anyfx,col=c("blue","red"),legend=F)
require(psych), error.bars.by(fnbmd,gender,bars=T,labels=c("Women","Men"))
Barplot
bmicat=bmi, bmicat[bmi<18]<-"underweight", bmicat[bmi>=18&bmi<=24]<-"Normal", counts=table(bmicat), barplot(counts, horiz=T, cex.axis=1,cex.names=0.5,cex=0.5,las=1,xlab="Math Mark")
Hàm đếm
counts=table(anyfx,gender)
So sánh hai nhóm: odds ratio
require(epiR)
epi.2by2(dat,method="cohort.count",conf.level=0.95)
Phương tích phương sai
Dương tinh giả cho t test khi có nhiều nhóm
Kiểm định CrossTable: chisq, fisher
result=aov(x~ as.factor(group)), summary(result)
Posthoc analysis
TukeyHSD(model), plot(TukeyHSD(model),ordered=T)
pairwise.t.test(gal,group,p.adjust="bonferroni",pool.sd=T)
Phương pháp với khoảng tin cậy gần nhất là tôi ưu nhất
Phân tích phương sai phi tham số
kruskal.test(x~as.factor(group), data=dat)
require(pgirmess), kruskalmc(x, group)
Mô hình tuyến tính cho phương tích phương sai
model=lm(y~group), summary(model), anova(model)
So sánh nhiều nhóm: kiểm định Chi square
Hai biến độc lập thì hệ số tương quan giữa chúng bằng không
So sánh giá trị kì vọng với giá trị thực tế
chisq.test(data)
Có thể dùng anova 2 way
cor.test(age,fnbmd)
require(psych), pairs.panels(vars), corr.test(vars)
Mô hình hồi quy tuyến tính
Giả định
X, Y là tuyến tính về tham số
X không có sai số ngẫu nhiên
Y độc lập với nhau
Sai số ngẫu nhiên có phân bố chuẩn, trung bình 0, phương sai bất biến
var(na.omit(fnbmd)): loại bỏ những giá trị N.A
Plot
plot(varY~varX), abline(lm(varY~varX),col="blue",lwd=3)
require(car), scatterplot(varY~varX|group, pch=16,col="blue",reg.line=F)
p = ggplot(data, aes(x,y)), p + geom_xxx() + theme_xxx()
p = ggplot(data, aes(x=bmi,y=pcfat,fill=gender)), p + geom_point(aes(col=gender)) + geom_smooth() + theme_mw() + xlab("body mass") + ylab("body fat") + ggtitile("distribution")+theme(legend.position="bottom")
summary(lm(lsbmd~age+gender+age:gender,data=example))
Đa cộng tuyến- multicolinearity
Hệ số VIF
Đánh giá tầm quan trọng của biến
require(relaimpo), metric=calc.relimp(model,type=c("lmg"))
boot = boot.relimp(model,b=1000,type=c("lmg"),fixed=F), booteval.relimp(boot,typesel=c("lmg"),level=0.9,bty="perc",nodiff=T)
Chọn mô hình tối ưu
Mô hình sai stepwise
null= lm(sysbp~1,data=ind), result=step(null,sysbp~age+yrmig+wt+ht+chin+fore+calf,direction="both",test="F", k=log(nrow(ind)),trace=0), step$anova, summary(step)
Phương pháp BMA
Cp Mallow
require(peaps),a=regsubset(sysbp!age+yrmig+wt+ht+chin+fore+calf,nbest=3,data=ind), plot(a,scale="Cp")
require(BMA), X=cbind(age, yemig, wt, ht, chin, fore, calf), Y= imd$sysbp, result = bicreg(X,Y,strict=FALSE, OR=20), summary(result), imageplot.bma(result)
Mô hình hồi quy logistic
Odds ratio
epi.2by2(dat,method="cohort.count",conf.level=0.95), #method: case.control, cohort.count, cross.sectional
require(epiR)
Biến phụ thuộc: chỉ có 2 giá trị. Biến tiên lượng: đa dạng
res = glm(cancer~smoking, family=binomial, weight=ntotal), summary(res)
require(rms), dat= datadist(fract), options(datadist="dat"), res=lrm(fx~sex), summary(res)
Chuyển đổi số liệu
require(epicalc), logistic.display(res)
lines(cons,res$fitted, type="1", col="red")
Diễn giải giá trị nhỏ (trên mỗi SD hay đơn vị dễ hiểu): res= glm(fx~I(fnbmd/0/15),family=binomial)
Tiêu chuẩn đánh giá
Deviance càng thấp
AIC càng thấp
LRT có ý nghĩa thống kê của mô hình 2 so với 1
lrt(model1, model2): 2 mô hình phải "nested"
AIC(model1), AIC(model2), AIC(model3)
Chọn mô hình tối ưu
Hàm step
temp = glm(facture~., family="binomial", data=fulldata), search = step(temp), summary(search)
Phương pháp BMA (chỉ số BIC)
require(BMA), xvars = fulldata[,3:7], y = facture
bma.search = bic.glm(xvars,y,strict=F, OR=20, glm.family="binomial"), summary(bma.search), imageplot.bma(bma.search)
Đánh giá độ phân định- Discrimination
require(pROC), p1=predict(model1, type="response"), p2=predict(model2, type="response"), r1= roc(fx, p1), r2= roc(fx,p2), ci(r1), auc(r1), polt.roc(r1), plot.roc(r2,add=T, col="red", lty=2) #vẽ roc2 chồng lên roc1
Đánh giá độ chính xác- Calibration
Đánh giá tái phân nhóm- Reclassification
Ước tính cỡ mẫu
(2/e)^2p(1-p) hoặc (z(a/2)/e)^2p(1-p) với e=(p*(1-p)/n)1/2
e=[khác nhau của khoảng tin cây]/(2*1.96)
(z(a/2)*s/e)^2
Effect size= (p2-p1)/(p1(1-p1)+p2(1-p2))^1/2, N= c(alpha,beta)/ES^2, alpha=falsely positive= 1- Specitivity
require(epiR), epi.studysize(treat=0.85, control=0.9, power=0.9, method="proportions", sigma=NA, n=NA)
require(epiR), epi.studysize(treat=0.84, control=0.8, sigma=0.12, power=0.9, n=NA, method="means")
Effect size= (m2-m1)/s, N= 2*c(alpha,beta)/ES^2
s: độ lệch chuẩn
RMarkdown
Dùng # tiêu đề lớn, ##, ### tiêu đề nhỏ
require(table1), table1(~mpg+factor(cyl)+hp+wt | gear, data=mtcars) #phân nhóm theo gear, #cyl là một biến phân nhóm
require(explore), explore_all(mrcars), explore(mtcars, target=gear), temp=mtcars[,c("mpg", "cyl", "disp", "wt", "gear")], explore_all(temp, target=gear)
{r}: đoạn code trong R
require(GGally), ggpairs(mtcars)
require(gridExtra), grid.arrange(p1, p2, p3 ncol=3): vẽ p1, p2, p3 trên cùng 1 biểu đồ
Knit & Pubish in Rpubs
Phương pháp Bayes
Xác suất có điều kiện
Sai lầm loại 1: xác suất giả thuyết không đúng nhưng bị bác bỏ bởi phép thử hệ thống, DƯƠNG TÍNH GIẢ
Sai lầm loại 2: xác xuất giả thuyết không sai nhưng được chấp nhận bởi phép thử hệ thống, chấp nhận H0 sai, ÂM TÍNH GIẢ
Bằng với mức ý nghĩa được đặt cho kiểm định giả thuyết
Bằng một trừ đi mức ý nghĩa của kiểm định, có thể tăng mức ý nghĩa bằng cách tăng cỡ mẫu
Prior, Likelihood, Posterior
Mô hình beta- binomial
par(mfrow=c(2,3)) - chia màn hình theo dòng [2,3]
theta=seq(from=0, to=1, by=0.01), plot(dbeta(theta,1,1)~theta, type="1", lwd=3, col="red")
Vai trò của xác suất tiền định và các thức xác định
theta=seq(from=0, to=1, by=0.01), alpha =1, beta=1, prior=dbeta(theta, alpha,beta), n=130, x=2, likelihood=theta^x(1-theta)^(n-x), prior.like= priorlikelihood, posterior= prior.like/sum(prior.like), post.prob=dbeta(theta, x+alpha, n-x+beta), mean=weigthed.mean(theta,posterior), var=sum(posterior*(theta-weigted. mean(theta, posterior))^2, sd= sqrt(var), cbind(mean, sd)
require(Bolstad), alpha=1, beta=1, x=2, n=130, res= biobp(x, n, alpha, beta), par(mfrow=c(3,1)), plot(res$prior-res$param.x, type="1", col="red"), plot(res$likelihood-res$param.x, type="1", col="blue"), plot(res$posterior-res$param.x, type="1", col="purple"),
Tìm alpha, beta từ suy luận: theta < 0.5 thường nằm trong khoảng 0.01- 0.2
mean =(0.01+0.2)/2, sd^2=((0.2-0.01)/2)^2, alpha =mean^2((1-mean)/sd^2-1/mean), beta =alpha(1/mean-1), cbind(alpha, beta)
Biên tập với pack tidyverse
dplyr, magrittr, tidyr
temp = bw %>% select(low, bwt, lwt, age)- chỉ lấy 4 biến
temp= filter(bw, race==1, bwt<2500) hoặc temp=bw %>%filter(race==1, bwt<2500)- chọn dòng phù hợp điều kiện
temp= bw%v% filter(race %in% c(1,2)) %>% select(id, age, low, race, bwt) - chọn race bằng 1 hoặc 2
temp= bw %>% mutate (mother.wt=lwt*0.453, weight =bwt/1000)- tạo 2 biến mới mother.wt và weight
arrane(temp, mother.wt, weight)- sắp xếp dữ liệu theo mother.wt và weight
bygroup = group_by(bw, race, smoke), temp= summarize(bygroup, count=n(), mean.age=mean(age, na.rm=T), mean.bw=mean(bwt, na.rm=T)- tạo nhóm, tóm tắt theo số lượng mẫu, mean của age, mean của bw, na.rm remove các giá trị na
temp= bw %>% count(race, smoke)- đếm số lượng
temp= sample_n(bw,10)- lấy mẫu ngẫu nhiên 10 dòng từ bw hoặc temp= sample_frac(bw, 0.05)- lẫy mẫu ngẫu nhiên 5%
temp= bw %>% mutate(group=recode(low, "0" = "No", "1" = "Yes", ethnic=recode(race, "1"= "White", "2"= "Black", "3"= "Hospanics))- tạo biến gropu từ low, ethnic từ race
Phương pháp Bootstrap
Suy luận dựa trên việc tái chọn mẫu có hoàn lại
require(boot), require(simpleboot), temp = one.boot(original_sample, median, R=10000), boot.ci(temp), hist(temp, breaks=20, col="blue", border="white")
n= length(original_sample), B=10000, median= numeric(B), for (i in 1:B) { bs = sample(original_sample,n,replace=T), median[i] =median(bs) }, quantile(median, probs = c(0.025, 0.05, 0.975))
Biến phụ thuộc
Biến độc lập
require(boot), require(simpleboot), temp= two.boot(placebo, coffee, mean, R=10000), boot.ci(temp), hist(temp, breaks=20, col="blue")
require(boot), require(simpleboot), temp= pairs_boot(genetic, foster, FUN=cor, R=1000), boot.ci(temp), plot(temp, breaks=20, col="blue")- function dùng hàm cor
Bootstrap cho tham mô hình hồi quy
x= seq(from=0, to=20, by=0.1), y= 200+ 2x -0.2(x-10)^2 +rnorm(legth(x), mean=10, sd=20), dat= data.frame(x,y), ggplot( data=dat, aes(x=x, y=y)) +geom_point()+ geom_smooth(formula=y~x+x^2) +labs(title="Quadratic model"), m = lm(y~x+I(x^2), data= dat), boot = lm.boot(m, R=1000), s= samples(boot, "coef"), par(mfrow=c(1,3)), hist(s[1,], breaks=20, col="blue"), hist(s[2,], breaks=20, col="blue"), hist(s[3,], breaks=20, col="blue")
Bayes cho 2 tỉ lệ
x1= 652, n1= 10343, x2= 1473, n2= 20173, aplha1= 1, beta1= 1, alpha2= 1, beta2= 1
p1=rbeta(10000, x1+alpha1, n1-x1+beta1), p2= rbeta(10000, x2+alpha2, n2-x2+beta2), diff= p1- p2, plot(density(diff), main="Posterior distribution of diff"), quantile(diff, c(0.025, 0.5, 0.975))- so sánh 2 nhóm
p1=rbeta(10000, x1+alpha1, n1-x1+beta1), p2= rbeta(10000, x2+alpha2, n2-x2+beta2), rr= p1/p2, plot(density(rr), main="Posterior distribution of diff"), quantile(rr, c(0.025, 0.5, 0.975))- phân tích rr
table(rr<0.9)- xác suất thuốc mới giảm nguy cơ tử vong 10% trở lên là bao nhiêu (trên 10000 mẫu)
require(compareGroups), createTable (compareGroups(cut ~ carat + price + color, data=diamonds)
require(compareGroups), createTable (compareGroups(cut ~ carat + price + x + y + z, method=c(price=2), data=diamonds) # xem biến price như biến phi tham số, dùng 1/4Q thay cho Mean
Bayes cho số trung bình- Normal normal model
n likelihood=1, meanp= (m1s0^2+m0s1^2)/(s0^2+s1^2), sp=sqrt((s0^2+s1^2)/(s0^2*s1^2), 1- pnorm(6.4, mean=meanp, sd=sp)
n likelihood>1, sp=sqrt(1/(1/s0^2+1/(s1^2/n))), meanp= sp^2* (m0/s0^2+x#/(s1^2/n)), 1- pnorm(6.4, mean=meanp, sd=sp)
Vẽ biểu đồ bằng R
require(ggplot2), worldmap=map_data("world"), ggplot(data=worldmap, aes(x=long, y=lat, group=group, fill=factor(group)))+geom_polygon(col="white")+ theme(legend.position="none")
countries=c("vietnam", "lao", "thai", "cam"), map= map_data("world", region=countries), gglot(data= map, aes(x=long, y=lat, group=group, fill=factor(group)))+geom_polygon(col="white") + them(legend.position="none")
require(raster, tidyverse, viridis, rgdal), vnm =getData("GADM", country="Vietnam", level=1), vn= fortify(vnm, regions="VARNAME_1"). ggplot(data=vn, aes(x=long, y=lat, group=group))+ geom_polygon(aes(fill=id))+theme(legend.position="none")
require(readxl), pop=read_excel("~provine_data.xlsx"), map=merge(vn, pop, by="id"), ggplot(data=map, aes(x=long, y-lat, group=group))+ geom_polygon(aes(fill=pop.area), show.legend=T)+scale_fill_viridis(direction=-1, option="virisis")
Mô hình hồi quy
Hồi quy tuyến tính
Hồi quy với ML
training data và testing data
k fold validation
Boostrap
Hồi quy tuyền tính
require(capet), control= trainControl(method="cv", number=10), fit= train(form= y~age+ sex+ bmi+ smoking+ children+ region, data=df, method="lm", trControl=control), fit
require(capet), control= trainControl(method="boot", number=1000), fit= train(form= y~age+ sex+ bmi+ smoking+ children+ region, data=df, method="lm", trControl=control), fit
Hồi quy logistic
require(rstanarm, broom, BAS), m= stan_glm(Bodyfat~Abdomen, data=bodyfat), summary(m), prior_summary(m), post.r2= bayes_R2(m)
require(MCMCpack), posterior= MCMCregress(Bodyfat~ Abdomen, data=bodyfat, b0=c(0,0,0), B0=c(0.01, 0.01, 0.01)), plot(posterior), summary(posterior) #b0, B0 lần lượt là mean và sd
require(caret), require(readxl), df= read.excel("insurance data.xlsx"), df$y=log(df$charges), model= lm(y~age+sex+bmi+children+smoking+region, data=df), summary(model)
require(caret), set.seed(123), index= createDataPartition(df$y, p=0.7, list=F), training= df[index, ], testing= df[-index, ], fit= lm(y~age+sex+bmi+children+smoking+region, data=training), summary(fit), testing$predicted= predict(fit, testing), error= testing$predicted- testing$y, sqrt(mean(error^2)), cor(testing$predicted, testing$y)^2
require(caret, readxl), df= read_excel("Diabetes.xlsx"), df$Y= as.factor(df$Outcome), levels(df$Y)= c("No", "Yes"), control= trainControl(method="cv", number=10, summaryFuntion=twoClassSummary, classProbs=T), fit= train(form=Y~Age+ bmi+ Pregnancies+ Glucose, data=df, method="glm", family="binomial, trControl=control), fit
control= trainControl(method="boot", number=10, summaryFuntion=twoClassSummary, classProbs=T), fit= train(form=Y~Age+ bmi+ Pregnancies+ Glucose, data=df, method="glm", family="binomial, trControl=control), fit
index= createDataPartition(df$Y, p=0.7, list=F), training= df[index,], testing=df[-index,], control= trainControl(method="boot", number=1000, summaryFunction=twoClassSummary, classProbs=T), fit= train(form=Y~Age+ bmi+ Pregnancies+ Glucose, data=training, method="glm", family="binomial", trControl=control), testing$predicted= predict(fit, testing), library(caTools), colAUC(testing$predicted, testing$Y, plot=T), confusionMatrix(testing$predicted, testing$Y)
Trong ML hàm glm chỉ dùng được cho biến Yes-No, có thể dùng df$Db= ifelse(df$Outcome==1, "Yes", "No"), ước tính cho Yes trong mô hình glm
testing$predicted=predict(fit, testing, type="prob", testing$predicted=factor(testing$predicted, levels=c("Yes","No") #quy định cho Yes trước
Phân tích Cluster
Tìm ra sự tương đồng dựa vào nhiều biến khác nhau để phân nhóm đối tượng
Chỉ số khoảng cách "distance"= 1- similarity
Chuẩn hóa chỉ số z =(xi-x#)/s
height=c(128, 158, 177, 131), zheight= scale(height)
Xác định chỉ số distance
Euclidean
Manhattan
Minkowski
df=data.frame(nam, height, weight), dist(df[,2:3])
dist(df[,2:3], method="manhattan")
dist(df[,2:3], method= "minkowski")
Tình trạng bất tái lập trong nghiên cứu
Dendrogram
require(datasets, tidyverse, ape, ggdendro,dendextend), data(USArrests), distance= dist(scale(USArrests), method="euclidean"), hc= hclust(distance, method= "ward.D2"), nclust= cutree(hc,4), colors= c("blue", green", "black", "red"), plot(as.phylo(hc), type="fan", tip.color=colors[nclust], label.offset=1, cex=0.7)
USArrests %>% select(Muder, Assault, UrbanPop, Rape) %>% dist() %>% hclust() %>% as.dendrogram() -> dend, plot(dend), dend %>% set(labels_col", value= c("blue", "orange", "purple"), k=3) %>% set("branches_k_color", value= c("blue", "orange", "purple"), k=3) %>% plot(axes=FALSE), rect.dendrogram(dend, k=3, lty=5, lwd=0, x=1, col=rgb(0.1, 0.2, 0.4, 0.1))
k mean clustering
G- statistic
require(factoextra, NbClust), df= scale(USArrests), fviz_nbclust(df, kmeans, nstart= 25, method="gap_stat", nboot=50) + labs(subtitle="Gap statistic method"), km= kmeans(df, centers=3, nstart= 20), fviz_cluster(km, df, ellipse.type="norm")
Eblbow
df= scale(USArrests), fviz_nbclust(df, kmeans, method="wss") + labs(subtitle="Eblbow method")
Average silhouette
df= scale(USArrests), fviz_nbclust(df, kmeans, method="silhouette") + labs(subtitle="Silhouette method")
Xác định k
Chỉ số Hopkin
c1= get_clust_tendency(df, n=50, graph=T), c1$hopkins, print(c1$plot) #n nhỏ hơn số quan sát, Hopkin càng gần 1 càng tốt
gap= clusGap(df, FUN=kmeans, nstart=25, K.max=10, B=100), fviz_gap_stat(gap), km= kmeans(df, 3, nstart=25), fviz_cluster(km, df), sil= silhouette(km$cluster, dist(df)), fviz_silhouette(sil), hcm= eclust(df, "hclust"), fviz_dend(hc, rect=TRUE)