潜在クラス分析[Stata]

Stataで潜在クラス分析ができるのでやってみた。
Stata15から標準機能で利用できるらしい。

Lanzaらのプラグインを使えばStata11から利用可能である。

www.methodology.psu.edu

Stataの標準機能だと機能が足りないように思うので、Lanzaらのパッケージを結局使うことになりそうだ。
LanzaらのパッケージはSASでも利用できる。SASといっても、無料で使えるSAS University Editionには対応していない。
ソフトウェアの価格の面から考えるとStataでの利用となりそうだ。

以下で実行しているデモはStata15から搭載されている標準機能を使ったものである。

サンプルデータ

Stataの用意しているサンプルデータを利用する。
http://www.stata-press.com/data/r16/gsem_lca1.dta

変数	ラベル
accident	would testify against friend in accident case 事故の場合に友人に対して証言するだろう
play	would give negative review of friend's play 友人の遊びに否定的な評価をするだろう
insurance	would disclose health concerns to friend's insurance company 友人の健康上の問題を保険会社に開示するだろう
stock	would keep company secret from friend 友人から会社の秘密を守るだろう

何のデータがいまいちわからないが、デモであれば特に問題ないだろう。

潜在クラス分析

Stataでは潜在クラスの式は最初に走らせる必要がある。とりあえず3クラスモデルで走らせている。

. use http://www.stata-press.com/data/r16/gsem_lca1.dta, clear
. gsem (accident play insurance stock <- ), logit lclass(C 3)

コマンドはgsemである。一般化SEMのバリエーションの一つだ。フォーマットは、()の中に潜在クラスに利用する変数を入れて最後に<-をつける。次に、コンマを打ちlogit lclass(C 3)と潜在クラスのオプションをつける。(C 3)は３クラスモデルを意味している。

クラス構成割合

クラス構成割合のコマンドと表示は以下である。

. estat lcprob

Latent class marginal probabilities             Number of obs     =        216

--------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
           C |
          1  |   .1560766   15.31178      2.1e-100           1
          2  |   .6346788   11.94001      2.55e-44           1
          3  |   .2092446   3.373463      1.17e-18           1
--------------------------------------------------------------

応答確率

条件付き応答確率は次のように出力する。

. estat lcmean

Latent class marginal means                     Number of obs     =        216

--------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
1            |
    accident |   .3935409   25.44072      1.20e-91           1
        play |   .0466199   12.58716      4.3e-243           1
   insurance |   .2541585   3.758665      4.52e-18           1
       stock |   .0656737   2.330868      3.26e-34           1
-------------+------------------------------------------------
2            |
    accident |   .8220824   4.909919      1.23e-28           1
        play |    .455712   9.045093      7.63e-32           1
   insurance |   .4224278    5.13116      9.18e-19           1
       stock |   .1818264    3.94581      5.89e-24           1
-------------+------------------------------------------------
3            |
    accident |   .9963741   .3046207      4.65e-70           1
        play |   .9725123   1.448169      2.73e-45           1
   insurance |   .9850408   3.020941      2.1e-173           1
       stock |    .881904   6.083824      1.41e-49           1
--------------------------------------------------------------

フィッティング指標

潜在クラス分析ではクラス数を検討しなければならない。 estat lcgofコマンドで尤度比、尤度比検定、AIC、BICが出力できる。

. estat lcgof
----------------------------------------------------------------------------
Fit statistic        |      Value   Description
---------------------+------------------------------------------------------
Likelihood ratio     |
          chi2_ms(1) |      0.387   model vs. saturated
            p > chi2 |      0.534
---------------------+------------------------------------------------------
Information criteria |
                 AIC |   1034.602   Akaike's information criterion
                 BIC |   1081.856   Bayesian information criterion
----------------------------------------------------------------------------

フィッティング指標の比較

. gsem (accident play insurance stock <- ), logit lclass(C 2)
estimates store twoclass
. gsem (accident play insurance stock <- ), logit lclass(C 3)
estimates store threeclass

それぞれの結果をtwoclassとthreeclassに格納する。この値は任意である。
この２つのフィッティング指標を比較することができる。

. estimates stats twoclass threeclass
Akaike's information criterion and Bayesian information criterion

-----------------------------------------------------------------------------
       Model |          N   ll(null)  ll(model)      df        AIC        BIC
