Bien le bonjour,
un petit doc intéressant (mon premier contacte avec R et la régression logistique :D):
-
http://pbil.univ-lyon1.fr/R/pdf/br5.pdfquelques traces de notre passé ... version à la paluche ... la convergence est assez rapide :D
Code : Tout sélectionner
#-----------------------------
# course on glm
# from esteves 2006
#-----------------------------
age <- c(25,32.5,37.5,42.5,47.5,52.5,57.5,65)
n <- c(100,150,120,150,130,80,170,100)
y <- c(10,20,30,50,60,50,130,80)
# proportions
f <- y/n
f
# link-transformation
g <- log(f/(1-f))
# weighting
w <- n*f*(1-f)
# weighted regression
r <- predict(lm(g~age,weights=w))
p <- exp(r)/(1+exp(r))
p
plot(age,f,ylim=c(0,1))
lines(age,p,col="green")
symbols(age,f,circles=w,add=TRUE,inc=0.5,col="darkgreen",lwd=2)
# after iterations
for(i in 1:5){
w <- n*p*(1-p)
gu <- r+(f-p)/p/(1-p)
r <- predict(lm(gu~age,weights=w))
r
p <- exp(r)/(1+exp(r))
p
}
coefficients(lm(gu~age,weights=w))
coefficients(glm(cbind(y,n-y)~age,family=binomial(link="logit")))
#-----------------------------
# binary data
#-----------------------------
age <- sort(rnorm(100,45,10))
y <- rbinom(100,1,1/(1+exp(-(-5+0.1*age))))
n <- length(y)
# proportions
f <- y
# link-transformation
g <- log(f/(1-f))
# weighting
w <- n*f*(1-f)
# weighted regression
r <- predict(lm(g~age,weights=w))
p <- exp(r)/(1+exp(r))
plot(age,f,ylim=c(0,1))
lines(age,p,col="green")
cols1 <- rainbow(5)
for(i in 1:5){
w <- p*(1-p)
gu <- r+(f-p)/p/(1-p)
r <- predict(lm(gu~age,weights=w))
p <- exp(r)/(1+exp(r))
lines(age,p,col=cols1[i])
}
glm1 <- glm(f~age, family=binomial(link="logit"))
lines(age,predict(glm1,type="response"),col="red",lty=2,lwd=2)
coefficients(lm(gu~age,weights=w))
coefficients(glm1)
####
On peut aussi s'amuser à jouer avec d'autres moteurs d'optimisation ...
Code : Tout sélectionner
#----------------------------------------------
# regression logistique (05/05/07 11:21:23)
# avec donn?es binaires sans consid?ration
# des poids
#----------------------------------------------
y <- c(0,0,0,1,0,1,1,1,1,1)
x <- sort(rnorm(10)+runif(10))
plot(x,y)
f <- function (y, x, par)
{
eta <- par[1] + x * par[2]
mu <- 1/(1+exp(-eta))
ll <- sum(y*log(mu)+(1-y)*log(1-mu))
return(-ll)
}
nlm(f, y=y,x=x,p=c(0,0))
glm1 <- glm(y~x,family=binomial)
glm1
logLik(glm1)
# regression logistic avec algorithme génétique
require(gafit)
invlogit <- function(eta) 1/(1+exp(-eta))
e <- expression( - sum(y*log(invlogit(a+b*x))+(1-y)*log(1-invlogit(a+b*x))))
gafit( e, list( a=0, b=0 ), step=0.0001, maxiter=1000)
approche similaire pour régression de poisson ...
Code : Tout sélectionner
#------------------------------------------------
# regression poissonienne (05/05/07 11:41:54)
#----------------------------------------------
y <- sort(rpois(10,1))
x <- sort(rnorm(10)+runif(10))
plot(x,y)
f <- function (y, x, par)
{
eta <- par[1] + x * par[2]
mu <- exp(eta)
# sum(w * (-mu + y * log(mu) - lgamma(y + 1))) avec w = poids
ll <- sum(y*log(mu)-mu-lgamma(y+1))
return(-ll)
}
nlm(f, y=y,x=x,p=c(0,0),steptol = 1e-09,gradtol = 1e-09)
glm1 <- glm(y~x,family=poisson)
glm1
logLik(glm1)
optim(par=c(0,0), f,x=x,y=y)
nlminb(c(0,0), f,x=x,y=y)
# regression de poisson avec algorithme génétique
require(gafit)
e <- expression( - sum(y*log(exp(a+b*x))-exp(a+b*x)-lgamma(y+1)))
gafit( e, list( a=0, b=0 ), step=0.0001, maxiter=1000)
je dois aussi avoir des versions en mode BUGS/JAGS quelques parts qui traînent sur un des mes disques :D
Il faut faire attention aux paramètres de contrôle, tolérance, d’itération,epsilon, ... (voir par exemple glm.control()).
Ces différentes fonctions utilisent rarement les mêmes options par défaut. On avait même eu quelques surprises entre l'implémentation
des GLM de S et de R pendant le projet EFI+ (clin d'oeil pour max :P).
pour le détail du fonctionnement de la function "glm", le manuel de GLIM (qui est quelques parts sur mon bureau ...mais où??? :'( ) est très informatif. Sinon voir aussi Le McCullagh and Nelder (Generalized linear model, 1989, cf post de Eric) ainsi que l'article fondateur (?) de Nelder and Wedderburn(1972). "Generalized Linear Models" (
http://www.jstor.org/stable/2344614?ori ... b_contents).
en espérant ne pas être trop hors sujet et en espérant avoir aidé un peu :)
@++
pierre