• 统计学数理统计

  • WinBUGS问题：Panel data模型中的Incidental Parameter Problem.有人知道如何解决吗？

[quote]引用第17楼hangover于2007-04-17 02:05发表的“”:





如果用Likelihood可以做的，何必再用MCMC？missing data 可否用prior补一下？posterior is propotional to the product of likelihood and prior。 likelihood principal 还是有效的。

如果spiegelhalter说不行，winbugs肯定不行，至少现在。那你只能自己写code了。你把模型放上来，我可以帮你想想。winbugs我用过，很久以前了，不过现在实在没有时间去查manual 去翻译模型。隔行如隔山， 我觉得你的模型很复杂， 但是也很想学习一下。[/quote]



谢谢楼上的指点.我打算用data imputation的方法来弥补缺失数据.所以这个目前不是大问题.至于你说的模型，我太不理解.我上面的code里就是一个标准的panel data probit model了.里面复杂的部分是关于如何 orthogonal reparametrizations.按照Tony Lancaster的说法，是为了“This is an attempt to reduce the dependence of inference about beta on the choice of prioior for the individual effects (alpha i).”说老实话他关于如何处理incidental parameter problenm的部分我没有看得太懂.我也不会更复杂的编程了.我打算模仿Tony的代码改写成包含多个自变量的codes.

就是这个部分

alpha ~ dnorm(nu, 1)

nu <- lam [i, z]

z ~ dcat (p[ ])

lam[i, 1]   <-   -beta * x[i, 1]

lam[i, 2]   <-   -beta * x[i, 2]

lam[i, 3]   <-   -beta * x[i, 3]

lam[i, 4]   <-   -beta * x[i, 4]

lam[i, 5]   <-   -beta * x[i, 5]

lam[i, 6]   <-   -beta * x[i, 6]

lam[i, 7]   <-   -beta * x[i, 7]

lam[i, 8]   <-   -beta * x[i, 8]

lam[i, 9]   <-   -beta * x[i, 9]

lam[i, 10] <-   -beta * x[i, 10]

lam[i, 11] <-   -beta * x[i, 11]

lam[i, 12] <-   -beta * x[i, 12]

lam[i, 13] <-   -beta * x[i, 13]
[quote]引用第20楼revell于2007-04-17 18:48发表的“”:

“This is an attempt to reduce the dependence of inference about beta on the choice of prioior for the individual effects (alpha i).”[/quote]



alpha is the fixed-effect? I guess your method is similar to the so-called conditional estimation. it is still not clear to me, i must read the book or the model. can not help you, sorry.
Fine. I will manage to upload related pages in Tony Lancaster's book.
看了一下，挺有帮助。把书传上来吧，大家可以讨论一下。
[quote]引用第0楼revell于2007-04-10 23:13发表的“WinBUGS问题：

但是问题是，Tony的书里面只给了一个自变量x[i,j]的代码.数据是二维的（N*T）.他在这个代码后面说道，如果自变量是一个向量，包含x=(x1,x2,x3,.......)的话，就是我的情况。必须对代码进行改写，数据也要变成3维的（N*T*K,K就是自变量的数量，比如=4）.我不知道该如何改写代码和数据，请大家帮我看看吧.

.......[/quote]



先谢谢你上传的论文.



我觉得x是什么东西并不是很重要,关键是要把likelihood写出来, pp.306 最后一个公式,决定因数是"H", the probability distribution function. 这个不应该是什么难题, 硬活, 可能增加几个loop就可以了. 有了likelihood之后再考虑如何specify prior. 他这里的问题是 objective information(likelihood) can not dominate subjective belief(prior)，就是likelihood不够强，数据量不够。

看他的两个例子pp.304,一个是informative prior (the precision of fixed/effect has a gamma prior) good result, 另外一个是noninformative prior(flat uniform) which produces inconsistent parameter estimates. 他是reparametrize (redefine) the fixed-effect in order to have orthogonalized likelihoods. 然后选择一个prior让这个新的fixed-effect is prior independent of beta. 也是只有在这种情况下, bayesian inference about beta can be independent of knowledge about the new fixed-effect. 觉得论文比书说得明白些, 可能是前几章没看的原因吧. 还是奇怪, beta和原始的fixed-effect应该有关联阿. 他的解释是"If we can find