-------------+---------------------------------------------------------------
    twoclass |        216          .  -504.4677       9   1026.935   1057.313
  threeclass |        216          .  -503.3011      14   1034.602   1081.856
-----------------------------------------------------------------------------
Note: BIC uses N = number of observations. See [R] BIC note.

AICとBICをみるといずれも2クラスモデルの方が値が小さいので、2クラスモデルの方がさそうである。

ベイジアンロジスティック回帰分析・その３

HardinとHilbeのStata本にあるベイジアンロジスティック回帰分析の続き。今回は事前分布としてコーシー分布を使う方法についてである。

１回目エントリはこちら
２回目エントリはこちら

Generalized Linear Models and Extensions: Fourth Edition

情報事前分布

先の２つのエントリーでは無情報事前分布(non-informative prior distribution)を用いてきたが、今回は情報事前分布(informative priors)を利用する。

センタリングした変数は医師の訪問cdocと年齢cageの２つで、この２つの変数をコーシー事前分布(Cauchy priors)を用いて推定する。ベイズ事前分布は堅実な外部情報に基づいている必要がある。ベイズ推定のように、既に同じようなものを計測していると、より確かな推定値が得られる。コーシー分布が堅実だと言い切れるかというとそうではない。実際に推定して推定結果を比較してどちらが良いか、比較することになる。

下記のケースでは推定の際に、任意の値をspecificationととして用いることで、デフォルトとは異なる推定を行う。 コーシー分布の理論的な説明は後ほどするとして、先にStataでの使用方法について述べていこう。

コーシー分布をベイジアンロジスティック回帰

Stata15以降の機能である。

. bayes,
> prior({outwork:female kids _cons}, normal(0, 1000))
> prior({outwork:cdoc}, cauchy(0,6))
> prior({outwork:cage}, cauchy(0,11))
> rseed (35841): glm outwork female kids cdoc cage, family(binomial)

Model summary
------------------------------------------------------------------------------
Likelihood:
  outwork ~ glm(xb_outwork)

Priors:
  {outwork:female kids _cons} ~ normal(0,1000)                             (1)
               {outwork:cdoc} ~ cauchy(0,6)                                (1)
               {outwork:cage} ~ cauchy(0,11)                               (1)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outwork.

Bayesian generalized linear models               MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
Family : Bernoulli                               Number of obs    =      3,874
Link   : logit                                   Scale parameter  =          1
                                                 Acceptance rate  =      .2057
                                                 Efficiency:  min =     .03496
                                                              avg =     .04486
Log marginal-likelihood = -1993.0855                          max =     .05717

------------------------------------------------------------------------------
             |                                                Equal-tailed
     outwork |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
      female |  2.260141   .0842407   .004506   2.254865   2.097708   2.431913
        kids |  .3540421   .0906731   .003792   .3560368   .1811122   .5379961
        cdoc |  .0248758   .0063243   .000318   .0247956   .0130655   .0382379
        cage |  .0544032   .0039354    .00017   .0543025   .0469003   .0626092
       _cons | -2.012612   .0799112   .004051  -2.011676  -2.171733  -1.862155
------------------------------------------------------------------------------

コーシー分布を事前分布ベイジアンロジスティック回帰・MH法

同じく、Stata15以降の機能である。

. bayesmh outwork female kids cdoc cage, likelihood (logit)
> prior({outwork:female kids _cons}, normal(0, 1000))
> prior({outwork:cdoc}, cauchy(0,6))
> prior({outwork:cage}, cauchy(0,11))
> rseed (35841)

Model summary
------------------------------------------------------------------------------
Likelihood:
  outwork ~ logit(xb_outwork)

Priors:
  {outwork:female kids _cons} ~ normal(0,1000)                             (1)
               {outwork:cdoc} ~ cauchy(0,6)                                (1)
               {outwork:cage} ~ cauchy(0,11)                               (1)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outwork.

Bayesian logistic regression                     MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
                                                 Number of obs    =      3,874
                                                 Acceptance rate  =      .1844
                                                 Efficiency:  min =     .03402
                                                              avg =     .05568
Log marginal-likelihood = -1993.0309                          max =     .07315

