Profil Konfidenzintervalle in R: mle2

Question

Ich versuche, den Befehl mle2 zu verwenden, in dem Paket bbmle. Ich schaue auf p2 von "Maximum likelihood Schätzung und Analyse mit dem bbmle Paket" von Bolker. Irgendwie kann ich die richtigen Startwerte nicht eingeben. Hier ist der reproduzierbaren Code:

l.lik.probit <-function(par, ivs, dv){ 
Y <- as.matrix(dv) 
X <- as.matrix(ivs) 
K <-ncol(X) 
b <- as.matrix(par[1:K]) 
phi <- pnorm(X %*% b) 
sum(Y * log(phi) + (1 - Y) * log(1 - phi)) 
} 

n=200 

set.seed(1000) 

x1 <- rnorm(n) 
x2 <- rnorm(n) 
x3 <- rnorm(n) 
x4 <- rnorm(n) 

latentz<- 1 + 2.0 * x1 + 3.0 * x2 + 5.0 * x3 + 8.0 * x4 + rnorm(n,0,5) 

y <- latentz 
y[latentz < 1] <- 0 
y[latentz >=1] <- 1 
x <- cbind(1,x1,x2,x3,x4) 
values.start <-c(1,1,1,1,1) 

foo2<-mle2(l.lik.probit, start=list(dv=0,ivs=values.start),method="BFGS",optimizer="optim", data=list(Y=y,X=x)) 

und das ist der Fehler, den ich bekommen:

Error in mle2(l.lik.probit, start = list(Y = 0, X = values.start), method = "BFGS", : 
    some named arguments in 'start' are not arguments to the specified log-likelihood function 

Jede Idee, warum? Danke für Ihre Hilfe!

Answer 1

Sie haben ein paar Dinge verpasst, aber das wichtigste ist, dass mle2 standardmäßig eine Liste der Parameter dauert; Sie können stattdessen einen Parameter Vektor nehmen, aber Sie müssen ein bisschen härter arbeiten.

Ich habe den Code an einigen Stellen leicht verändert. (Änderte ich die Log-Likelihood-Funktion zu einem negativen Log-Likelihood-Funktion, ohne das dies würde nie funktionieren!)

l.lik.probit <-function(par, ivs, dv){ 
    K <- ncol(ivs) 
    b <- as.matrix(par[1:K]) 
    phi <- pnorm(ivs %*% b) 
    -sum(dv * log(phi) + (1 - dv) * log(1 - phi)) 
} 

n <- 200 

set.seed(1000) 

dat <- data.frame(x1=rnorm(n), 
        x2=rnorm(n), 
        x3=rnorm(n), 
        x4=rnorm(n)) 

beta <- c(1,2,3,5,8) 
mm <- model.matrix(~x1+x2+x3+x4,data=dat) 
latentz<- rnorm(n,mean=mm%*%beta,sd=5) 

y <- latentz 
y[latentz < 1] <- 0 
y[latentz >=1] <- 1 
x <- mm 
values.start <- rep(1,5) 

Jetzt machen wir die Passform. Die Hauptsache ist vecpar=TRUE zu spezifizieren und parnames zu verwenden mle2 kennen die Namen der Elemente im Parametervektor zu lassen ...

library("bbmle") 
names(values.start) <- parnames(l.lik.probit) <- paste0("b",0:4) 
m1 <- mle2(l.lik.probit, start=values.start, 
      vecpar=TRUE, 
      method="BFGS",optimizer="optim", 
      data=list(dv=y,ivs=x)) 

Wie oben ausgeführt für dieses Beispiel haben Sie gerade die Probit-re umgesetzt Regression (obwohl ich verstehe, dass Sie jetzt diese erweitern wollen für heteroscedasticity in irgendeiner Art und Weise zu ermöglichen, ...)

dat2 <- data.frame(dat,y) 
m2 <- glm(y~x1+x2+x3+x4,family=binomial(link="probit"), 
    data=dat2) 

Als abschließende Bemerkung, würde ich sagen, dass Sie das parameters Argument sollten überprüfen, die Sie erlaubt ein sublineares Modell für einen der Parameter anzugeben, und die formula Schnittstelle:

m3 <- mle2(y~dbinom(prob=pnorm(eta),size=1), 
      parameters=list(eta~x1+x2+x3+x4), 
      start=list(eta=0), 
      data=dat2) 

PS confint(foo2) scheint gut zu funktionieren (was Profil CIs wie gewünscht) mit diesem Set-up.

ae <- function(x,y) all.equal(unname(coef(x)),unname(coef(y)),tol=5e-5) 
ae(m1,m2) && ae(m2,m3)