Asymptotes de régression binomiale

La régression logistique binomiale a des asymptotes supérieures et inférieures de 1 et 0 respectivement. Cependant, les données de précision (à titre d'exemple) peuvent avoir des asymptotes supérieures et inférieures très différentes de 1 et / ou 0. Je peux voir trois solutions potentielles à cela:

1. Ne vous inquiétez pas si vous obtenez de bons ajustements dans la zone d'intérêt. Si vous n'obtenez pas de bons ajustements, alors:
2. Transformez les données de sorte que le nombre minimum et maximum de réponses correctes dans l'échantillon donne des proportions de 0 et 1 (au lieu de dire 0 et 0,15).
   ou
3. Utilisez une régression non linéaire pour pouvoir spécifier les asymptotes ou demander à l'installateur de le faire pour vous.

Il me semble que les options 1 et 2 seraient préférées à l'option 3 en grande partie pour des raisons de simplicité, auquel cas l'option 3 est peut-être la meilleure option car elle peut fournir plus d'informations?

modifier
Voici un exemple. La correction totale possible pour la précision est de 100, mais la précision maximale dans ce cas est de ~ 15.

accuracy <- c(0,0,0,0,0,1,3,5,9,13,14,15,14,15,16,15,14,14,15)
x<-1:length(accuracy)
glmx<-glm(cbind(accuracy, 100-accuracy) ~ x, family=binomial)
ndf<- data.frame(x=x)
ndf$fit<-predict(glmx, newdata=ndf, type="response")
plot(accuracy/100 ~ x)
with(ndf, lines(fit ~ x))

L'option 2 (selon les commentaires et pour clarifier ma signification) serait alors le modèle

glmx2<-glm(cbind(accuracy, 16-accuracy) ~ x, family=binomial)

L'option 3 (par souci d'exhaustivité) s'apparenterait à:

fitnls<-nls(accuracy ~ upAsym + (y0 - upAsym)/(1 + (x/midPoint)^slope), 
  start = list("upAsym" = max(accuracy), "y0" = 0, "midPoint" = 10, "slope" = 5), 
  lower = list("upAsym" = 0, "y0" = 0, "midPoint" = 1, "slope" = 0), 
  upper = list("upAsym" = 100, "y0" = 0, "midPoint" = 19, hillslope = Inf), 
  control = nls.control(warnOnly = TRUE, maxiter=1000),
  algorithm = "port")

Question interessante. Une possibilité qui me vient à l'esprit est d'inclure un paramètre supplémentaire afin de contrôler la limite supérieure de la fonction 'link'.$p\in[0,1]$

Soit , des observations indépendantes, où , , est un vecteur de variables explicatives, est un vecteur de coefficients de régression et est la fonction de lien. Ensuite, la fonction de vraisemblance est donnée par$\{{\bf x}_j,y_j,n_j\}$$j=1,...,n$$y_j\sim \text{Binomial}\{n_i,pF({\bf x}_j^T\beta)\}$$p\in[0,1]$${\bf x}_j=(1,x_{j1},...,x_{jk})^T$$\beta=(\beta_0,...,\beta_k)$$F^{-1}$

L'étape suivante consiste à choisir un lien, par exemple la distribution logistique et à trouver le MLE correspondant de .$(\beta,p)$

Considérons l'exemple de jouet simulé suivant utilisant un modèle dose-réponse avec et$(\beta_0,\beta_1,p)=(0.5,0.5,0.25)$$n=31$

dose = seq(-15,15,1)
a = 0.5
b = 0.5
n=length(dose)
sim = rep(0,n)
for(i in 1:n) sim[i] = rbinom(1,100,0.25*plogis(a+b*dose[i]))

plot(dose,sim/100)

lp = function(par){
if(par[3]>0& par[3]<1) return(-(n*mean(sim)*log(par[3]) +  sum(sim*log(plogis(par[1]+par[2]*dose)))  + sum((100-sim)*log(1-par[3]*plogis(par[1]+par[2]*dose))) ))
else return(-Inf)
}

optim(c(0.5,0.5,0.25),lp)

L'un des résultats que j'ai obtenus est . Par conséquent, cela semble exact. Bien sûr, une exploration plus détaillée de ce modèle serait nécessaire car l'inclusion de paramètres dans un modèle de régression binaire peut être délicate et des problèmes d'identification ou d'existence du MLE peuvent sauter sur l'étape 1 2 .$(\hat\beta_0,\hat\beta_1,\hat p)=( 0.4526650, 0.4589112, 0.2395564)$

Éditer

Compte tenu de la modification (qui modifie considérablement le problème), la méthode que j'ai proposée précédemment peut être modifiée pour ajuster les données que vous avez fournies. Considérez le modèle

où est le CDF logistique, est un paramètre d'emplacement, est un paramètre d'échelle et le paramètre contrôle la hauteur de la courbe de la même manière que dans l'ancien modèle. Ce modèle peut être ajusté en utilisant des moindres carrés non linéaires . Le code R suivant montre comment procéder pour vos données.$F$$\mu$$\sigma$$p$

rm(list=ls())
y = c(0,0,0,0,0,1,3,5,9,13,14,15,14,15,16,15,14,14,15)/100
x = 1:length(y)
N = length(y)

plot(y ~ x)

Data = data.frame(x,y)

nls_fit = nls(y ~ p*plogis(x,m,s), Data, start = list(m = 10, s = 1,  p = 0.2) )

lines(Data$x, predict(nls_fit), col = "red")

J'utiliserais le maximum du vecteur X comme le nombre total possible de succès. (Il s'agit d'une estimation biaisée du véritable nombre maximal de succès, mais cela devrait fonctionner assez bien si vous avez suffisamment de données).

accuracy <- c(0,0,0,0,0,1,3,5,9,13,14,15,14,15,16,15,14,14,15)
x<-1:length(accuracy)
glmx<-glm(cbind(accuracy, max(accuracy)-accuracy) ~ x, family=binomial)
ndf<- data.frame(x=x)
ndf$fit<-predict(glmx, newdata=ndf, type="response")
plot(accuracy/max(accuracy) ~ x)
with(ndf, lines(fit ~ x))

Cela crée un tracé qui ressemble à:

Notez que la régression binomiale est basée sur une réponse binaire pour chaque cas individuel. chaque réponse individuelle doit pouvoir prendre l'une des deux valeurs. S'il y a une certaine limite à la proportion, il doit également y avoir eu des cas qui ne pouvaient prendre qu'une seule valeur.

Il semble que vous n'ayez pas affaire à des données binaires mais à des données sur une plage finie. si tel est le cas, la régression bêta semble plus appropriée. Nous pouvons écrire la distribution bêta comme suit:

Vous définissez ensuite comme n'importe quelle fonction de lien qui mappe l'intervalle dans les réels. Il existe un package R qui peut être utilisé pour s'adapter à ces modèles, bien que je pense que vous devez connaître les limites. Si vous le faites, redéfinissez la nouvelle variable .$g(\mu_i)=x_i^T\beta$$[L,U]$$y_i=\frac{d_i-L}{U-L}$