------------------------------------------------------------------------------
             |                                                Equal-tailed
     outwork |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
      female |  2.258607   .0839686   .004552   2.261023     2.0928   2.423481
        kids |  .3565683   .0899923   .003727   .3605532   .1768716   .5279933
        cdoc |   .024811   .0064852   .000282   .0244694   .0122816   .0378834
        cage |  .0547071   .0039411   .000146   .0547988   .0469502     .06197
       _cons | -2.010109   .0794933   .003246  -2.009283  -2.171436  -1.854177
------------------------------------------------------------------------------

情報基準量

情報基準量はいつもの通り。

. bayesstats ic
Bayesian information criteria

----------------------------------------------
             |       DIC    log(ML)    log(BF)
-------------+--------------------------------
      active |  3928.181  -1993.031          .
----------------------------------------------
Note: Marginal likelihood (ML) is computed
      using Laplace-Metropolis approximation.

分布の比較

無情報事前分布(noninformative)とコーシー分布(informative)を用いた推定値の比較をしてみよう。

	noninformative	informative
female	2.60	2.259
kids	0.356	0.356
cdoc	0.024	0.025
cage	0.055	0.055
xcons	-2.01	-2.01

Stataにはモデルの事後確率を基に仮説検定ができるコマンドがあるとのこと。bayestest modelコマンドを利用すると、モデルの限界尤度と事前確率と事後確率が表示できるようだ。レファレンスを参照とのこと。

プロット

cageのグラフに改善点があったため、それが前回の推定より改善されているかが問題となる。改善をしていれば、コーシー事前分布は有用ということになる。

. bayesgraph diagnostics {outwork:cdoc}

. bayesgraph diagnostics {outwork:cage}

このグラフは以前よりも改善されている。センタリングされた年齢の変数は以前のグラフではほぼ二峰性だった(参照)。 また、このモデルでは、無情報事前分布のモデルよりもcageの確信区間が狭くなっている。その理由は事前確率があった方が情報量がわずかに多いためである。Cauchy事前分布のモデルは、フラット事前分布のモデルよりも、センタリングされた年齢変数に適したものになっている。

コーシー分布

コーシー事前分布では、位置パラメーターとスケールパラメーター指定する。

上の式をStataのコマンドにすると以下のようになる。

prior({outwork: cage}, logdensity(ln(11)-ln(11^2+({outwork: cage})^2-ln(_pi))))

Stata15から導入されたprefixであるcauchy()はこの式をlogdensity()オプションに入れなくてもコマンドを走らせることができる。

ベイズ統計でコーシー分布といえばGelmanの論文である。

projecteuclid.org PDF
http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf

本としては以下のもの。

Bayesian Data Analysis (Chapman & Hall/CRC Texts in Statistical Science)

コーシー分布を事前分布としたベイジアンロジスティック回帰

logdensity()のオプション指定をすると、Stata15から搭載されたcauchy()のprefixを使わずに、同様の意味を持つ式が書けるのでStata14でも走るように思うのだが、Stata14は持ってないのでよくわからない。

