Tuesday, September 23, 2014

Các hàm R phổ biến 4: Mô hình hồi qui logistic, binomial, Poisson

Dưới đây là một số hàm có thể dùng cho mô hình hồi qui logistic. Hàm glm trong R được dùng chủ yếu. Hàm glm không cung cấp OR và khoảng tin cậy 95%, nên chúng ta cần dùng hàm logistic.display trong package epicalc.  Ngoài glm, còn có hàm lrm trong package chuyên dụng rms (Frank Harrell).



Phân tích 
Hàm R
Mô hình hồi qui logistic: biến y là nhị phân, biến x có thể là nhị phân hay liên tục

m = glm(y ~ x, family="binomial");
summary(m)
library(epicalc)
logistic.display(m)


Tính toán giá trị tiên lượng;vẽ biểu đồ 

predict(m, type="response");
plot(x, fitted(glm(y ~ x, family="binomial")))

Mô hình hồi qui logistic đa biến với 3 biến x1, x2, x3

m = glm(y ~ x1 + x2 + x3, family="binomial");
summary(m)

Chọn mô hình “tối ưu” bằng phương pháp Bayesian Model Average
library(BMA);
xvars = cbind(x1,x2,x3,x4,x5);
bma = bic.glm(xvars, y, strict=F, OR=20,glm.family="binomial");
imageplot.bma(bma)

Mô hình hồi qui Poisson
age=c(19.5,29.5,39.5,49.5,59.5,69.5,79.5);
cases = c(1, 16, 30, 71, 102, 130, 133);
pop = c(172675, 123065, 96216, 92051, 72159, 54722, 32185);
dat = data.frame(age, cases, pop)
m = glm(cases ~ age + offset(log(pop)), family=poisson,data=dat)
summary(fit)

Mô hình hồi qui Binomial
m = glm.nb(daysabs ~ math + prog, data = dat))
p = predict(m, newdata, type="response")




Mô hình hồi qui Poisson

p = read.csv("http://www.ats.ucla.edu/stat/data/poisson_sim.csv")

p = within(p, {
   prog <- factor(prog, levels=1:3, labels=c("General", "Academic", "Vocational"))
id <- factor(id)
})

m1 = glm(num_awards ~ prog + math, family="poisson", data=p))
 
p$phat = predict(m1, type="response")
p = p[with(p, order(prog, math)), ]
 
ggplot(p, aes(x = math, y = phat, colour = prog)) +
geom_point(aes(y = num_awards), alpha=.5, position=position_jitter(h=.2)) +
geom_line(size = 1) +
labs(x = "Math Score", y = "Expected number of awards")
 


No comments:

Post a Comment