Regression logistique

[Anne Quillevere]

Bonjour,

J'étudie la vitesse de recouvrement d'un couvert végétal.
Pour cela, j'utilise la fonction logistique suivante:

logit<-function(a1,b1,c1,t)(a1/(1+((a1/3)-1)*exp((b1-t)/c1)))
avec a1 qui doit être compris entre 0 et 100 (pourcentage de recouvrement max)

Lorsque j'utilise la fonction nlsList pour réaliser un fit de la fonction selon les variétés, pour chaque année, les valeurs de a1 sont parfois supérieures à 100.

Je souhaiterais donc savoir comment borner les valeurs de a1 [0,100] dans mon équation.

Merci d'avance pour votre aide

[Eric Casellas]

Salut,

Je n'utilise pas cette fonction donc je vais peut-être dire des bêtises, mais de ce que j'ai vu il semblerait que ce ne soit pas possible directement avec cette fonction nlsList (c'est par contre possible avec des fonctions comme optim ou optimize).
Mais si tu veux garder nlsList je suis tombé sur ce post (non testé) qui pourrait être une piste : https://stackoverflow.com/questions/27570758/r-how-to-use-bounds-and-the-port-algorithm-in-nlslist

Eric

[François Bonnot]

Bonjour,
Je n'utilise pas non plus cette fonction, lui préférant nlm ou optim.
Mais quelle que soit la fonction de régression, une solution consiste à remplacer le paramètre sous contrainte (ici a1) par f(a2) tel que 0 < a1 < 100 si -Inf < a2 <+Inf.
Par exemple :

Code : Tout sélectionner

f <- function(x) (atan(x)+pi/2)*100/pi
curve(f,-10,10)

Votre fonction devient alors :

Code : Tout sélectionner

logit<-function(a2,b1,c1,t)(f(a2)/(1+((f(a2)/3)-1)*exp((b1-t)/c1)))

et le paramètre d'intérêt final est f(a2).

[Logez Maxime]

Bonjour,

J'ai essayé la solution de François sur des données simulées et j'arrive à un problème de "gradient singulier" en utilisant la fonction f ou une autre pour ramener a1 entre 0 et 100. Peut-être est-ce du à mon exemple, mais peut-être que le cas va se représenter.

Une autre solution plus artisanale peut-être tout simplement de récupérer les modèles pour lesquels les valeurs de a1 sortent de la gamme attendue et de ré-estimer les paramètres en fixant a1 à 100 si son estimation est supérieure à 100 et 0 si <0.

Code : Tout sélectionner

nls(y ~ logit(100, b1, c1, t), tab, start = list(b1= 10, c1 = 10))

Comme ceci tu fixes la valeur de a1 et tu n'estimes plus que ce qui se passe pour b1 et c1. Attention toutefois aux répercussions sur les degrés de libertés, matrice de variance-covariance, estimation des écarts-types des paramètres, etc.
Tout dépend de ce que tu cherches à faire.

La deuxième solution s'est de réfléchir au pourquoi du comment, de se demander si le modèle est adapté et à la pertinence des résultats quand tu bornes a1 et à la signification de tout ça, mais là on sort largement du cadre de ce forum.

Cordialement,
Maxime

[Eric Wajnberg]

Je n'utilise pas non plus cette fonction, lui préférant nlm ou optim.

Effectivement, je préfère également optim() qui a des arguments lower et upper qui permettent de faire ceci.

HTH, Eric.

[Anne Quillevere]

Merci pour vos retours.

J'ai essayé la solution de François sur des données simulées et j'arrive à un problème de "gradient singulier" en utilisant la fonction f ou une autre pour ramener a1 entre 0 et 100. Peut-être est-ce du à mon exemple, mais peut-être que le cas va se représenter.

Je tombe également sur le même problème en utilisant la fonction f.

J'ai opté pour la méthode artisanale de Maxime qui fonctionne très bien :)

[François Bonnot]

Bonjour,

Sur un jeu de données simulé dont la régression donne une valeur de a1 > 100, les méthodes (1) optim avec reparamétrage de a1, (2) optim avec a1 fixé à 100, et (3) optim avec upper=100 donnent des résultats quasi identiques (je peux poster le code si ça intéresse quelqu'un).

Les essais ci-dessus ont été réalisés (pour aller vite) par minimisation des sommes de carrés résiduels (comme le fait d'ailleurs la fonction nlsList), mais ce n'est pas la meilleure méthode pour ajuster des proportions. Si les données ne présentent pas un plateau évident inférieur à 100, il est plus pertinent d'effectuer une vraie régression logistique avec par exemple un modèle glm(y/100~t,family=binomial) qui comporte alors 2 paramètres et une asymptote à 100. Le test réalisé sur le jeu d'essai donne des résultats voisins des 3 méthodes ci-dessus.

[Anne Quillevere]

Bonjour,

Sur un jeu de données simulé dont la régression donne une valeur de a1 > 100, les méthodes (1) optim avec reparamétrage de a1, (2) optim avec a1 fixé à 100, et (3) optim avec upper=100 donnent des résultats quasi identiques (je peux poster le code si ça intéresse quelqu'un)

Effectivement, cela serait intéressant d'avoir le code.

Merci d'avance
Anne

[François Bonnot]

Voici le code R.
Les données ont été simulées avec la méthode figurant en fin de script et ont été volontairement tronquées à t=50 pour auroriser des valeurs estimées de a1 tantôt inférieures à 100 et tantôt supérieures.

Code : Tout sélectionner

## données simulées
df <- structure(list(t = 0:50, y = c(3L, 3L, 2L, 2L, 4L, 2L, 2L, 3L,
8L, 3L, 5L, 7L, 6L, 8L, 11L, 4L, 8L, 5L, 10L, 8L, 12L, 16L, 11L,
21L, 16L, 23L, 23L, 20L, 20L, 19L, 22L, 31L, 32L, 27L, 38L, 30L,
38L, 39L, 42L, 42L, 47L, 48L, 56L, 47L, 67L, 63L, 67L, 65L, 54L,
72L, 71L)), .Names = c("t", "y"), row.names = c(NA, -51L), class = "data.frame")

plot(df,xlim=c(0,100),ylim=c(0,110))

## fonction de reparamétrage
f <- function(x) (atan(x)+pi/2)*100/pi

## fonctions correspondant ux 3 méthodes
logit1 <- function(a2,b1,c1,t)(f(a2)/(1+((f(a2)/3)-1)*exp((b1-t)/c1))) ## fonction reparamétrée
logit2 <- function(b1,c1,t)(100/(1+((100/3)-1)*exp((b1-t)/c1))) ## fonction bornée à 100
logit3 <- function(a1,b1,c1,t)(a1/(1+((a1/3)-1)*exp((b1-t)/c1))) ## fonction d'origine

## sommes de carrés à minimiser pur les 3 fonctions logit
s1 <- function(p,t,y) sum((y-logit1(p[1],p[2],p[3],t))^2)
s2 <- function(p,t,y) sum((y-logit2(p[1],p[2],t))^2)
s3 <- function(p,t,y) sum((y-logit3(p[1],p[2],p[3],t))^2)

## paramètres initaux (différents de ceux qui ont servi à simuler)
a2 <- 0 ; a1 <- f(a2) ; b1 <- 0 ; c1 <- 5
init1 <- c(a2,b1,c1)
init2 <- c(b1,c1)
init3 <- c(a1,b1,c1)

## t100 sera utilisé pour tracer les courbes complétées
t100 <- 0:100

## régression sans contrainte (avec fonction d'origine logit3)
z0 <- optim(init3,s3,t=df$t,y=df$y)
lines(t100,logit3(z0$par[1],z0$par[2],z0$par[3],t100),col="blue",lwd=2)
abline(h=z0$par[1],col="blue",lwd=2)

## méthode 1
z1 <- optim(init1,s1,t=df$t,y=df$y)
lines(t100,logit1(z1$par[1],z1$par[2],z1$par[3],t100),col="green",lwd=2)
abline(h=f(z1$par[1]),col="green",lwd=2)

## méthode 2
z2 <- optim(init2,s2,t=df$t,y=df$y)
lines(t100,logit2(z2$par[1],z2$par[2],t100),col="red",lwd=2,lty=2)
abline(h=100,col="red",lwd=2,lty=2)

## méthode 3
z3 <- optim(init3,s3,upper=100,t=df$t,y=df$y,method="L-BFGS-B") ## warning sans importance
lines(t100,logit3(z3$par[1],z3$par[2],z3$par[3],t100),col="black",lwd=3,lty=3)
abline(h=z3$par[1],col="black",lwd=3,lty=3)

## régression logistique avec glm
model <- glm(y/100~t,data=df,family=binomial)
yy1 <- model$coef[1]+model$coef[2]*t100
yy <- exp(yy1)/(1+exp(yy1))*100
## (à la fin le graphe n'est plus très lisible)
points(t100,yy,pch=2)

## méhode de simulation des données
a1 <- f(8) ; b1 <- 5 ; c1 <- 10
t <- 0:50
prob <- logit3(a1,b1,c1,t)/100
df <- data.frame(t=t,y=rbinom(length(t),100,prob))