変数cdocおよびcageについては、a = 0を設定する。cdocはscaleパラメーターb = 6を設定し、cageはb = 11とする。説明が前後しているが、これらの任意の値は冒頭のスクリプトにも入れてある。変数femaleとkidsはnormal(0,1000）の事前分布を使用。これはベイズ統計を使用するときのデフォルトである。

スクリプトは次のようになる。

. bayes, prior({outwork:female kids _cons}, normal(0, 1000)) 
> prior({outwork:cdoc},
> logdensity(ln(6)-ln(6^2+({outwork:cdoc})^2)-ln(_pi)))
> prior({outwork:cage},
> logdensity(ln(11)-ln(11^2+{outwork:cage}^2)-ln(_pi)))
> rseed (35841): glm outwork female kids cdoc cage, family(binomial)

コーシー分布を事前分布としたベイジアンロジスティック回帰・MH法

MH法のbayesmhコマンドでは次のように指定する。

. bayesmh outwork female kids cdoc cage, likelihood(logit)
> prior({outwork:female kids _cons}, normal(0, 1000))
> prior({outwork:cdoc},
> logdensity(ln(6)-ln(6^2+({outwork:cdoc})^2)-ln(_pi)))
> prior({outwork:cage},
> logdensity(ln(11)-ln(11^2+({outwork:cage})^2)-ln(_pi)))
> rseed (35841)

logdensity()でコーシー分布の指定をしなくてもいいので、Stata15以降の方が多少楽になった気はする。
HardinとHilbeのテキストのコーシー分布自力計算のスクリプト内logdensity()オプションの中に誤植があったので、教科書を書く人間が間違うくらいだから、Stata15以降の改善は良いことなのだろうと思う。

ベイジアンロジスティック回帰分析・その２

HardinとHilbeのStata本にあるベイジアンロジスティック回帰分析の続き。前回のエントリはこちら

Generalized Linear Models and Extensions: Fourth Edition

メトロポリス・ヘイスティングス法

メトロポリス・ヘイスティングス法のやり方である。

前回も同じことを行っているが、データの読み込みと変数加工(センタリング)もう一度行っておこう。

. use http://www.stata-press.com/data/hh4/rwm1984
. quietly summarize docvis, meanonly
. generate cdoc = docvis - r(mean)
. quietly summarize age, meanonly
. generate cage = age - r(mean)

医者の往診の変数docvisと年齢ageをセンタリングしている。

MH法では変数の設定が必要である。priorは従属変数を指定し、その後にコンマをうち、分布とパラメーターを指定する構文が必要だ。以下のコマンドの2行目のところである。

. bayesmh outwork female kids cdoc cage, likelihood(logit)
> prior({outwork:}, normal (0,10000)) rseed (2017)
Model summary
------------------------------------------------------------------------------
Likelihood:
  outwork ~ logit(xb_outwork)

Prior:
  {outwork:female kids cdoc cage _cons} ~ normal(0,10000)                  (1)
------------------------------------------------------------------------------
(1) Parameters are elements of the linear form xb_outwork.

Bayesian logistic regression                     MCMC iterations  =     12,500
Random-walk Metropolis-Hastings sampling         Burn-in          =      2,500
                                                 MCMC sample size =     10,000
                                                 Number of obs    =      3,874
                                                 Acceptance rate  =      .2363
                                                 Efficiency:  min =     .02665
                                                              avg =     .05244
Log marginal-likelihood = -2001.0895                          max =       .101

------------------------------------------------------------------------------
             |                                                Equal-tailed
     outwork |      Mean   Std. Dev.     MCSE     Median  [95% Cred. Interval]
-------------+----------------------------------------------------------------
      female |  2.260566   .0803699   .002529    2.25819   2.103102   2.421029
        kids |  .3562706   .0893508   .004221   .3583318   .1783027   .5281715
        cdoc |   .024483   .0061758   .000329   .0244374   .0130205   .0369231
        cage |  .0545722   .0040688   .000174   .0547603   .0464829   .0622289
       _cons |  -2.01027    .082349   .005044   -2.00927  -2.172552  -1.852391
------------------------------------------------------------------------------

プロット

プロットは非収束の徴候が確認できる。

従属変数と各独立変数でプロットを作成する。
従属変数であるoutworkとcdocでグラフを作成した。

Stataでは引数なしでグラフ診断を走らせると、４タイプのグラフが描画される。

. bayesgraph diagnostics {outwork: cage}

トレースプロットは各反復におけるパラメータ推定値をグラフ化する。ラグ項との自己相関である。ヒストグラムは事前分布と一致しているか確認する。outworkとcdocでのグラフ診断は正規分布に近いが、世正規分布が必須というわけではないようだ。

outworkとセンタリングした年齢の変数であるcageのグラフ診断も走らせてみる。

cdocほどは良いグラフになっていないので、サンプリングに反復を追加し、値の選択を3から5減らすと自己相関を減らすことができるとのこと。

ESSコマンド

ESSコマンドでは、モデル内の指定された各パラメータについて、有効なサンプルサイズ、相関時間、および効率を表示できる。

. bayesstats ess

Efficiency summaries    MCMC sample size =    10,000
                        Efficiency:  min =    .02665
                                     avg =    .05244
                                     max =      .101

----------------------------------------------------
    outwork  |        ESS   Corr. time    Efficiency
-------------+--------------------------------------
      female |    1009.74         9.90        0.1010
        kids |     448.19        22.31        0.0448
        cdoc |     352.37        28.38        0.0352
        cage |     545.17        18.34        0.0545
       _cons |     266.54        37.52        0.0267
----------------------------------------------------

ベイズモデルでは、モデルパラメータの事後平均は、パラメータの事後分布のMCMCサンプリングの平均によって決定される。推定されたサンプル平均の標準誤差はに比例する。

ESS値が低いと注意が必要である。理想的には、ESS値はMCMCサンプルサイズに近い値をとる。効率は、ESS値にまたはを乗算することによって決定される。女性のESS値は1009.74であり、効率は0.101である。

Gibbs法によるサンプリングを行わない限り、平均の効率値0.01が許容できるモデルの最小効率である。0.1以上の値はランダムウォーク・メトロポリス・アルゴリズム(random-walk Metropolis algorithm)にとって効率的であるとされているとのこと。

相関時間(correlation time)は効率の逆数である。自己相関が低い推定遅延と考えられ、値は小さい方が良い。

献本御礼：宮本みち子編『すべての若者が生きられる未来を』

すべての若者が生きられる未来を――家族・教育・仕事からの排除に抗して

宮本先生に講演会でお会いして頂戴いたしました。
ありがとうございます。

www.city.shinagawa.tokyo.jp

インドでひきこもり、いや、Netflix依存症だ

dramanavi.net

米CNBCによると、昨年度Netflix加入者の1日平均視聴時間は50分。しかし、インド出身のある男性は、1日7時間もNetflixの作品を見続けるという生活を半年以上も続けた結果、睡眠障害や身体的疲労、眼精疲労などの症状が出てしまい、インドのバンガロールにある同施設に入所したという。

施設にある、テクノロジーの健全な使用推進を奨励している部署の精神科医、マノジュ・クマール・シャルマ医師によると、その男性は半年以上無職で、家族から職を見つけるようプレッシャーをかけられており、その重圧から逃れるためにNetflixで視聴を続けていたとのこと。

この男性は日本でいうとひきこもりであったり、精神疾患の分類ではうつ病あたりが妥当なところだろう。
Netflixを見たため睡眠障害や身体的疲労などの症状が出たというのではなく、ストレス下でそれらの精神症状が出たと捉える方が説明として無理がない。
Netflixは、ストレスから逃れるため、つまりコーピングの手段として使われていたのだから、Netflixはむしろ彼を救っていた部分の方が多いのではないだろうか。

有病率の計算[R]

こちらのエントリでStataで行ったことをRでもやってみたい。 有病率の計算はprevalenceパッケージを用いる。

prevalence.cbra.be

prevalenceパッケージはベイズ統計学の計算にJAGSとrjagsを利用する。
JAGSを使ったことがない場合には、インストールが必要である。

JAGSのインストール

sourceforgeからダウンロードする。

sourceforge.net

ほとんどの人は/files/JAGS/4.x/Windows/からダウンロードすることになるだろう。
通常のアプリケーションと同じようにダウンロードをして、ダイアログボックスに従ってインストールする。

ほとんどの人はJAGS/4.x/Windows/からダウンロードすることになるだろう。ちなみに僕は/files/JAGS/4.x/Windows/JAGS-4.3.0.exeをインストールした。
通常のアプリケーションと同じようにダウンロードをして、ダイアログボックスに従ってインストールする。Pathは通す必要がない。

rjagsパッケージのインストール

rjagsはRからインストールする。

install.packages("rjags")

prevalenceパッケージのインストール

いつもは書かないが、一応prevalenceパッケージのインストールもコードを書いておこう。

install.packages("prevalence")

内閣府ひきこもり調査の例

生活状況に関する調査 https://www8.cao.go.jp/youth/kenkyu/life/h30/pdf/kekka_gaiyo.pdf

この調査の結果は下記であった。
- 有効回答数 3248人 - ひきこもり群 47人

library(prevalence)
propCI(47, 3248, method = "all", level = 0.95, sortby = "level")

結果は以下。

   x    n          p        method level      lower      upper
1 47 3248 0.01447044 agresti.coull  0.95 0.01086018 0.01922783
2 47 3248 0.01447044         exact  0.95 0.01065120 0.01919662
3 47 3248 0.01447044      jeffreys  0.95 0.01078351 0.01902143
4 47 3248 0.01447044          wald  0.95 0.01036353 0.01857736
5 47 3248 0.01447044        wilson  0.95 0.01089964 0.01918838

Rのprevalenceパッケージの場合はどの推定方法であっても、信頼区間のレベルを指定できるとこがStataより優れている。Stataはすべての推定法で区間レベルの変更がオプションで用意されていない。
通常は、95%しか使わないので、不要かもしれないが、かゆいところに手が届くのがRおよび有志の作ったパッケージのよいところである。

推定方法が決まっている場合にはmethodで指定する。agresti.coull法を使うのであればmethod = agresti.coullにする。

ちなみに、僕の環境ではpropCI()はRStudioのコンソールからしか実行できなかった。 いつもの回避法だが、R Markdownでは、一度結果を格納してからprint()をすると表示できる。

res<-propCI(47, 3248, method = "all", level = 0.95, sortby = "level")
print(res)

各推定方法についてはStataのエントリを参考にして欲しい → こちら

有病率の計算[Stata]

疫学調査で有病率を出す方法である。専門用語では母比率の推定という。日本語で母比率と言うとピンとこないが、英語だとpopulation rateなので、こちらの方が直感的に理解できるはずである。

Stataで有病率を推定する

Stataではcii propコマンドを使用する。

データは内閣府のひきこもり調査のものを利用しよう。

生活状況に関する調査 https://www8.cao.go.jp/youth/kenkyu/life/h30/pdf/kekka_gaiyo.pdf

この調査の結果は下記であった。

• 有効回答者数 3248人
• ひきこもり者数 47人

. cii prop 3248 47 , level(95)

                                                         -- Binomial Exact --
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
             |      3,248    .0144704    .0020954        .0106512    .0191966

日本の満40歳から満64歳までの人口は4235万人である。それぞれの値を人口に掛けると推定値と95%信頼区間が出てくる。

Proportion: 42350000 * .0144704 = 612821.4
95% CI Upper: 42350000 * .0106512 = 451078.3
95% CI Lower: 42350000 * .0191966 = 812976

よって推定値61.3万人、95%信頼区間45.1～81.3万人となる。

その他の推定方法

母比率の推定にはいくつか方法がある。

コマンド	説明
exact	calculate exact confidence intervals; the default
wald	calculate Wald confidence intervals
wilson	wilson calculate Wilson confidence intervals
agresti	calculate Agresti–Coull confidence intervals
jeffreys	calculate Jeffreys confidence intervals

clopper-pearson正確法

デフォルトの設定で有病率である。F分布を使用する。 スクリプトは上段で示したものである。

wald法

本やウェブに掲載されているのはだいたいこの推定方法である。
wald法は正規分布を使用する。Stataでは wald法はデフォルトではないので、注意が必要かもしれない。

. cii prop 3248 47, wald

                                                         -- Binomial Wald ---
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
             |      3,248    .0144704    .0020954        .0103635    .0185774

wilsonのスコア法

. cii prop 3248 47, wilson

                                                         ------ Wilson ------
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
             |      3,248    .0144704    .0020954        .0108996    .0191884

agresti-coull法

. cii prop 3248 47, agresti

                                                         -- Agresti-Coull ---
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
             |      3,248    .0144704    .0020954        .0108602    .0192278

元論文
https://amstat.tandfonline.com/doi/abs/10.1080/00031305.1998.10480550#.Xbk-muj7Shc

本文(PDF)
http://www.uvm.edu/~rsingle/stat380/F04/possible/Agresti+Coull-Amstat-1998_ApproxVsExactCIfoP.pdf

jeffrey's法

JeffreyとはRichard Jeffreyのことでベイズ的な方法とされる。

定義式を見るとclopper-pearson法の改良版と位置付けてもよさそうだ。

. cii prop 3248 47, jeffreys

                                                         ----- Jeffreys -----
    Variable |        Obs  Proportion    Std. Err.       [95% Conf. Interval]
-------------+---------------------------------------------------------------
             |      3,248    .0144704    .0020954        .0107835    .0190214

考察

これらの方法は95%信頼区間の計算方法なので、比率(Stataの結果Proportionのところ、wald法の中央値)には影響を与えない。

例題のように一般人口を対象にしたサンプルサイズが十分に大きい調査で、有病率を計算する場合、agresti-coull法がreasonableであろう。

母比率の計算は様々なケースで行うので、当然、分布も様々なものが想定される。その際には正規分布よりF分布の方が適切という場合もあるだろうから、どんな場合でもagresti-coull法がよいというわけではない。その場合はF分布のclopper-pearson正確法かjeffrey's法を検討することになるだろう。ただ、だいたいのケースではagresti-coull法が最適なのではないかと思う。