such a reparametrization there are incidental parameters but there is, in general, no incidental

parameter problem." 不解??????



data augmentation对于mixture distribution是非常通用的办法,在前面的帖子里已经讨论过了.
谢谢hangover的解释。我对likelihood之类的概念掌握的不是很好。没办法写出likelihood,再根据likelihood写出WinBUGS程序.所以基本上是改写别人的代码，当然在概念上尽量搞懂为什么这样做.：（



我最近改写了Tony的代码，试着把多个自变量放进去，而且还是二维的数据.我用WinBUGS试过了可以运行.但是不知道代码对不对.如果你看得懂pp.307页的代码的话，能不能帮我看看我改写的代码是否符合orthogonal reparametrization的思想.下面的代码包含两个自变量（x和num），以及数据和inits.简单的解释一下，pz是指定每一时期的lam都是1/13的概率，而pw是指定每个自变量对lam赋值的时候都有1/2的概率.也许是比较stupid的方法？



我另外还将二维数据改写成了三维数据，但是很怀疑是否正确。所以如果下面的应用于二维数据的代码可以用的话，我情愿不去碰三维数据.



MODEL2.1         treating incidental parameter problem, two independent variables

model{

for( i in 1 : N ) {

for( j in 1 : T ) {

y[i, j] ~ dbern(mu[i, j])

mu[i, j] <- phi(alpha +beta * x[i, j] + beta.num * num[i, j] + beta.yr[j])}

alpha ~ dnorm(nu, 1)

nu <- lam[i, z, w]

z ~ dcat (pz[])

w ~ dcat (pw[])

lam[i, 1, 1] <- -beta * x[i, 1]

lam[i, 2, 1] <- -beta * x[i, 2]

lam[i, 3, 1] <- -beta * x[i, 3]

lam[i, 4, 1] <- -beta * x[i, 4]

lam[i, 5, 1] <- -beta * x[i, 5]

lam[i, 6, 1] <- -beta * x[i, 6]

lam[i, 7, 1] <- -beta * x[i, 7]

lam[i, 8, 1] <- -beta * x[i, 8]

lam[i, 9, 1] <- -beta * x[i, 9]

lam[i, 10, 1] <- -beta * x[i, 10]

lam[i, 11, 1] <- -beta * x[i, 11]

lam[i, 12, 1] <- -beta * x[i, 12]

lam[i, 13, 1] <- -beta * x[i, 13]

lam[i, 1, 2] <- -beta.num * num[i, 1]

lam[i, 2, 2] <- -beta.num * num[i, 2]

lam[i, 3, 2] <- -beta.num * num[i, 3]

lam[i, 4, 2] <- -beta.num * num[i, 4]

lam[i, 5, 2] <- -beta.num * num[i, 5]

lam[i, 6, 2] <- -beta.num * num[i, 6]

lam[i, 7, 2] <- -beta.num * num[i, 7]

lam[i, 8, 2] <- -beta.num * num[i, 8]

lam[i, 9, 2] <- -beta.num * num[i, 9]

lam[i, 10, 2] <- -beta.num * num[i, 10]

lam[i, 11, 2] <- -beta.num * num[i, 11]

lam[i, 12, 2] <- -beta.num * num[i, 12]

lam[i, 13, 2] <- -beta.num * num[i, 13]

}

for(j in 1:T){

beta.yr[j] ~ dnorm(0,0.001)

}

beta ~ dnorm(0, 0.001)

beta.num ~ dnorm(0, 0.001)

}



init for model 1 and 2 only

list(alpha=c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0),

    beta = 0,

    beta.num = 0,

    beta.yr = c(NA,0,0,0,0,0,0,0,0,0,0,0,0),

)



Data 2 for Model 2 and model 2.1 only



