7

mle2パッケージでコマンドを使用しようとしていますbbmle。bbmleBolker による「パッケージによる最尤推定と分析」の p2 を見ています。どういうわけか、正しい開始値を入力できません。再現可能なコードは次のとおりです。

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))

そして、これは私が得るエラーです:

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

理由はありますか？ご協力いただきありがとうございます！

4

1 に答える 1

6

いくつかのことを見逃していますが、最も重要なことは、デフォルトでパラメーターmle2のリストを受け取ることです。代わりにパラメーターベクトルを使用するようにすることもできますが、少し手間がかかります。

所々でコードを微調整しました。(対数尤度関数を負の対数尤度関数に変更しました。これがないと機能しません!)

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)

次に、フィットを行います。主なことは、パラメータ ベクトル内の要素の名前を指定して知らせるvecpar=TRUEために使用することです...parnamesmle2

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))

この特定の例について上で指摘したように、プロビット回帰を再実装したばかりです (ただし、何らかの方法で不均一分散を可能にするためにこれを拡張したいことは理解しています...)

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

最後の注意として、パラメーターのいずれか 1 つとインターフェイスparametersに対して準線形モデルを指定できるようにする引数を確認する必要があると思います。formula

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

PSconfint(foo2)は、このセットアップで正常に動作しているように見えます (要求に応じてプロファイル CI を提供します)。

ae <- function(x,y) all.equal(unname(coef(x)),unname(coef(y)),tol=5e-5)
ae(m1,m2) && ae(m2,m3)
于 2012-05-25T16:43:17.467 に答える