list(

  N=10,T=13,

  pz=c(0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923, 0.076923),

  pw=c(0.5, 0.5),

  y=structure(

      .Data=c(

0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,

1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,

0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 1,

0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0,

0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0,

0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,

0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,

0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0),

          .Dim=c(10,13)

),

      x=structure(

        .Data=c(

0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0,

0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,

1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0,

1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1,

1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,

0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0,

1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0,

0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,

0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 0,

1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0),

          .Dim=c(10,13)

),    

    num=structure(

            .Data=c(

0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 1.0986, 1.0986, 1.0986, 1.0986,

1.6094, 1.6094, 1.6094, 1.6094, 1.6094, 1.6094, 1.6094, 1.7918, 1.7918, 1.7918, 1.7918, 1.7918, 1.7918,

2.3979, 2.4849, 2.5649, 2.5649, 2.5649, 2.5649, 2.5649, 2.6391, 2.6391, 2.6391, 2.6391, 2.6391, 2.6391,

2.6391, 2.6391, 2.6391, 2.6391, 2.7081, 2.7081, 2.7081, 2.7081, 2.7726, 2.7726, 2.9444, 2.9444, 2.9957,

2.7081, 2.7081, 2.7081, 2.7726, 2.8904, 2.9444, 3.0445, 3.0445, 3.1355, 3.1355, 3.1781, 3.1781, 3.1781,

2.1972, 2.1972, 2.3026, 2.3026, 2.3026, 2.3979, 2.3979, 2.3979, 2.3979, 2.3979, 2.3979, 2.3979, 2.3979,

1.7918, 1.7918, 1.7918, 1.7918, 1.7918, 1.7918, 1.7918, 2.0794, 2.1972, 2.3026, 2.3026, 2.3026, 2.3026,

0.0000, 0.0000, 0.0000, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931, 0.6931,

1.9459, 1.9459, 2.0794, 2.0794, 2.0794, 2.0794, 2.0794, 2.0794, 2.0794, 2.0794, 2.1972, 2.1972, 2.1972,

2.3026, 2.3026, 2.3026, 2.4849, 2.4849, 2.5649, 2.8332, 2.8332, 2.8332, 2.9444, 3.0445, 3.0445, 3.0445),  

          .Dim=c(10, 13)            

    )

)

END
[quote]引用第28楼revell于2007-04-19 18:12发表的“”:

.如果你看得懂pp.307页的代码的话，能不能帮我看看我改写的代码是否符合orthogonal reparametrization的思想.

.......[/quote]



不好意思的说，我根本没看他的代码。建议楼主用几个小时，一天，或者几天时间把他的论文推倒一下，相信剩下的东西会迎刃而解的。
Fine.Thanks for the help and I will try.
1 年 后
请问楼主：你的面板数据导入问题解决了吗？我现在用winbugs解一个面板数据的随即前沿生产函数模型，也遇到和楼主一样的问题，就是有多个解释变量，不知道该怎么把这些数据导入程序，老外的原文这里也是语焉不详，该怎么办啊？
这里给出的是随机前沿的成本函数，模型如下：

model

{

for ( k in 1:K ) {

firm[k] <- data[k, p + 1]

t<-data[k,p+2]

mu[k] <- alpha + exp(eta*(t-T))*u[-firm[k]] + inprod(beta[],data[k,1:p])

y[k] ~ dnorm(mu[k], prec)

}

for (i in 1:N) {

eta ~ dnorm(0.0,4)

u ~ dexp(lambda)

eff <-exp(- u)

}

lambda0 <--log(rstar)

lambda ~ dexp(lambda0)

alpha~ dnorm(0.0, 1.0E-6)

for (i in 1:p) {

beta ~ dnorm(0.0, 1.0E-6)

}

prec~ dgamma(0.001, 0.001)

sigmasq <-1 / prec

}
首先，如何将它从成本函数变为随机前沿生产函数呢？，按照Griffin与Steel的原文，我已经将技术效率项的符号改过来了，其他还需要如何改动呢；第二点，与楼主一样，我不知道该如何将数据导入，因为是面板数据，解释变量有好几个，这样就成了三维的数据了。现在困在这里了，还请论坛里的高人给予指点啊。