🕷️ Crawler Inspector

URL Lookup

Direct Parameter Lookup

Raw Queries and Responses

1. Shard Calculation

Query:

Response:

Calculated Shard: 143 (from laksa000)

2. Crawled Status Check

Query:

curl -X POST \
  'http://laksa143.int.ahrefs:8124/' \
  -H 'Content-Type: text/plain' \
  -H 'X-ClickHouse-Database: crawler3' \
  -H 'Authorization: Basic YXBpOg==' \
  -d 'SELECT getAhrefsURLFromUnparsed(src_unparsed) AS found_url, ifNull(toUnixTimestamp(download_stamp), 0) AS crawl_time, ifNull(toUnixTimestamp(props_url_first_seen), 0) AS first_indexed_time, download_http_code AS http_code, src_unparsed AS src_unparsed, src_root_hash AS src_root_hash, history_drop_reason AS history_drop_reason, meta_title AS meta_title, meta_descriptions AS meta_descriptions, attrs_boilerpipe_text AS attrs_boilerpipe_text, attrs_markdown AS attrs_markdown, attrs_readable_markdown AS attrs_readable_markdown, meta_canonical AS meta_canonical FROM crawler3.page_info_local FINAL PREWHERE (src_root_hash, src_unparsed) IN ((getAhrefsRootHashFromUnparsed(getAhrefsUnparsedNoserviceFromURL(\'https://m-clark.github.io/models-by-example/zi.html\')), getAhrefsUnparsedNoserviceFromURL(\'https://m-clark.github.io/models-by-example/zi.html\'))) FORMAT JSONEachRow'

Response:

{"found_url":"https:\/\/m-clark.github.io\/models-by-example\/zi.html","crawl_time":1773948913,"first_indexed_time":1605125796,"http_code":200,"src_unparsed":"io,github!m-clark,\/models-by-example\/zi.html s443","src_root_hash":"2566890010099092343","history_drop_reason":null,"meta_title":"Zero-Inflated Model | Model Estimation by Example","meta_descriptions":["This document provides ‘by-hand’ demonstrations of various models and algorithms. The goal is to take away some of the mystery by providing clean code examples that are easy to run and compare with other tools."],"attrs_boilerpipe_text":"Log likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples\nand comparison to\npscl\npackage output. Also includes approach based on Hilbe GLM text.\nZero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory.\nPoisson\nData Setup\nGet the data.\nlibrary\n(tidyverse)\nfish\n=\nhaven\n::\nread_dta\n(\n\"http:\/\/www.stata-press.com\/data\/r11\/fish.dta\"\n)\nFunction\nThe log likelihood function.\nzip_ll\n<-\nfunction\n(y, X, par) {\n# arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois'\n# Extract parameters\nlogitpars\n=\npar[\ngrep\n(\n'logit'\n,\nnames\n(par))]\npoispars\n=\npar[\ngrep\n(\n'pois'\n,\nnames\n(par))]\n# Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example\nXlogit\n=\nX\nLPlogit\n=\nXlogit\n%*%\nlogitpars\nlogi0\n=\nplogis\n(LPlogit)\n# alternative 1\/(1+exp(-LPlogit))\n# Poisson part\nXpois\n=\nX\nmupois\n=\nexp\n(Xpois\n%*%\npoispars)\n# LLs\nlogliklogit\n=\nlog\n( logi0\n+\nexp\n(\nlog\n(\n1\n-\nlogi0)\n-\nmupois) )\nloglikpois\n=\nlog\n(\n1\n-\nlogi0)\n+\ndpois\n(y,\nlambda =\nmupois,\nlog =\nTRUE\n)\n# Hilbe formulation\n# logliklogit = log(logi0 + (1 - logi0)*exp(- mupois) )\n# loglikpois = log(1-logi0) -mupois + log(mupois)*y     #not necessary: - log(gamma(y+1))\ny0\n=\ny\n==\n0\n# 0 values\nyc\n=\ny\n>\n0\n# Count part\nloglik\n=\nsum\n(logliklogit[y0])\n+\nsum\n(loglikpois[yc])\n-\nloglik\n}\nEstimation\nGet starting values or simply do zeros.\n# for zip: need 'logit', 'pois'\ninitial_model\n=\nglm\n(\ncount\n~\npersons\n+\nlivebait,\ndata =\nfish,\nx =\nTRUE\n,\ny =\nTRUE\n,\n\"poisson\"\n)\n# starts = c(logit = coef(initial_model), pois = coef(initial_model))\nstarts\n=\nc\n(\nrep\n(\n0\n,\n3\n),\nrep\n(\n0\n,\n3\n))\nnames\n(starts)\n=\nc\n(\npaste0\n(\n'pois.'\n,\nnames\n(\ncoef\n(initial_model))),\npaste0\n(\n'logit.'\n,\nnames\n(\ncoef\n(initial_model))))\nEstimate with\noptim\n.\nfit\n=\noptim\n(\npar =\nstarts ,\nfn  =\nzip_ll,\nX   =\ninitial_model\n$\nx,\ny   =\ninitial_model\n$\ny,\nmethod  =\n\"BFGS\"\n,\ncontrol =\nlist\n(\nmaxit =\n5000\n,\nreltol =\n1e-12\n),\nhessian =\nTRUE\n)\n# fit\nComparison\nExtract for clean display.\nB\n=\nfit\n$\npar\nse\n=\nsqrt\n(\ndiag\n(\nsolve\n((fit\n$\nhessian))))\nZ\n=\nB\n\/\nse\np\n=\npnorm\n(\nabs\n(Z),\nlower =\nFALSE\n)\n*\n2\nResults from\npscl\n.\nlibrary\n(pscl)\nfit_pscl\n=\nzeroinfl\n(count\n~\npersons\n+\nlivebait,\ndata =\nfish,\ndist =\n\"poisson\"\n)\nCompare.\ncoef\nB\nse\nZ\np\npscl\nX.Intercept.\n-2.006\n0.324\n-6.196\n0.000\npscl\npersons\n0.747\n0.043\n17.516\n0.000\npscl\nlivebait\n1.809\n0.292\n6.195\n0.000\npscl\nX.Intercept..1\n0.303\n0.674\n0.449\n0.654\npscl\npersons.1\n-0.069\n0.129\n-0.537\n0.591\npscl\nlivebait.1\n-0.031\n0.558\n-0.056\n0.956\nzip_poisson_ll\npois.(Intercept)\n-2.006\n0.324\n-6.196\n0.000\nzip_poisson_ll\npois.persons\n0.747\n0.043\n17.516\n0.000\nzip_poisson_ll\npois.livebait\n1.809\n0.292\n6.195\n0.000\nzip_poisson_ll\nlogit.(Intercept)\n0.303\n0.674\n0.449\n0.654\nzip_poisson_ll\nlogit.persons\n-0.069\n0.129\n-0.537\n0.591\nzip_poisson_ll\nlogit.livebait\n-0.031\n0.558\n-0.056\n0.956\nNegative Binomial\nFunction\nzinb_ll\n<-\nfunction\n(y, X, par) {\n# arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta'\n# Extract parameters\nlogitpars\n=\npar[\ngrep\n(\n'logit'\n,\nnames\n(par))]\nnegbinpars\n=\npar[\ngrep\n(\n'negbin'\n,\nnames\n(par))]\ntheta\n=\nexp\n(par[\ngrep\n(\n'theta'\n,\nnames\n(par))])\n# Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example\nXlogit\n=\nX\nLPlogit\n=\nXlogit\n%*%\nlogitpars\nlogi0\n=\nplogis\n(LPlogit)\n# Negbin part\nXnegbin\n=\nX\nmunb\n=\nexp\n(Xnegbin\n%*%\nnegbinpars)\n# LLs\nlogliklogit\n=\nlog\n(\nlogi0\n+\nexp\n(\nlog\n(\n1\n-\nlogi0)\n+\nsuppressWarnings\n(\ndnbinom\n(\n0\n,\nsize =\ntheta,\nmu   =\nmunb,\nlog  =\nTRUE\n))\n)\n)\nlogliknegbin\n=\nlog\n(\n1\n-\nlogi0)\n+\nsuppressWarnings\n(\ndnbinom\n(y,\nsize =\ntheta,\nmu   =\nmunb,\nlog  =\nTRUE\n))\n# Hilbe formulation\n# theta part\n# alpha = 1\/theta\n# m = 1\/alpha\n# p = 1\/(1 + alpha*munb)\n# logliklogit = log( logi0 + (1 - logi0)*(p^m) )\n# logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + m*log(p) + y*log(1-p)   # gamma(y+1) not needed\ny0\n=\ny\n==\n0\n# 0 values\nyc\n=\ny\n>\n0\n# Count part\nloglik\n=\nsum\n(logliklogit[y0])\n+\nsum\n(logliknegbin[yc])\n-\nloglik\n}\nEstimation\nGet starting values or simply do zeros.\n# for zinb: 'logit', 'negbin', 'theta'\ninitial_model\n=\nmodel.matrix\n(count\n~\npersons\n+\nlivebait,\ndata =\nfish)\n# to get X matrix\nstartlogi\n=\nglm\n(count\n==\n0\n~\npersons\n+\nlivebait,\ndata =\nfish,\nfamily =\n\"binomial\"\n)\nstartcount\n=\nglm\n(count\n~\npersons\n+\nlivebait,\ndata =\nfish,\nfamily =\n\"poisson\"\n)\nstarts\n=\nc\n(\nnegbin =\ncoef\n(startcount),\nlogit =\ncoef\n(startlogi),\ntheta =\n1\n)\n# starts = c(negbin = rep(0, 3),\n#            logit = rep(0, 3),\n#            theta = log(1))\nEstimate with\noptim\n.\nfit_nb\n=\noptim\n(\npar =\nstarts ,\nfn  =\nzinb_ll,\nX   =\ninitial_model,\ny   =\nfish\n$\ncount,\nmethod  =\n\"BFGS\"\n,\ncontrol =\nlist\n(\nmaxit =\n5000\n,\nreltol =\n1e-12\n),\nhessian =\nTRUE\n)\n# fit_nb\nComparison\nExtract for clean display.\nB\n=\nfit_nb\n$\npar\nse\n=\nsqrt\n(\ndiag\n(\nsolve\n((fit_nb\n$\nhessian))))\nZ\n=\nB\n\/\nse\np\n=\npnorm\n(\nabs\n(Z),\nlower =\nFALSE\n)\n*\n2\nResults from\npscl\n.\n# pscl results\nlibrary\n(pscl)\nfit_pscl\n=\nzeroinfl\n(count\n~\npersons\n+\nlivebait,\ndata =\nfish,\ndist =\n\"negbin\"\n)\nCompare.\ncoef\nB\nse\nZ\np\npscl\nX.Intercept.\n-2.803\n0.558\n-5.026\n0.000\npscl\npersons\n0.849\n0.124\n6.833\n0.000\npscl\nlivebait\n1.791\n0.511\n3.504\n0.000\npscl\nLog.theta.\n-0.969\n0.302\n-3.206\n0.001\npscl\nX.Intercept..1\n-4.276\n4.278\n-1.000\n0.318\npscl\npersons.1\n0.560\n0.517\n1.084\n0.279\npscl\nlivebait.1\n1.168\n3.661\n0.319\n0.750\nzinb_ll\nnegbin.(Intercept)\n-2.803\n0.558\n-5.026\n0.000\nzinb_ll\nnegbin.persons\n0.849\n0.124\n6.834\n0.000\nzinb_ll\nnegbin.livebait\n1.791\n0.511\n3.504\n0.000\nzinb_ll\nlogit.(Intercept)\n-4.276\n4.278\n-1.000\n0.318\nzinb_ll\nlogit.persons\n0.560\n0.517\n1.084\n0.279\nzinb_ll\nlogit.livebait\n1.168\n3.661\n0.319\n0.750\nzinb_ll\ntheta\n-0.969\n0.302\n-3.206\n0.001","attrs_markdown":"- [Model Estimation by Example](https:\/\/m-clark.github.io\/models-by-example\/)\n- [Introduction](https:\/\/m-clark.github.io\/models-by-example\/introduction.html)\n- **Models**\n- [Linear Regression](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html#data-setup)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html#functions)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html#estimation)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html#comparison)\n  - [Python](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html#python)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/linear-regression.html#source)\n- [Logistic Regression](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html#data-setup-1)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html#functions-1)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html#estimation-1)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html#comparison-1)\n  - [Python](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html#python-1)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/logistic-regression.html#source-1)\n- [Mixed Models](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html)\n  - [One-factor Mixed Model](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#one-factor-mixed-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#data-setup-2)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#function)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#estimation-2)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#comparison-2)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#source-2)\n  - [Two-factor Mixed Model](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#two-factor-mixed-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#data-setup-3)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#function-1)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#estimation-3)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#comparison-3)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#source-3)\n  - [Mixed Model via ML](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#mixed-model-via-ml)\n    - [Introduction](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#introduction-1)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#data-setup-4)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#function-2)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#estimation-4)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#comparison-4)\n    - [Link with penalized regression](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#link-with-penalized-regression)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/mixed-models.html#source-4)\n- [Probit & Bivariate Probit](https:\/\/m-clark.github.io\/models-by-example\/probit.html)\n  - [Standard Probit](https:\/\/m-clark.github.io\/models-by-example\/probit.html#standard-probit)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/probit.html#function-3)\n    - [Examples](https:\/\/m-clark.github.io\/models-by-example\/probit.html#examples)\n  - [Bivariate Probit](https:\/\/m-clark.github.io\/models-by-example\/probit.html#bivariate-probit)\n    - [Example](https:\/\/m-clark.github.io\/models-by-example\/probit.html#example)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/probit.html#source-5)\n- [Heckman Selection](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#data-setup-5)\n  - [Two step approach](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#two-step-approach)\n    - [Step 1: Probit Model](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#step-1-probit-model)\n    - [Step 2: Estimate via Linear Regression](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#step-2-estimate-via-linear-regression)\n  - [Maximum Likelihood](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#maximum-likelihood)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#comparison-5)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/heckman-selection.html#source-6)\n- [Marginal Structural Model](https:\/\/m-clark.github.io\/models-by-example\/marginal-structural.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/marginal-structural.html#data-setup-6)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/marginal-structural.html#function-4)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/marginal-structural.html#estimation-5)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/marginal-structural.html#comparison-6)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/marginal-structural.html#source-7)\n- [Tobit Regression](https:\/\/m-clark.github.io\/models-by-example\/tobit.html)\n  - [Censoring with an Upper Limit](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#censoring-with-an-upper-limit)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#data-setup-7)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#function-5)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#estimation-6)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#comparison-7)\n  - [Censoring with a Lower Limit](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#censoring-with-a-lower-limit)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#comparison-8)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/tobit.html#source-8)\n- [Cox Survival](https:\/\/m-clark.github.io\/models-by-example\/cox.html)\n  - [Standard Proportional Hazards](https:\/\/m-clark.github.io\/models-by-example\/cox.html#standard-proportional-hazards)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/cox.html#data-setup-8)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/cox.html#function-6)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/cox.html#estimation-7)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/cox.html#comparison-9)\n  - [Time-varying coefficients](https:\/\/m-clark.github.io\/models-by-example\/cox.html#time-varying-coefficients)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/cox.html#data-setup-9)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/cox.html#function-7)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/cox.html#estimation-8)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/cox.html#comparison-10)\n  - [Stratified Cox Model](https:\/\/m-clark.github.io\/models-by-example\/cox.html#stratified-cox-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/cox.html#data-setup-10)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/cox.html#function-8)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/cox.html#estimation-9)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/cox.html#comparison-11)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/cox.html#source-9)\n- [Hurdle Model](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html)\n  - [Poisson](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#poisson)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#data-setup-11)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#function-9)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#estimation-10)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#comparison-12)\n  - [Negative Binomial](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#negative-binomial)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#function-10)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#estimation-11)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#comparison-13)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/hurdle.html#source-10)\n- [Zero-Inflated Model](https:\/\/m-clark.github.io\/models-by-example\/zi.html)\n  - [Poisson](https:\/\/m-clark.github.io\/models-by-example\/zi.html#poisson-1)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/zi.html#data-setup-12)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/zi.html#function-11)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/zi.html#estimation-12)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/zi.html#comparison-14)\n  - [Negative Binomial](https:\/\/m-clark.github.io\/models-by-example\/zi.html#negative-binomial-1)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/zi.html#function-12)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/zi.html#estimation-13)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/zi.html#comparison-15)\n  - [Supplemental Example](https:\/\/m-clark.github.io\/models-by-example\/zi.html#supplemental-example)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/zi.html#source-11)\n- [Naive Bayes](https:\/\/m-clark.github.io\/models-by-example\/naive-bayes.html)\n  - [Initialization](https:\/\/m-clark.github.io\/models-by-example\/naive-bayes.html#initialization)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/naive-bayes.html#comparison-16)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/naive-bayes.html#source-12)\n- [Multinomial Regression](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html)\n  - [Standard (Categorical) Model](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#standard-categorical-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#data-setup-13)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#function-13)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#comparison-17)\n  - [Alternative specific and constant variables](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#alternative-specific-and-constant-variables)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#comparison-18)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/multinomial.html#source-13)\n- [Ordinal Regression](https:\/\/m-clark.github.io\/models-by-example\/ordinal.html)\n  - [Data](https:\/\/m-clark.github.io\/models-by-example\/ordinal.html#data)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/ordinal.html#function-14)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/ordinal.html#estimation-14)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/ordinal.html#comparison-19)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/ordinal.html#source-14)\n- [Markov Model](https:\/\/m-clark.github.io\/models-by-example\/markov.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/markov.html#data-setup-14)\n    - [Data Functions](https:\/\/m-clark.github.io\/models-by-example\/markov.html#data-functions)\n    - [Two states Demo](https:\/\/m-clark.github.io\/models-by-example\/markov.html#two-states-demo)\n    - [Three states demo](https:\/\/m-clark.github.io\/models-by-example\/markov.html#three-states-demo)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/markov.html#function-15)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/markov.html#estimation-15)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/markov.html#comparison-20)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/markov.html#source-15)\n- [Hidden Markov Model](https:\/\/m-clark.github.io\/models-by-example\/hmm.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/hmm.html#data-setup-15)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/hmm.html#function-16)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/hmm.html#estimation-16)\n  - [Supplemental demo](https:\/\/m-clark.github.io\/models-by-example\/hmm.html#supplemental-demo)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/hmm.html#source-16)\n- [Quantile Regression](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#data-setup-16)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#function-17)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#estimation-17)\n    - [Other quantiles](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#other-quantiles)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#comparison-21)\n  - [Visualize](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#visualize)\n  - [Python](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#python-2)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/quantile-regression.html#source-17)\n- [Cubic Spline Model](https:\/\/m-clark.github.io\/models-by-example\/cubic-spline.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/cubic-spline.html#data-setup-17)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/cubic-spline.html#functions-2)\n  - [Example 1](https:\/\/m-clark.github.io\/models-by-example\/cubic-spline.html#example-1)\n  - [Example 2](https:\/\/m-clark.github.io\/models-by-example\/cubic-spline.html#example-2)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/cubic-spline.html#source-18)\n- [Gaussian Processes](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html)\n  - [Noise-Free Demonstration](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#noise-free-demonstration)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#data-setup-18)\n    - [Functions](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#functions-3)\n    - [Visualize the prior distribution](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#visualize-the-prior-distribution)\n    - [Generate the posterior predictive distribution](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#generate-the-posterior-predictive-distribution)\n    - [Visualize the Posterior Predictive Distribution](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#visualize-the-posterior-predictive-distribution)\n  - [Noisy Demonstration](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#noisy-demonstration)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#data-setup-19)\n    - [Functions](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#functions-4)\n    - [Visualize the prior distribution](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#visualize-the-prior-distribution-1)\n    - [Generate the posterior predictive distribution](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#generate-the-posterior-predictive-distribution-1)\n    - [Visualize the Posterior Predictive Distribution](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#visualize-the-posterior-predictive-distribution-1)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/gaussian-process.html#source-19)\n- [Neural Network](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html)\n  - [Example 1](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#example-1-1)\n    - [Python](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#python-3)\n    - [R](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#r)\n  - [Example 2](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#example-2-1)\n    - [Python](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#python-4)\n    - [R](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#r-1)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#comparison-22)\n  - [Example 3](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#example-3)\n    - [Python](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#python-5)\n    - [R](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#r-2)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#comparison-23)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/neural-network.html#source-20)\n- [Extreme Learning Machine](https:\/\/m-clark.github.io\/models-by-example\/elm.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/elm.html#data-setup-20)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/elm.html#function-18)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/elm.html#estimation-18)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/elm.html#comparison-24)\n  - [Supplemental Example](https:\/\/m-clark.github.io\/models-by-example\/elm.html#supplemental-example-1)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/elm.html#source-21)\n- [Reproducing Kernel Hilbert Space Regression](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#data-setup-21)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#functions-5)\n    - [Function to find the inverse of a matrix](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#function-to-find-the-inverse-of-a-matrix)\n    - [Reproducing Kernel](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#reproducing-kernel)\n    - [Gram matrix](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#gram-matrix)\n    - [Ridge regression](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#ridge-regression)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#estimation-19)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#comparison-25)\n  - [Example: Cubic Spline](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#example-cubic-spline)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#data-setup-22)\n    - [Functions](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#functions-6)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#estimation-20)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#comparison-26)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#source-22)\n    - [Current Supplemental Code](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#current-supplemental-code)\n    - [Original Supplemental Code](https:\/\/m-clark.github.io\/models-by-example\/rkhs.html#original-supplemental-code)\n- [Confirmatory Factor Analysis](https:\/\/m-clark.github.io\/models-by-example\/cfa.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#data-setup-23)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#functions-7)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#estimation-21)\n    - [Raw](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#raw)\n    - [Standardized](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#standardized)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#comparison-27)\n    - [Mplus](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#mplus)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/cfa.html#source-23)\n- **Bayesian**\n- [Introduction to Bayesian Methods](https:\/\/m-clark.github.io\/models-by-example\/bayesian.html)\n- [Basics](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html)\n  - [Prior, likelihood, & posterior distributions](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html#prior-likelihood-posterior-distributions)\n  - [Prior](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html#prior)\n  - [Likelihood](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html#likelihood)\n  - [Posterior](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html#posterior)\n  - [Posterior predictive](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html#posterior-predictive)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-basics.html#source-24)\n- [Bayesian t-test](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html#data-setup-24)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html#model-code)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html#estimation-22)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html#comparison-28)\n  - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html#visualization)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-t-test.html#source-25)\n- [Bayesian Linear Regression](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html#data-setup-25)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html#model-code-1)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html#estimation-23)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html#comparison-29)\n  - [Visualize](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html#visualize-1)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-linear-regression.html#source-26)\n- [Bayesian Beta Regression](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html#data-setup-26)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html#model-code-2)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html#estimation-24)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html#comparison-30)\n  - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html#visualization-1)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-beta-regression.html#source-27)\n- [Bayesian Mixed Model](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html#data-setup-27)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html#model-code-3)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html#estimation-25)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html#comparison-31)\n  - [Visualize](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html#visualize-2)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-model.html#source-28)\n- [Bayesian Multilevel Mediation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html#data-setup-28)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html#model-code-4)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html#estimation-26)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html#comparison-32)\n  - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html#visualization-2)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-mixed-mediation.html#source-29)\n- [Bayesian IRT](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html)\n  - [One Parameter IRT](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#one-parameter-irt)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#data-setup-29)\n    - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#model-code-5)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#estimation-27)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#comparison-33)\n  - [Two Parameter IRT](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#two-parameter-irt)\n    - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#model-code-6)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#estimation-28)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#comparison-34)\n  - [Three Parameter IRT](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#three-parameter-irt)\n  - [Four Parameter IRT](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#four-parameter-irt)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-irt.html#source-30)\n- [Bayesian CFA](https:\/\/m-clark.github.io\/models-by-example\/bayesian-cfa.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-cfa.html#data-setup-30)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-cfa.html#model-code-7)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-cfa.html#estimation-29)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-cfa.html#comparison-35)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-cfa.html#source-31)\n- [Bayesian Nonparametric Models](https:\/\/m-clark.github.io\/models-by-example\/bayesian-non-parametric.html)\n  - [Chinese Restaurant Process](https:\/\/m-clark.github.io\/models-by-example\/bayesian-non-parametric.html#chinese-restaurant-process)\n  - [Indian Buffet Process](https:\/\/m-clark.github.io\/models-by-example\/bayesian-non-parametric.html#indian-buffet-process)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-non-parametric.html#source-32)\n- [Bayesian Stochastic Volatility Model](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html#data-setup-31)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html#model-code-8)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html#estimation-30)\n  - [Results](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html#results)\n  - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html#visualization-3)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-stochastic-volatility.html#source-33)\n- [Bayesian Multinomial Models](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html#data-setup-32)\n  - [Model Code](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html#model-code-9)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html#estimation-31)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html#comparison-36)\n  - [Adding Complexity](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html#adding-complexity)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-multinomial-model.html#source-34)\n- [Variational Bayes Regression](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#data-setup-33)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#function-19)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#estimation-32)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#comparison-37)\n  - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#visualization-4)\n  - [Supplemental Example](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#supplemental-example-2)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-variational.html#source-35)\n- [Topic Model](https:\/\/m-clark.github.io\/models-by-example\/bayesian-topic-model.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/bayesian-topic-model.html#data-setup-34)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/bayesian-topic-model.html#function-20)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/bayesian-topic-model.html#estimation-33)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/bayesian-topic-model.html#comparison-38)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/bayesian-topic-model.html#source-36)\n- **Estimation**\n- [Maximum Likelihood](https:\/\/m-clark.github.io\/models-by-example\/maximum-likelihood.html)\n  - [Example](https:\/\/m-clark.github.io\/models-by-example\/maximum-likelihood.html#example-4)\n  - [Linear Model](https:\/\/m-clark.github.io\/models-by-example\/maximum-likelihood.html#linear-model)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/maximum-likelihood.html#source-37)\n- [Penalized Maximum Likelihood](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html)\n  - [Introduction](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#introduction-2)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#data-setup-35)\n    - [Functions](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#functions-8)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#estimation-34)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#comparison-39)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#source-38)\n  - [L1 (lasso) regularization](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#l1-lasso-regularization)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#data-setup-36)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#function-21)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#estimation-35)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#comparison-40)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#source-39)\n  - [L2 (ridge) regularization](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#l2-ridge-regularization)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#data-setup-37)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#function-22)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#estimation-36)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#comparison-41)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/penalized-maximum-likelihood.html#source-40)\n- [Newton and IRLS](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#data-setup-38)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#functions-9)\n    - [Newton’s Method](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#newtons-method)\n    - [IRLS](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#irls)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#estimation-37)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#comparison-42)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/newton-irls.html#source-41)\n- [Nelder-Mead](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#functions-10)\n    - [First Version](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#first-version)\n    - [Example](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#example-5)\n    - [A Regression Model](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#a-regression-model)\n  - [Second Version](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#second-version)\n    - [Example](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#example-6)\n    - [A Regression Model](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#a-regression-model-1)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/nelder-mead.html#source-42)\n- [Expectation-Maximization](https:\/\/m-clark.github.io\/models-by-example\/em.html)\n  - [Mixture Model](https:\/\/m-clark.github.io\/models-by-example\/em.html#mixture-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/em.html#data-setup-40)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/em.html#function-23)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/em.html#estimation-38)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/em.html#comparison-45)\n    - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/em.html#visualization-5)\n    - [Supplemental Example](https:\/\/m-clark.github.io\/models-by-example\/em.html#supplemental-example-3)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/em.html#source-43)\n  - [Multivariate Mixture Model](https:\/\/m-clark.github.io\/models-by-example\/em.html#multivariate-mixture-model)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/em.html#function-24)\n    - [Example 1: Old Faithful](https:\/\/m-clark.github.io\/models-by-example\/em.html#example-1-old-faithful)\n    - [Example 2: Iris](https:\/\/m-clark.github.io\/models-by-example\/em.html#example-2-iris)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/em.html#source-44)\n  - [Probit Model](https:\/\/m-clark.github.io\/models-by-example\/em.html#probit-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/em.html#data-setup-43)\n    - [Probit via Maximum Likelihood](https:\/\/m-clark.github.io\/models-by-example\/em.html#probit-via-maximum-likelihood)\n    - [EM for Latent Variable Approach to Probit](https:\/\/m-clark.github.io\/models-by-example\/em.html#em-for-latent-variable-approach-to-probit)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/em.html#source-45)\n  - [PCA](https:\/\/m-clark.github.io\/models-by-example\/em.html#pca)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/em.html#data-setup-44)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/em.html#function-27)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/em.html#estimation-43)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/em.html#comparison-50)\n    - [Visualize](https:\/\/m-clark.github.io\/models-by-example\/em.html#visualize-4)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/em.html#source-46)\n  - [Probabilistic PCA](https:\/\/m-clark.github.io\/models-by-example\/em.html#probabilistic-pca)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/em.html#data-setup-45)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/em.html#function-28)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/em.html#estimation-44)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/em.html#comparison-51)\n    - [Visualize](https:\/\/m-clark.github.io\/models-by-example\/em.html#visualize-5)\n    - [Missing Data Example](https:\/\/m-clark.github.io\/models-by-example\/em.html#missing-data-example)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/em.html#source-47)\n  - [State Space Model](https:\/\/m-clark.github.io\/models-by-example\/em.html#state-space-model)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/em.html#data-setup-47)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/em.html#function-30)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/em.html#estimation-46)\n    - [Visualization](https:\/\/m-clark.github.io\/models-by-example\/em.html#visualization-6)\n    - [Source](https:\/\/m-clark.github.io\/models-by-example\/em.html#source-48)\n- [Gradient Descent](https:\/\/m-clark.github.io\/models-by-example\/gradient-descent.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/gradient-descent.html#data-setup-48)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/gradient-descent.html#function-31)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/gradient-descent.html#estimation-47)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/gradient-descent.html#comparison-53)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/gradient-descent.html#source-49)\n- [Stochastic Gradient Descent](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#data-setup-49)\n  - [Function](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#function-32)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#estimation-48)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#comparison-54)\n  - [Visualize Estimates](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#visualize-estimates)\n  - [Data Set Shift](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#data-set-shift)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#data-setup-50)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#estimation-49)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#comparison-55)\n    - [Visualize Estimates](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#visualize-estimates-1)\n  - [SGD Variants](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#sgd-variants)\n    - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#data-setup-51)\n    - [Function](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#function-33)\n    - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#estimation-50)\n    - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#comparison-56)\n    - [Visualize Estimates](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#visualize-estimates-2)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/stochastic-gradient-descent.html#source-50)\n- [Metropolis Hastings](https:\/\/m-clark.github.io\/models-by-example\/metropolis-hastings.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/metropolis-hastings.html#data-setup-52)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/metropolis-hastings.html#functions-11)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/metropolis-hastings.html#estimation-51)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/metropolis-hastings.html#comparison-57)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/metropolis-hastings.html#source-51)\n- [Hamiltonian Monte Carlo](https:\/\/m-clark.github.io\/models-by-example\/hamiltonian-monte-carlo.html)\n  - [Data Setup](https:\/\/m-clark.github.io\/models-by-example\/hamiltonian-monte-carlo.html#data-setup-53)\n  - [Functions](https:\/\/m-clark.github.io\/models-by-example\/hamiltonian-monte-carlo.html#functions-12)\n  - [Estimation](https:\/\/m-clark.github.io\/models-by-example\/hamiltonian-monte-carlo.html#estimation-52)\n  - [Comparison](https:\/\/m-clark.github.io\/models-by-example\/hamiltonian-monte-carlo.html#comparison-58)\n  - [Source](https:\/\/m-clark.github.io\/models-by-example\/hamiltonian-monte-carlo.html#source-52)\n- **Other**\n- [Supplemental](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html)\n  - [Other Languages](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#other-languages)\n  - [Python Demos](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-demos)\n    - [Linear Regression](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-linreg)\n    - [Logistic Regression](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-logreg)\n    - [Quantile Regression](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-qreg)\n    - [Stochastic Gradient Descent](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-sgd)\n    - [Nelder-Mead](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-nelder)\n    - [HMM](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#python-hmm)\n  - [Julia Demos](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#julia-demos)\n    - [Mixed Models](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#mixed-models-1)\n  - [Matlab Demos](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#matlab-demos)\n    - [Mixed Models](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#mixed-models-2)\n    - [Gaussian Processes](https:\/\/m-clark.github.io\/models-by-example\/supplemental.html#matlab-gp)\n- [![MC logo](https:\/\/m-clark.github.io\/models-by-example\/img\/mc_logo.png)](https:\/\/m-clark.github.io\/)\n\n# [Model Estimation by Example](https:\/\/m-clark.github.io\/models-by-example\/)\n# Zero-Inflated Model\nLog likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples and comparison to pscl package output. Also includes approach based on Hilbe GLM text.\n\nZero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory.\n\n## Poisson\n### Data Setup\nGet the data.\n```\nlibrary(tidyverse)\n\nfish = haven::read_dta(\"http:\/\/www.stata-press.com\/data\/r11\/fish.dta\")\n```\n\n### Function\nThe log likelihood function.\n```\nzip_ll <- function(y, X, par) {\n  # arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois'\n  \n  # Extract parameters\n  logitpars = par[grep('logit', names(par))]   \n  poispars  = par[grep('pois', names(par))]     \n  \n  # Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example\n  Xlogit  = X\n  LPlogit = Xlogit %*% logitpars\n  logi0   = plogis(LPlogit)  # alternative 1\/(1+exp(-LPlogit))\n    \n  # Poisson part\n  Xpois  = X\n  mupois = exp(Xpois %*% poispars)\n  \n  # LLs\n  logliklogit = log( logi0 + exp(log(1 - logi0) - mupois) )\n  loglikpois  = log(1 - logi0) + dpois(y, lambda = mupois, log = TRUE)\n  \n  # Hilbe formulation\n  # logliklogit = log(logi0 + (1 - logi0)*exp(- mupois) )\n  # loglikpois = log(1-logi0) -mupois + log(mupois)*y     #not necessary: - log(gamma(y+1))\n    \n  y0 = y == 0  # 0 values\n  yc = y > 0   # Count part\n\n  loglik = sum(logliklogit[y0]) + sum(loglikpois[yc])\n  -loglik\n}\n```\n\n### Estimation\nGet starting values or simply do zeros.\n```\n# for zip: need 'logit', 'pois'\ninitial_model = glm(\n  count ~ persons + livebait,\n  data = fish,\n  x = TRUE,\n  y = TRUE,\n  \"poisson\"\n)\n\n# starts = c(logit = coef(initial_model), pois = coef(initial_model))  \nstarts = c(rep(0, 3), rep(0, 3))\n\nnames(starts) = c(paste0('pois.', names(coef(initial_model))),\n                  paste0('logit.', names(coef(initial_model))))\n```\nEstimate with optim.\n```\nfit = optim(\n  par = starts ,\n  fn  = zip_ll,\n  X   = initial_model$x,\n  y   = initial_model$y,\n  method  = \"BFGS\",\n  control = list(maxit = 5000, reltol = 1e-12),\n  hessian = TRUE\n)\n\n# fit\n```\n\n### Comparison\nExtract for clean display.\n```\nB  = fit$par\nse = sqrt(diag(solve((fit$hessian))))\nZ  = B\/se\np  = pnorm(abs(Z), lower = FALSE)*2\n```\nResults from pscl.\n```\nlibrary(pscl)\n\nfit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = \"poisson\")\n```\nCompare.\n\n|   | coef | B | se | Z | p |\n|---|---|---|---|---|---|\n| pscl | X.Intercept. | \\-2.006 | 0\\.324 | \\-6.196 | 0\\.000 |\n| pscl | persons | 0\\.747 | 0\\.043 | 17\\.516 | 0\\.000 |\n| pscl | livebait | 1\\.809 | 0\\.292 | 6\\.195 | 0\\.000 |\n| pscl | X.Intercept..1 | 0\\.303 | 0\\.674 | 0\\.449 | 0\\.654 |\n| pscl | persons.1 | \\-0.069 | 0\\.129 | \\-0.537 | 0\\.591 |\n| pscl | livebait.1 | \\-0.031 | 0\\.558 | \\-0.056 | 0\\.956 |\n| zip\\_poisson\\_ll | pois.(Intercept) | \\-2.006 | 0\\.324 | \\-6.196 | 0\\.000 |\n| zip\\_poisson\\_ll | pois.persons | 0\\.747 | 0\\.043 | 17\\.516 | 0\\.000 |\n| zip\\_poisson\\_ll | pois.livebait | 1\\.809 | 0\\.292 | 6\\.195 | 0\\.000 |\n| zip\\_poisson\\_ll | logit.(Intercept) | 0\\.303 | 0\\.674 | 0\\.449 | 0\\.654 |\n| zip\\_poisson\\_ll | logit.persons | \\-0.069 | 0\\.129 | \\-0.537 | 0\\.591 |\n| zip\\_poisson\\_ll | logit.livebait | \\-0.031 | 0\\.558 | \\-0.056 | 0\\.956 |\n\n## Negative Binomial\n### Function\n```\nzinb_ll <- function(y, X, par) {\n  # arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta'\n  \n  # Extract parameters\n  logitpars  = par[grep('logit', names(par))]\n  negbinpars = par[grep('negbin', names(par))]\n  theta = exp(par[grep('theta', names(par))])\n  \n  # Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example\n  Xlogit  = X\n  LPlogit = Xlogit %*% logitpars\n  logi0   = plogis(LPlogit) \n  \n  # Negbin part\n  Xnegbin = X\n  munb = exp(Xnegbin %*% negbinpars)\n  \n  # LLs\n  logliklogit  = \n    log( \n      logi0 + \n        exp(\n          log(1 - logi0) + \n            suppressWarnings(dnbinom(0, size = theta, mu   = munb, log  = TRUE))\n        )\n    )\n  \n  logliknegbin = log(1 - logi0) + \n    suppressWarnings(dnbinom(y, \n                             size = theta, \n                             mu   = munb, \n                             log  = TRUE))\n  \n  # Hilbe formulation\n  # theta part \n  # alpha = 1\/theta  \n  # m = 1\/alpha\n  # p = 1\/(1 + alpha*munb)\n  \n  # logliklogit = log( logi0 + (1 - logi0)*(p^m) )\n  # logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + m*log(p) + y*log(1-p)   # gamma(y+1) not needed\n  \n  y0 = y == 0   # 0 values\n  yc = y > 0    # Count part\n  \n  loglik = sum(logliklogit[y0]) + sum(logliknegbin[yc])\n  -loglik\n}\n```\n\n### Estimation\nGet starting values or simply do zeros.\n```\n# for zinb: 'logit', 'negbin', 'theta'\ninitial_model = model.matrix(count ~ persons + livebait, data = fish) # to get X matrix\n\nstartlogi  = glm(count == 0 ~ persons + livebait, data = fish, family = \"binomial\")\nstartcount = glm(count ~ persons + livebait, data = fish, family = \"poisson\")\n\nstarts = c(\n  negbin = coef(startcount),\n  logit = coef(startlogi),\n  theta = 1\n)  \n# starts = c(negbin = rep(0, 3),\n#            logit = rep(0, 3),\n#            theta = log(1))\n```\nEstimate with optim.\n```\nfit_nb = optim(\n  par = starts ,\n  fn  = zinb_ll,\n  X   = initial_model,\n  y   = fish$count,\n  method  = \"BFGS\",\n  control = list(maxit = 5000, reltol = 1e-12),\n  hessian = TRUE\n)\n# fit_nb\n```\n\n### Comparison\nExtract for clean display.\n```\nB  = fit_nb$par\nse = sqrt(diag(solve((fit_nb$hessian))))\nZ  = B\/se\np  = pnorm(abs(Z), lower = FALSE)*2\n```\nResults from pscl.\n```\n# pscl results\nlibrary(pscl)\n\nfit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = \"negbin\")\n```\nCompare.\n\n|   | coef | B | se | Z | p |\n|---|---|---|---|---|---|\n| pscl | X.Intercept. | \\-2.803 | 0\\.558 | \\-5.026 | 0\\.000 |\n| pscl | persons | 0\\.849 | 0\\.124 | 6\\.833 | 0\\.000 |\n| pscl | livebait | 1\\.791 | 0\\.511 | 3\\.504 | 0\\.000 |\n| pscl | Log.theta. | \\-0.969 | 0\\.302 | \\-3.206 | 0\\.001 |\n| pscl | X.Intercept..1 | \\-4.276 | 4\\.278 | \\-1.000 | 0\\.318 |\n| pscl | persons.1 | 0\\.560 | 0\\.517 | 1\\.084 | 0\\.279 |\n| pscl | livebait.1 | 1\\.168 | 3\\.661 | 0\\.319 | 0\\.750 |\n| zinb\\_ll | negbin.(Intercept) | \\-2.803 | 0\\.558 | \\-5.026 | 0\\.000 |\n| zinb\\_ll | negbin.persons | 0\\.849 | 0\\.124 | 6\\.834 | 0\\.000 |\n| zinb\\_ll | negbin.livebait | 1\\.791 | 0\\.511 | 3\\.504 | 0\\.000 |\n| zinb\\_ll | logit.(Intercept) | \\-4.276 | 4\\.278 | \\-1.000 | 0\\.318 |\n| zinb\\_ll | logit.persons | 0\\.560 | 0\\.517 | 1\\.084 | 0\\.279 |\n| zinb\\_ll | logit.livebait | 1\\.168 | 3\\.661 | 0\\.319 | 0\\.750 |\n| zinb\\_ll | theta | \\-0.969 | 0\\.302 | \\-3.206 | 0\\.001 |\n\n## Supplemental Example\nThis supplemental example uses the bioChemists data. It contains a sample of 915 biochemistry graduate students with the following information:\n\n- **art**: count of articles produced during last 3 years of Ph.D.\n- **fem**: factor indicating gender of student, with levels Men and Women\n- **mar**: factor indicating marital status of student, with levels Single and Married\n- **kid5**: number of children aged 5 or younger\n- **phd**: prestige of Ph.D. department\n- **ment**: count of articles produced by Ph.D. mentor during last 3 years\n```\ndata(\"bioChemists\", package = \"pscl\")\n\ninitial_model   = model.matrix(art ~ fem + mar + kid5 + phd + ment, data = bioChemists) # to get X matrix\nstartlogi  = glm(art==0 ~ fem + mar + kid5 + phd + ment, data = bioChemists, family = \"binomial\")\nstartcount = glm(art ~ fem + mar + kid5 + phd + ment, data = bioChemists, family = \"quasipoisson\")\n\nstarts = c(\n  negbin = coef(startcount),\n  logit  = coef(startlogi),\n  theta  = summary(startcount)$dispersion\n)  \n\n# starts = c(negbin = rep(0, 6),\n#            logit = rep(0, 6),\n#            theta = 1)\n\n\nfit_nb_pub = optim(\n  par = starts ,\n  fn  = zinb_ll,\n  X   = initial_model,\n  y   = bioChemists$art,\n  method  = \"BFGS\",\n  control = list(maxit = 5000, reltol = 1e-12),\n  hessian = TRUE\n)\n# fit_nb_pub\n\n\nB  = fit_nb_pub$par\nse = sqrt(diag(solve((fit_nb_pub$hessian))))\nZ  = B\/se\np  = pnorm(abs(Z), lower = FALSE)*2\n\n\nlibrary(pscl)\nfit_pscl = zeroinfl(art ~ . | ., data = bioChemists, dist = \"negbin\")\n\nsummary(fit_pscl)$coefficients\n```\n```\n$count\n                 Estimate  Std. Error     z value     Pr(>|z|)\n(Intercept)  0.4167465901 0.143596450  2.90220678 3.705439e-03\nfemWomen    -0.1955076374 0.075592558 -2.58633447 9.700275e-03\nmarMarried   0.0975826042 0.084451953  1.15548073 2.478936e-01\nkid5        -0.1517320709 0.054206071 -2.79917119 5.123397e-03\nphd         -0.0006997593 0.036269674 -0.01929323 9.846072e-01\nment         0.0247861500 0.003492672  7.09661548 1.278491e-12\nLog(theta)   0.9763577454 0.135469554  7.20721163 5.710921e-13\n\n$zero\n               Estimate Std. Error    z value    Pr(>|z|)\n(Intercept) -0.19160645  1.3227962 -0.1448496 0.884829645\nfemWomen     0.63587048  0.8488959  0.7490559 0.453823498\nmarMarried  -1.49943716  0.9386562 -1.5974296 0.110169987\nkid5         0.62840922  0.4427746  1.4192531 0.155825245\nphd         -0.03773288  0.3080059 -0.1225070 0.902497523\nment        -0.88227364  0.3162186 -2.7900755 0.005269575\n```\n```\nround(data.frame(B, se, Z, p), 4)\n```\n```\n                         B     se       Z      p\nnegbin.(Intercept)  0.4167 0.1436  2.9021 0.0037\nnegbin.femWomen    -0.1955 0.0756 -2.5863 0.0097\nnegbin.marMarried   0.0976 0.0845  1.1555 0.2479\nnegbin.kid5        -0.1517 0.0542 -2.7992 0.0051\nnegbin.phd         -0.0007 0.0363 -0.0195 0.9845\nnegbin.ment         0.0248 0.0035  7.0967 0.0000\nlogit.(Intercept)  -0.1916 1.3229 -0.1449 0.8848\nlogit.femWomen      0.6359 0.8490  0.7491 0.4538\nlogit.marMarried   -1.4995 0.9387 -1.5974 0.1102\nlogit.kid5          0.6284 0.4428  1.4192 0.1558\nlogit.phd          -0.0377 0.3080 -0.1225 0.9025\nlogit.ment         -0.8823 0.3162 -2.7901 0.0053\ntheta               0.9763 0.1355  7.2071 0.0000\n```\n\n## Source\nOriginal code for zip\\_ll found at <https:\/\/github.com\/m-clark\/Miscellaneous-R-Code\/blob\/master\/ModelFitting\/poiszeroinfl.R>\n\nOriginal code for zinb\\_ll found at <https:\/\/github.com\/m-clark\/Miscellaneous-R-Code\/blob\/master\/ModelFitting\/NBzeroinfl.R>","attrs_readable_markdown":"Log likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples and comparison to pscl package output. Also includes approach based on Hilbe GLM text.\n\nZero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory.\n\n## Poisson\n### Data Setup\nGet the data.\n```\nlibrary(tidyverse)\n\nfish = haven::read_dta(\"http:\/\/www.stata-press.com\/data\/r11\/fish.dta\")\n```\n\n### Function\nThe log likelihood function.\n```\nzip_ll <- function(y, X, par) {\n  # arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois'\n  \n  # Extract parameters\n  logitpars = par[grep('logit', names(par))]   \n  poispars  = par[grep('pois', names(par))]     \n  \n  # Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example\n  Xlogit  = X\n  LPlogit = Xlogit %*% logitpars\n  logi0   = plogis(LPlogit)  # alternative 1\/(1+exp(-LPlogit))\n    \n  # Poisson part\n  Xpois  = X\n  mupois = exp(Xpois %*% poispars)\n  \n  # LLs\n  logliklogit = log( logi0 + exp(log(1 - logi0) - mupois) )\n  loglikpois  = log(1 - logi0) + dpois(y, lambda = mupois, log = TRUE)\n  \n  # Hilbe formulation\n  # logliklogit = log(logi0 + (1 - logi0)*exp(- mupois) )\n  # loglikpois = log(1-logi0) -mupois + log(mupois)*y     #not necessary: - log(gamma(y+1))\n    \n  y0 = y == 0  # 0 values\n  yc = y > 0   # Count part\n\n  loglik = sum(logliklogit[y0]) + sum(loglikpois[yc])\n  -loglik\n}\n```\n\n### Estimation\nGet starting values or simply do zeros.\n```\n# for zip: need 'logit', 'pois'\ninitial_model = glm(\n  count ~ persons + livebait,\n  data = fish,\n  x = TRUE,\n  y = TRUE,\n  \"poisson\"\n)\n\n# starts = c(logit = coef(initial_model), pois = coef(initial_model))  \nstarts = c(rep(0, 3), rep(0, 3))\n\nnames(starts) = c(paste0('pois.', names(coef(initial_model))),\n                  paste0('logit.', names(coef(initial_model))))\n```\nEstimate with optim.\n```\nfit = optim(\n  par = starts ,\n  fn  = zip_ll,\n  X   = initial_model$x,\n  y   = initial_model$y,\n  method  = \"BFGS\",\n  control = list(maxit = 5000, reltol = 1e-12),\n  hessian = TRUE\n)\n\n# fit\n```\n\n### Comparison\nExtract for clean display.\n```\nB  = fit$par\nse = sqrt(diag(solve((fit$hessian))))\nZ  = B\/se\np  = pnorm(abs(Z), lower = FALSE)*2\n```\nResults from pscl.\n```\nlibrary(pscl)\n\nfit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = \"poisson\")\n```\nCompare.\n\n|   | coef | B | se | Z | p |\n|---|---|---|---|---|---|\n| pscl | X.Intercept. | \\-2.006 | 0\\.324 | \\-6.196 | 0\\.000 |\n| pscl | persons | 0\\.747 | 0\\.043 | 17\\.516 | 0\\.000 |\n| pscl | livebait | 1\\.809 | 0\\.292 | 6\\.195 | 0\\.000 |\n| pscl | X.Intercept..1 | 0\\.303 | 0\\.674 | 0\\.449 | 0\\.654 |\n| pscl | persons.1 | \\-0.069 | 0\\.129 | \\-0.537 | 0\\.591 |\n| pscl | livebait.1 | \\-0.031 | 0\\.558 | \\-0.056 | 0\\.956 |\n| zip\\_poisson\\_ll | pois.(Intercept) | \\-2.006 | 0\\.324 | \\-6.196 | 0\\.000 |\n| zip\\_poisson\\_ll | pois.persons | 0\\.747 | 0\\.043 | 17\\.516 | 0\\.000 |\n| zip\\_poisson\\_ll | pois.livebait | 1\\.809 | 0\\.292 | 6\\.195 | 0\\.000 |\n| zip\\_poisson\\_ll | logit.(Intercept) | 0\\.303 | 0\\.674 | 0\\.449 | 0\\.654 |\n| zip\\_poisson\\_ll | logit.persons | \\-0.069 | 0\\.129 | \\-0.537 | 0\\.591 |\n| zip\\_poisson\\_ll | logit.livebait | \\-0.031 | 0\\.558 | \\-0.056 | 0\\.956 |\n\n## Negative Binomial\n### Function\n```\nzinb_ll <- function(y, X, par) {\n  # arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta'\n  \n  # Extract parameters\n  logitpars  = par[grep('logit', names(par))]\n  negbinpars = par[grep('negbin', names(par))]\n  theta = exp(par[grep('theta', names(par))])\n  \n  # Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example\n  Xlogit  = X\n  LPlogit = Xlogit %*% logitpars\n  logi0   = plogis(LPlogit) \n  \n  # Negbin part\n  Xnegbin = X\n  munb = exp(Xnegbin %*% negbinpars)\n  \n  # LLs\n  logliklogit  = \n    log( \n      logi0 + \n        exp(\n          log(1 - logi0) + \n            suppressWarnings(dnbinom(0, size = theta, mu   = munb, log  = TRUE))\n        )\n    )\n  \n  logliknegbin = log(1 - logi0) + \n    suppressWarnings(dnbinom(y, \n                             size = theta, \n                             mu   = munb, \n                             log  = TRUE))\n  \n  # Hilbe formulation\n  # theta part \n  # alpha = 1\/theta  \n  # m = 1\/alpha\n  # p = 1\/(1 + alpha*munb)\n  \n  # logliklogit = log( logi0 + (1 - logi0)*(p^m) )\n  # logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + m*log(p) + y*log(1-p)   # gamma(y+1) not needed\n  \n  y0 = y == 0   # 0 values\n  yc = y > 0    # Count part\n  \n  loglik = sum(logliklogit[y0]) + sum(logliknegbin[yc])\n  -loglik\n}\n```\n\n### Estimation\nGet starting values or simply do zeros.\n```\n# for zinb: 'logit', 'negbin', 'theta'\ninitial_model = model.matrix(count ~ persons + livebait, data = fish) # to get X matrix\n\nstartlogi  = glm(count == 0 ~ persons + livebait, data = fish, family = \"binomial\")\nstartcount = glm(count ~ persons + livebait, data = fish, family = \"poisson\")\n\nstarts = c(\n  negbin = coef(startcount),\n  logit = coef(startlogi),\n  theta = 1\n)  \n# starts = c(negbin = rep(0, 3),\n#            logit = rep(0, 3),\n#            theta = log(1))\n```\nEstimate with optim.\n```\nfit_nb = optim(\n  par = starts ,\n  fn  = zinb_ll,\n  X   = initial_model,\n  y   = fish$count,\n  method  = \"BFGS\",\n  control = list(maxit = 5000, reltol = 1e-12),\n  hessian = TRUE\n)\n# fit_nb\n```\n\n### Comparison\nExtract for clean display.\n```\nB  = fit_nb$par\nse = sqrt(diag(solve((fit_nb$hessian))))\nZ  = B\/se\np  = pnorm(abs(Z), lower = FALSE)*2\n```\nResults from pscl.\n```\n# pscl results\nlibrary(pscl)\n\nfit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = \"negbin\")\n```\nCompare.\n\n|   | coef | B | se | Z | p |\n|---|---|---|---|---|---|\n| pscl | X.Intercept. | \\-2.803 | 0\\.558 | \\-5.026 | 0\\.000 |\n| pscl | persons | 0\\.849 | 0\\.124 | 6\\.833 | 0\\.000 |\n| pscl | livebait | 1\\.791 | 0\\.511 | 3\\.504 | 0\\.000 |\n| pscl | Log.theta. | \\-0.969 | 0\\.302 | \\-3.206 | 0\\.001 |\n| pscl | X.Intercept..1 | \\-4.276 | 4\\.278 | \\-1.000 | 0\\.318 |\n| pscl | persons.1 | 0\\.560 | 0\\.517 | 1\\.084 | 0\\.279 |\n| pscl | livebait.1 | 1\\.168 | 3\\.661 | 0\\.319 | 0\\.750 |\n| zinb\\_ll | negbin.(Intercept) | \\-2.803 | 0\\.558 | \\-5.026 | 0\\.000 |\n| zinb\\_ll | negbin.persons | 0\\.849 | 0\\.124 | 6\\.834 | 0\\.000 |\n| zinb\\_ll | negbin.livebait | 1\\.791 | 0\\.511 | 3\\.504 | 0\\.000 |\n| zinb\\_ll | logit.(Intercept) | \\-4.276 | 4\\.278 | \\-1.000 | 0\\.318 |\n| zinb\\_ll | logit.persons | 0\\.560 | 0\\.517 | 1\\.084 | 0\\.279 |\n| zinb\\_ll | logit.livebait | 1\\.168 | 3\\.661 | 0\\.319 | 0\\.750 |\n| zinb\\_ll | theta | \\-0.969 | 0\\.302 | \\-3.206 | 0\\.001 |","meta_canonical":null}

3. Robots.txt Check

Query:

Response:

4. Spam/Ban Check

Query:

Response:

5. Seen Status Check

ℹ️ Skipped - page is already crawled

📄

INDEXABLE

✅

CRAWLED

18 days ago

🤖

ROBOTS ALLOWED

Page Info Filters

Filter	Status	Condition	Details
HTTP status	PASS	`download_http_code = 200`	HTTP 200
Age cutoff	PASS	`download_stamp > now() - 6 MONTH`	0.6 months ago
History drop	PASS	`isNull(history_drop_reason)`	No drop reason
Spam/ban	PASS	`fh_dont_index != 1 AND ml_spam_score = 0`	ml_spam_score=0
Canonical	PASS	`meta_canonical IS NULL OR = '' OR = src_unparsed`	Not set

Page Details

Property	Value
URL	https://m-clark.github.io/models-by-example/zi.html
Last Crawled	2026-03-19 19:35:13 (18 days ago)
First Indexed	2020-11-11 20:16:36 (5 years ago)
HTTP Status Code	200
Meta Title	Zero-Inflated Model \| Model Estimation by Example
Meta Description	This document provides ‘by-hand’ demonstrations of various models and algorithms. The goal is to take away some of the mystery by providing clean code examples that are easy to run and compare with other tools.
Meta Canonical	null
Boilerpipe Text	Log likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples and comparison to pscl package output. Also includes approach based on Hilbe GLM text. Zero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory. Poisson Data Setup Get the data. library (tidyverse) fish = haven :: read_dta ( "http://www.stata-press.com/data/r11/fish.dta" ) Function The log likelihood function. zip_ll <- function (y, X, par) { # arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois' # Extract parameters logitpars = par[ grep ( 'logit' , names (par))] poispars = par[ grep ( 'pois' , names (par))] # Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example Xlogit = X LPlogit = Xlogit %% logitpars logi0 = plogis (LPlogit) # alternative 1/(1+exp(-LPlogit)) # Poisson part Xpois = X mupois = exp (Xpois %% poispars) # LLs logliklogit = log ( logi0 + exp ( log ( 1 - logi0) - mupois) ) loglikpois = log ( 1 - logi0) + dpois (y, lambda = mupois, log = TRUE ) # Hilbe formulation # logliklogit = log(logi0 + (1 - logi0)exp(- mupois) ) # loglikpois = log(1-logi0) -mupois + log(mupois)y #not necessary: - log(gamma(y+1)) y0 = y == 0 # 0 values yc = y > 0 # Count part loglik = sum (logliklogit[y0]) + sum (loglikpois[yc]) - loglik } Estimation Get starting values or simply do zeros. # for zip: need 'logit', 'pois' initial_model = glm ( count ~ persons + livebait, data = fish, x = TRUE , y = TRUE , "poisson" ) # starts = c(logit = coef(initial_model), pois = coef(initial_model)) starts = c ( rep ( 0 , 3 ), rep ( 0 , 3 )) names (starts) = c ( paste0 ( 'pois.' , names ( coef (initial_model))), paste0 ( 'logit.' , names ( coef (initial_model)))) Estimate with optim . fit = optim ( par = starts , fn = zip_ll, X = initial_model $ x, y = initial_model $ y, method = "BFGS" , control = list ( maxit = 5000 , reltol = 1e-12 ), hessian = TRUE ) # fit Comparison Extract for clean display. B = fit $ par se = sqrt ( diag ( solve ((fit $ hessian)))) Z = B / se p = pnorm ( abs (Z), lower = FALSE ) * 2 Results from pscl . library (pscl) fit_pscl = zeroinfl (count ~ persons + livebait, data = fish, dist = "poisson" ) Compare. coef B se Z p pscl X.Intercept. -2.006 0.324 -6.196 0.000 pscl persons 0.747 0.043 17.516 0.000 pscl livebait 1.809 0.292 6.195 0.000 pscl X.Intercept..1 0.303 0.674 0.449 0.654 pscl persons.1 -0.069 0.129 -0.537 0.591 pscl livebait.1 -0.031 0.558 -0.056 0.956 zip_poisson_ll pois.(Intercept) -2.006 0.324 -6.196 0.000 zip_poisson_ll pois.persons 0.747 0.043 17.516 0.000 zip_poisson_ll pois.livebait 1.809 0.292 6.195 0.000 zip_poisson_ll logit.(Intercept) 0.303 0.674 0.449 0.654 zip_poisson_ll logit.persons -0.069 0.129 -0.537 0.591 zip_poisson_ll logit.livebait -0.031 0.558 -0.056 0.956 Negative Binomial Function zinb_ll <- function (y, X, par) { # arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta' # Extract parameters logitpars = par[ grep ( 'logit' , names (par))] negbinpars = par[ grep ( 'negbin' , names (par))] theta = exp (par[ grep ( 'theta' , names (par))]) # Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example Xlogit = X LPlogit = Xlogit %% logitpars logi0 = plogis (LPlogit) # Negbin part Xnegbin = X munb = exp (Xnegbin %% negbinpars) # LLs logliklogit = log ( logi0 + exp ( log ( 1 - logi0) + suppressWarnings ( dnbinom ( 0 , size = theta, mu = munb, log = TRUE )) ) ) logliknegbin = log ( 1 - logi0) + suppressWarnings ( dnbinom (y, size = theta, mu = munb, log = TRUE )) # Hilbe formulation # theta part # alpha = 1/theta # m = 1/alpha # p = 1/(1 + alphamunb) # logliklogit = log( logi0 + (1 - logi0)(p^m) ) # logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + mlog(p) + ylog(1-p) # gamma(y+1) not needed y0 = y == 0 # 0 values yc = y > 0 # Count part loglik = sum (logliklogit[y0]) + sum (logliknegbin[yc]) - loglik } Estimation Get starting values or simply do zeros. # for zinb: 'logit', 'negbin', 'theta' initial_model = model.matrix (count ~ persons + livebait, data = fish) # to get X matrix startlogi = glm (count == 0 ~ persons + livebait, data = fish, family = "binomial" ) startcount = glm (count ~ persons + livebait, data = fish, family = "poisson" ) starts = c ( negbin = coef (startcount), logit = coef (startlogi), theta = 1 ) # starts = c(negbin = rep(0, 3), # logit = rep(0, 3), # theta = log(1)) Estimate with optim . fit_nb = optim ( par = starts , fn = zinb_ll, X = initial_model, y = fish $ count, method = "BFGS" , control = list ( maxit = 5000 , reltol = 1e-12 ), hessian = TRUE ) # fit_nb Comparison Extract for clean display. B = fit_nb $ par se = sqrt ( diag ( solve ((fit_nb $ hessian)))) Z = B / se p = pnorm ( abs (Z), lower = FALSE ) * 2 Results from pscl . # pscl results library (pscl) fit_pscl = zeroinfl (count ~ persons + livebait, data = fish, dist = "negbin" ) Compare. coef B se Z p pscl X.Intercept. -2.803 0.558 -5.026 0.000 pscl persons 0.849 0.124 6.833 0.000 pscl livebait 1.791 0.511 3.504 0.000 pscl Log.theta. -0.969 0.302 -3.206 0.001 pscl X.Intercept..1 -4.276 4.278 -1.000 0.318 pscl persons.1 0.560 0.517 1.084 0.279 pscl livebait.1 1.168 3.661 0.319 0.750 zinb_ll negbin.(Intercept) -2.803 0.558 -5.026 0.000 zinb_ll negbin.persons 0.849 0.124 6.834 0.000 zinb_ll negbin.livebait 1.791 0.511 3.504 0.000 zinb_ll logit.(Intercept) -4.276 4.278 -1.000 0.318 zinb_ll logit.persons 0.560 0.517 1.084 0.279 zinb_ll logit.livebait 1.168 3.661 0.319 0.750 zinb_ll theta -0.969 0.302 -3.206 0.001
Markdown	- [Model Estimation by Example](https://m-clark.github.io/models-by-example/) - [Introduction](https://m-clark.github.io/models-by-example/introduction.html) - Models - [Linear Regression](https://m-clark.github.io/models-by-example/linear-regression.html) - [Data Setup](https://m-clark.github.io/models-by-example/linear-regression.html#data-setup) - [Functions](https://m-clark.github.io/models-by-example/linear-regression.html#functions) - [Estimation](https://m-clark.github.io/models-by-example/linear-regression.html#estimation) - [Comparison](https://m-clark.github.io/models-by-example/linear-regression.html#comparison) - [Python](https://m-clark.github.io/models-by-example/linear-regression.html#python) - [Source](https://m-clark.github.io/models-by-example/linear-regression.html#source) - [Logistic Regression](https://m-clark.github.io/models-by-example/logistic-regression.html) - [Data Setup](https://m-clark.github.io/models-by-example/logistic-regression.html#data-setup-1) - [Functions](https://m-clark.github.io/models-by-example/logistic-regression.html#functions-1) - [Estimation](https://m-clark.github.io/models-by-example/logistic-regression.html#estimation-1) - [Comparison](https://m-clark.github.io/models-by-example/logistic-regression.html#comparison-1) - [Python](https://m-clark.github.io/models-by-example/logistic-regression.html#python-1) - [Source](https://m-clark.github.io/models-by-example/logistic-regression.html#source-1) - [Mixed Models](https://m-clark.github.io/models-by-example/mixed-models.html) - [One-factor Mixed Model](https://m-clark.github.io/models-by-example/mixed-models.html#one-factor-mixed-model) - [Data Setup](https://m-clark.github.io/models-by-example/mixed-models.html#data-setup-2) - [Function](https://m-clark.github.io/models-by-example/mixed-models.html#function) - [Estimation](https://m-clark.github.io/models-by-example/mixed-models.html#estimation-2) - [Comparison](https://m-clark.github.io/models-by-example/mixed-models.html#comparison-2) - [Source](https://m-clark.github.io/models-by-example/mixed-models.html#source-2) - [Two-factor Mixed Model](https://m-clark.github.io/models-by-example/mixed-models.html#two-factor-mixed-model) - [Data Setup](https://m-clark.github.io/models-by-example/mixed-models.html#data-setup-3) - [Function](https://m-clark.github.io/models-by-example/mixed-models.html#function-1) - [Estimation](https://m-clark.github.io/models-by-example/mixed-models.html#estimation-3) - [Comparison](https://m-clark.github.io/models-by-example/mixed-models.html#comparison-3) - [Source](https://m-clark.github.io/models-by-example/mixed-models.html#source-3) - [Mixed Model via ML](https://m-clark.github.io/models-by-example/mixed-models.html#mixed-model-via-ml) - [Introduction](https://m-clark.github.io/models-by-example/mixed-models.html#introduction-1) - [Data Setup](https://m-clark.github.io/models-by-example/mixed-models.html#data-setup-4) - [Function](https://m-clark.github.io/models-by-example/mixed-models.html#function-2) - [Estimation](https://m-clark.github.io/models-by-example/mixed-models.html#estimation-4) - [Comparison](https://m-clark.github.io/models-by-example/mixed-models.html#comparison-4) - [Link with penalized regression](https://m-clark.github.io/models-by-example/mixed-models.html#link-with-penalized-regression) - [Source](https://m-clark.github.io/models-by-example/mixed-models.html#source-4) - [Probit & Bivariate Probit](https://m-clark.github.io/models-by-example/probit.html) - [Standard Probit](https://m-clark.github.io/models-by-example/probit.html#standard-probit) - [Function](https://m-clark.github.io/models-by-example/probit.html#function-3) - [Examples](https://m-clark.github.io/models-by-example/probit.html#examples) - [Bivariate Probit](https://m-clark.github.io/models-by-example/probit.html#bivariate-probit) - [Example](https://m-clark.github.io/models-by-example/probit.html#example) - [Source](https://m-clark.github.io/models-by-example/probit.html#source-5) - [Heckman Selection](https://m-clark.github.io/models-by-example/heckman-selection.html) - [Data Setup](https://m-clark.github.io/models-by-example/heckman-selection.html#data-setup-5) - [Two step approach](https://m-clark.github.io/models-by-example/heckman-selection.html#two-step-approach) - [Step 1: Probit Model](https://m-clark.github.io/models-by-example/heckman-selection.html#step-1-probit-model) - [Step 2: Estimate via Linear Regression](https://m-clark.github.io/models-by-example/heckman-selection.html#step-2-estimate-via-linear-regression) - [Maximum Likelihood](https://m-clark.github.io/models-by-example/heckman-selection.html#maximum-likelihood) - [Comparison](https://m-clark.github.io/models-by-example/heckman-selection.html#comparison-5) - [Source](https://m-clark.github.io/models-by-example/heckman-selection.html#source-6) - [Marginal Structural Model](https://m-clark.github.io/models-by-example/marginal-structural.html) - [Data Setup](https://m-clark.github.io/models-by-example/marginal-structural.html#data-setup-6) - [Function](https://m-clark.github.io/models-by-example/marginal-structural.html#function-4) - [Estimation](https://m-clark.github.io/models-by-example/marginal-structural.html#estimation-5) - [Comparison](https://m-clark.github.io/models-by-example/marginal-structural.html#comparison-6) - [Source](https://m-clark.github.io/models-by-example/marginal-structural.html#source-7) - [Tobit Regression](https://m-clark.github.io/models-by-example/tobit.html) - [Censoring with an Upper Limit](https://m-clark.github.io/models-by-example/tobit.html#censoring-with-an-upper-limit) - [Data Setup](https://m-clark.github.io/models-by-example/tobit.html#data-setup-7) - [Function](https://m-clark.github.io/models-by-example/tobit.html#function-5) - [Estimation](https://m-clark.github.io/models-by-example/tobit.html#estimation-6) - [Comparison](https://m-clark.github.io/models-by-example/tobit.html#comparison-7) - [Censoring with a Lower Limit](https://m-clark.github.io/models-by-example/tobit.html#censoring-with-a-lower-limit) - [Comparison](https://m-clark.github.io/models-by-example/tobit.html#comparison-8) - [Source](https://m-clark.github.io/models-by-example/tobit.html#source-8) - [Cox Survival](https://m-clark.github.io/models-by-example/cox.html) - [Standard Proportional Hazards](https://m-clark.github.io/models-by-example/cox.html#standard-proportional-hazards) - [Data Setup](https://m-clark.github.io/models-by-example/cox.html#data-setup-8) - [Function](https://m-clark.github.io/models-by-example/cox.html#function-6) - [Estimation](https://m-clark.github.io/models-by-example/cox.html#estimation-7) - [Comparison](https://m-clark.github.io/models-by-example/cox.html#comparison-9) - [Time-varying coefficients](https://m-clark.github.io/models-by-example/cox.html#time-varying-coefficients) - [Data Setup](https://m-clark.github.io/models-by-example/cox.html#data-setup-9) - [Function](https://m-clark.github.io/models-by-example/cox.html#function-7) - [Estimation](https://m-clark.github.io/models-by-example/cox.html#estimation-8) - [Comparison](https://m-clark.github.io/models-by-example/cox.html#comparison-10) - [Stratified Cox Model](https://m-clark.github.io/models-by-example/cox.html#stratified-cox-model) - [Data Setup](https://m-clark.github.io/models-by-example/cox.html#data-setup-10) - [Function](https://m-clark.github.io/models-by-example/cox.html#function-8) - [Estimation](https://m-clark.github.io/models-by-example/cox.html#estimation-9) - [Comparison](https://m-clark.github.io/models-by-example/cox.html#comparison-11) - [Source](https://m-clark.github.io/models-by-example/cox.html#source-9) - [Hurdle Model](https://m-clark.github.io/models-by-example/hurdle.html) - [Poisson](https://m-clark.github.io/models-by-example/hurdle.html#poisson) - [Data Setup](https://m-clark.github.io/models-by-example/hurdle.html#data-setup-11) - [Function](https://m-clark.github.io/models-by-example/hurdle.html#function-9) - [Estimation](https://m-clark.github.io/models-by-example/hurdle.html#estimation-10) - [Comparison](https://m-clark.github.io/models-by-example/hurdle.html#comparison-12) - [Negative Binomial](https://m-clark.github.io/models-by-example/hurdle.html#negative-binomial) - [Function](https://m-clark.github.io/models-by-example/hurdle.html#function-10) - [Estimation](https://m-clark.github.io/models-by-example/hurdle.html#estimation-11) - [Comparison](https://m-clark.github.io/models-by-example/hurdle.html#comparison-13) - [Source](https://m-clark.github.io/models-by-example/hurdle.html#source-10) - [Zero-Inflated Model](https://m-clark.github.io/models-by-example/zi.html) - [Poisson](https://m-clark.github.io/models-by-example/zi.html#poisson-1) - [Data Setup](https://m-clark.github.io/models-by-example/zi.html#data-setup-12) - [Function](https://m-clark.github.io/models-by-example/zi.html#function-11) - [Estimation](https://m-clark.github.io/models-by-example/zi.html#estimation-12) - [Comparison](https://m-clark.github.io/models-by-example/zi.html#comparison-14) - [Negative Binomial](https://m-clark.github.io/models-by-example/zi.html#negative-binomial-1) - [Function](https://m-clark.github.io/models-by-example/zi.html#function-12) - [Estimation](https://m-clark.github.io/models-by-example/zi.html#estimation-13) - [Comparison](https://m-clark.github.io/models-by-example/zi.html#comparison-15) - [Supplemental Example](https://m-clark.github.io/models-by-example/zi.html#supplemental-example) - [Source](https://m-clark.github.io/models-by-example/zi.html#source-11) - [Naive Bayes](https://m-clark.github.io/models-by-example/naive-bayes.html) - [Initialization](https://m-clark.github.io/models-by-example/naive-bayes.html#initialization) - [Comparison](https://m-clark.github.io/models-by-example/naive-bayes.html#comparison-16) - [Source](https://m-clark.github.io/models-by-example/naive-bayes.html#source-12) - [Multinomial Regression](https://m-clark.github.io/models-by-example/multinomial.html) - [Standard (Categorical) Model](https://m-clark.github.io/models-by-example/multinomial.html#standard-categorical-model) - [Data Setup](https://m-clark.github.io/models-by-example/multinomial.html#data-setup-13) - [Function](https://m-clark.github.io/models-by-example/multinomial.html#function-13) - [Comparison](https://m-clark.github.io/models-by-example/multinomial.html#comparison-17) - [Alternative specific and constant variables](https://m-clark.github.io/models-by-example/multinomial.html#alternative-specific-and-constant-variables) - [Comparison](https://m-clark.github.io/models-by-example/multinomial.html#comparison-18) - [Source](https://m-clark.github.io/models-by-example/multinomial.html#source-13) - [Ordinal Regression](https://m-clark.github.io/models-by-example/ordinal.html) - [Data](https://m-clark.github.io/models-by-example/ordinal.html#data) - [Function](https://m-clark.github.io/models-by-example/ordinal.html#function-14) - [Estimation](https://m-clark.github.io/models-by-example/ordinal.html#estimation-14) - [Comparison](https://m-clark.github.io/models-by-example/ordinal.html#comparison-19) - [Source](https://m-clark.github.io/models-by-example/ordinal.html#source-14) - [Markov Model](https://m-clark.github.io/models-by-example/markov.html) - [Data Setup](https://m-clark.github.io/models-by-example/markov.html#data-setup-14) - [Data Functions](https://m-clark.github.io/models-by-example/markov.html#data-functions) - [Two states Demo](https://m-clark.github.io/models-by-example/markov.html#two-states-demo) - [Three states demo](https://m-clark.github.io/models-by-example/markov.html#three-states-demo) - [Function](https://m-clark.github.io/models-by-example/markov.html#function-15) - [Estimation](https://m-clark.github.io/models-by-example/markov.html#estimation-15) - [Comparison](https://m-clark.github.io/models-by-example/markov.html#comparison-20) - [Source](https://m-clark.github.io/models-by-example/markov.html#source-15) - [Hidden Markov Model](https://m-clark.github.io/models-by-example/hmm.html) - [Data Setup](https://m-clark.github.io/models-by-example/hmm.html#data-setup-15) - [Function](https://m-clark.github.io/models-by-example/hmm.html#function-16) - [Estimation](https://m-clark.github.io/models-by-example/hmm.html#estimation-16) - [Supplemental demo](https://m-clark.github.io/models-by-example/hmm.html#supplemental-demo) - [Source](https://m-clark.github.io/models-by-example/hmm.html#source-16) - [Quantile Regression](https://m-clark.github.io/models-by-example/quantile-regression.html) - [Data Setup](https://m-clark.github.io/models-by-example/quantile-regression.html#data-setup-16) - [Function](https://m-clark.github.io/models-by-example/quantile-regression.html#function-17) - [Estimation](https://m-clark.github.io/models-by-example/quantile-regression.html#estimation-17) - [Other quantiles](https://m-clark.github.io/models-by-example/quantile-regression.html#other-quantiles) - [Comparison](https://m-clark.github.io/models-by-example/quantile-regression.html#comparison-21) - [Visualize](https://m-clark.github.io/models-by-example/quantile-regression.html#visualize) - [Python](https://m-clark.github.io/models-by-example/quantile-regression.html#python-2) - [Source](https://m-clark.github.io/models-by-example/quantile-regression.html#source-17) - [Cubic Spline Model](https://m-clark.github.io/models-by-example/cubic-spline.html) - [Data Setup](https://m-clark.github.io/models-by-example/cubic-spline.html#data-setup-17) - [Functions](https://m-clark.github.io/models-by-example/cubic-spline.html#functions-2) - [Example 1](https://m-clark.github.io/models-by-example/cubic-spline.html#example-1) - [Example 2](https://m-clark.github.io/models-by-example/cubic-spline.html#example-2) - [Source](https://m-clark.github.io/models-by-example/cubic-spline.html#source-18) - [Gaussian Processes](https://m-clark.github.io/models-by-example/gaussian-process.html) - [Noise-Free Demonstration](https://m-clark.github.io/models-by-example/gaussian-process.html#noise-free-demonstration) - [Data Setup](https://m-clark.github.io/models-by-example/gaussian-process.html#data-setup-18) - [Functions](https://m-clark.github.io/models-by-example/gaussian-process.html#functions-3) - [Visualize the prior distribution](https://m-clark.github.io/models-by-example/gaussian-process.html#visualize-the-prior-distribution) - [Generate the posterior predictive distribution](https://m-clark.github.io/models-by-example/gaussian-process.html#generate-the-posterior-predictive-distribution) - [Visualize the Posterior Predictive Distribution](https://m-clark.github.io/models-by-example/gaussian-process.html#visualize-the-posterior-predictive-distribution) - [Noisy Demonstration](https://m-clark.github.io/models-by-example/gaussian-process.html#noisy-demonstration) - [Data Setup](https://m-clark.github.io/models-by-example/gaussian-process.html#data-setup-19) - [Functions](https://m-clark.github.io/models-by-example/gaussian-process.html#functions-4) - [Visualize the prior distribution](https://m-clark.github.io/models-by-example/gaussian-process.html#visualize-the-prior-distribution-1) - [Generate the posterior predictive distribution](https://m-clark.github.io/models-by-example/gaussian-process.html#generate-the-posterior-predictive-distribution-1) - [Visualize the Posterior Predictive Distribution](https://m-clark.github.io/models-by-example/gaussian-process.html#visualize-the-posterior-predictive-distribution-1) - [Source](https://m-clark.github.io/models-by-example/gaussian-process.html#source-19) - [Neural Network](https://m-clark.github.io/models-by-example/neural-network.html) - [Example 1](https://m-clark.github.io/models-by-example/neural-network.html#example-1-1) - [Python](https://m-clark.github.io/models-by-example/neural-network.html#python-3) - [R](https://m-clark.github.io/models-by-example/neural-network.html#r) - [Example 2](https://m-clark.github.io/models-by-example/neural-network.html#example-2-1) - [Python](https://m-clark.github.io/models-by-example/neural-network.html#python-4) - [R](https://m-clark.github.io/models-by-example/neural-network.html#r-1) - [Comparison](https://m-clark.github.io/models-by-example/neural-network.html#comparison-22) - [Example 3](https://m-clark.github.io/models-by-example/neural-network.html#example-3) - [Python](https://m-clark.github.io/models-by-example/neural-network.html#python-5) - [R](https://m-clark.github.io/models-by-example/neural-network.html#r-2) - [Comparison](https://m-clark.github.io/models-by-example/neural-network.html#comparison-23) - [Source](https://m-clark.github.io/models-by-example/neural-network.html#source-20) - [Extreme Learning Machine](https://m-clark.github.io/models-by-example/elm.html) - [Data Setup](https://m-clark.github.io/models-by-example/elm.html#data-setup-20) - [Function](https://m-clark.github.io/models-by-example/elm.html#function-18) - [Estimation](https://m-clark.github.io/models-by-example/elm.html#estimation-18) - [Comparison](https://m-clark.github.io/models-by-example/elm.html#comparison-24) - [Supplemental Example](https://m-clark.github.io/models-by-example/elm.html#supplemental-example-1) - [Source](https://m-clark.github.io/models-by-example/elm.html#source-21) - [Reproducing Kernel Hilbert Space Regression](https://m-clark.github.io/models-by-example/rkhs.html) - [Data Setup](https://m-clark.github.io/models-by-example/rkhs.html#data-setup-21) - [Functions](https://m-clark.github.io/models-by-example/rkhs.html#functions-5) - [Function to find the inverse of a matrix](https://m-clark.github.io/models-by-example/rkhs.html#function-to-find-the-inverse-of-a-matrix) - [Reproducing Kernel](https://m-clark.github.io/models-by-example/rkhs.html#reproducing-kernel) - [Gram matrix](https://m-clark.github.io/models-by-example/rkhs.html#gram-matrix) - [Ridge regression](https://m-clark.github.io/models-by-example/rkhs.html#ridge-regression) - [Estimation](https://m-clark.github.io/models-by-example/rkhs.html#estimation-19) - [Comparison](https://m-clark.github.io/models-by-example/rkhs.html#comparison-25) - [Example: Cubic Spline](https://m-clark.github.io/models-by-example/rkhs.html#example-cubic-spline) - [Data Setup](https://m-clark.github.io/models-by-example/rkhs.html#data-setup-22) - [Functions](https://m-clark.github.io/models-by-example/rkhs.html#functions-6) - [Estimation](https://m-clark.github.io/models-by-example/rkhs.html#estimation-20) - [Comparison](https://m-clark.github.io/models-by-example/rkhs.html#comparison-26) - [Source](https://m-clark.github.io/models-by-example/rkhs.html#source-22) - [Current Supplemental Code](https://m-clark.github.io/models-by-example/rkhs.html#current-supplemental-code) - [Original Supplemental Code](https://m-clark.github.io/models-by-example/rkhs.html#original-supplemental-code) - [Confirmatory Factor Analysis](https://m-clark.github.io/models-by-example/cfa.html) - [Data Setup](https://m-clark.github.io/models-by-example/cfa.html#data-setup-23) - [Functions](https://m-clark.github.io/models-by-example/cfa.html#functions-7) - [Estimation](https://m-clark.github.io/models-by-example/cfa.html#estimation-21) - [Raw](https://m-clark.github.io/models-by-example/cfa.html#raw) - [Standardized](https://m-clark.github.io/models-by-example/cfa.html#standardized) - [Comparison](https://m-clark.github.io/models-by-example/cfa.html#comparison-27) - [Mplus](https://m-clark.github.io/models-by-example/cfa.html#mplus) - [Source](https://m-clark.github.io/models-by-example/cfa.html#source-23) - Bayesian - [Introduction to Bayesian Methods](https://m-clark.github.io/models-by-example/bayesian.html) - [Basics](https://m-clark.github.io/models-by-example/bayesian-basics.html) - [Prior, likelihood, & posterior distributions](https://m-clark.github.io/models-by-example/bayesian-basics.html#prior-likelihood-posterior-distributions) - [Prior](https://m-clark.github.io/models-by-example/bayesian-basics.html#prior) - [Likelihood](https://m-clark.github.io/models-by-example/bayesian-basics.html#likelihood) - [Posterior](https://m-clark.github.io/models-by-example/bayesian-basics.html#posterior) - [Posterior predictive](https://m-clark.github.io/models-by-example/bayesian-basics.html#posterior-predictive) - [Source](https://m-clark.github.io/models-by-example/bayesian-basics.html#source-24) - [Bayesian t-test](https://m-clark.github.io/models-by-example/bayesian-t-test.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-t-test.html#data-setup-24) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-t-test.html#model-code) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-t-test.html#estimation-22) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-t-test.html#comparison-28) - [Visualization](https://m-clark.github.io/models-by-example/bayesian-t-test.html#visualization) - [Source](https://m-clark.github.io/models-by-example/bayesian-t-test.html#source-25) - [Bayesian Linear Regression](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html#data-setup-25) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html#model-code-1) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html#estimation-23) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html#comparison-29) - [Visualize](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html#visualize-1) - [Source](https://m-clark.github.io/models-by-example/bayesian-linear-regression.html#source-26) - [Bayesian Beta Regression](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html#data-setup-26) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html#model-code-2) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html#estimation-24) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html#comparison-30) - [Visualization](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html#visualization-1) - [Source](https://m-clark.github.io/models-by-example/bayesian-beta-regression.html#source-27) - [Bayesian Mixed Model](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html#data-setup-27) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html#model-code-3) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html#estimation-25) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html#comparison-31) - [Visualize](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html#visualize-2) - [Source](https://m-clark.github.io/models-by-example/bayesian-mixed-model.html#source-28) - [Bayesian Multilevel Mediation](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html#data-setup-28) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html#model-code-4) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html#estimation-26) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html#comparison-32) - [Visualization](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html#visualization-2) - [Source](https://m-clark.github.io/models-by-example/bayesian-mixed-mediation.html#source-29) - [Bayesian IRT](https://m-clark.github.io/models-by-example/bayesian-irt.html) - [One Parameter IRT](https://m-clark.github.io/models-by-example/bayesian-irt.html#one-parameter-irt) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-irt.html#data-setup-29) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-irt.html#model-code-5) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-irt.html#estimation-27) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-irt.html#comparison-33) - [Two Parameter IRT](https://m-clark.github.io/models-by-example/bayesian-irt.html#two-parameter-irt) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-irt.html#model-code-6) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-irt.html#estimation-28) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-irt.html#comparison-34) - [Three Parameter IRT](https://m-clark.github.io/models-by-example/bayesian-irt.html#three-parameter-irt) - [Four Parameter IRT](https://m-clark.github.io/models-by-example/bayesian-irt.html#four-parameter-irt) - [Source](https://m-clark.github.io/models-by-example/bayesian-irt.html#source-30) - [Bayesian CFA](https://m-clark.github.io/models-by-example/bayesian-cfa.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-cfa.html#data-setup-30) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-cfa.html#model-code-7) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-cfa.html#estimation-29) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-cfa.html#comparison-35) - [Source](https://m-clark.github.io/models-by-example/bayesian-cfa.html#source-31) - [Bayesian Nonparametric Models](https://m-clark.github.io/models-by-example/bayesian-non-parametric.html) - [Chinese Restaurant Process](https://m-clark.github.io/models-by-example/bayesian-non-parametric.html#chinese-restaurant-process) - [Indian Buffet Process](https://m-clark.github.io/models-by-example/bayesian-non-parametric.html#indian-buffet-process) - [Source](https://m-clark.github.io/models-by-example/bayesian-non-parametric.html#source-32) - [Bayesian Stochastic Volatility Model](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html#data-setup-31) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html#model-code-8) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html#estimation-30) - [Results](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html#results) - [Visualization](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html#visualization-3) - [Source](https://m-clark.github.io/models-by-example/bayesian-stochastic-volatility.html#source-33) - [Bayesian Multinomial Models](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html#data-setup-32) - [Model Code](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html#model-code-9) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html#estimation-31) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html#comparison-36) - [Adding Complexity](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html#adding-complexity) - [Source](https://m-clark.github.io/models-by-example/bayesian-multinomial-model.html#source-34) - [Variational Bayes Regression](https://m-clark.github.io/models-by-example/bayesian-variational.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-variational.html#data-setup-33) - [Function](https://m-clark.github.io/models-by-example/bayesian-variational.html#function-19) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-variational.html#estimation-32) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-variational.html#comparison-37) - [Visualization](https://m-clark.github.io/models-by-example/bayesian-variational.html#visualization-4) - [Supplemental Example](https://m-clark.github.io/models-by-example/bayesian-variational.html#supplemental-example-2) - [Source](https://m-clark.github.io/models-by-example/bayesian-variational.html#source-35) - [Topic Model](https://m-clark.github.io/models-by-example/bayesian-topic-model.html) - [Data Setup](https://m-clark.github.io/models-by-example/bayesian-topic-model.html#data-setup-34) - [Function](https://m-clark.github.io/models-by-example/bayesian-topic-model.html#function-20) - [Estimation](https://m-clark.github.io/models-by-example/bayesian-topic-model.html#estimation-33) - [Comparison](https://m-clark.github.io/models-by-example/bayesian-topic-model.html#comparison-38) - [Source](https://m-clark.github.io/models-by-example/bayesian-topic-model.html#source-36) - Estimation - [Maximum Likelihood](https://m-clark.github.io/models-by-example/maximum-likelihood.html) - [Example](https://m-clark.github.io/models-by-example/maximum-likelihood.html#example-4) - [Linear Model](https://m-clark.github.io/models-by-example/maximum-likelihood.html#linear-model) - [Source](https://m-clark.github.io/models-by-example/maximum-likelihood.html#source-37) - [Penalized Maximum Likelihood](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html) - [Introduction](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#introduction-2) - [Data Setup](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#data-setup-35) - [Functions](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#functions-8) - [Estimation](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#estimation-34) - [Comparison](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#comparison-39) - [Source](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#source-38) - [L1 (lasso) regularization](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#l1-lasso-regularization) - [Data Setup](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#data-setup-36) - [Function](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#function-21) - [Estimation](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#estimation-35) - [Comparison](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#comparison-40) - [Source](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#source-39) - [L2 (ridge) regularization](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#l2-ridge-regularization) - [Data Setup](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#data-setup-37) - [Function](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#function-22) - [Estimation](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#estimation-36) - [Comparison](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#comparison-41) - [Source](https://m-clark.github.io/models-by-example/penalized-maximum-likelihood.html#source-40) - [Newton and IRLS](https://m-clark.github.io/models-by-example/newton-irls.html) - [Data Setup](https://m-clark.github.io/models-by-example/newton-irls.html#data-setup-38) - [Functions](https://m-clark.github.io/models-by-example/newton-irls.html#functions-9) - [Newton’s Method](https://m-clark.github.io/models-by-example/newton-irls.html#newtons-method) - [IRLS](https://m-clark.github.io/models-by-example/newton-irls.html#irls) - [Estimation](https://m-clark.github.io/models-by-example/newton-irls.html#estimation-37) - [Comparison](https://m-clark.github.io/models-by-example/newton-irls.html#comparison-42) - [Source](https://m-clark.github.io/models-by-example/newton-irls.html#source-41) - [Nelder-Mead](https://m-clark.github.io/models-by-example/nelder-mead.html) - [Functions](https://m-clark.github.io/models-by-example/nelder-mead.html#functions-10) - [First Version](https://m-clark.github.io/models-by-example/nelder-mead.html#first-version) - [Example](https://m-clark.github.io/models-by-example/nelder-mead.html#example-5) - [A Regression Model](https://m-clark.github.io/models-by-example/nelder-mead.html#a-regression-model) - [Second Version](https://m-clark.github.io/models-by-example/nelder-mead.html#second-version) - [Example](https://m-clark.github.io/models-by-example/nelder-mead.html#example-6) - [A Regression Model](https://m-clark.github.io/models-by-example/nelder-mead.html#a-regression-model-1) - [Source](https://m-clark.github.io/models-by-example/nelder-mead.html#source-42) - [Expectation-Maximization](https://m-clark.github.io/models-by-example/em.html) - [Mixture Model](https://m-clark.github.io/models-by-example/em.html#mixture-model) - [Data Setup](https://m-clark.github.io/models-by-example/em.html#data-setup-40) - [Function](https://m-clark.github.io/models-by-example/em.html#function-23) - [Estimation](https://m-clark.github.io/models-by-example/em.html#estimation-38) - [Comparison](https://m-clark.github.io/models-by-example/em.html#comparison-45) - [Visualization](https://m-clark.github.io/models-by-example/em.html#visualization-5) - [Supplemental Example](https://m-clark.github.io/models-by-example/em.html#supplemental-example-3) - [Source](https://m-clark.github.io/models-by-example/em.html#source-43) - [Multivariate Mixture Model](https://m-clark.github.io/models-by-example/em.html#multivariate-mixture-model) - [Function](https://m-clark.github.io/models-by-example/em.html#function-24) - [Example 1: Old Faithful](https://m-clark.github.io/models-by-example/em.html#example-1-old-faithful) - [Example 2: Iris](https://m-clark.github.io/models-by-example/em.html#example-2-iris) - [Source](https://m-clark.github.io/models-by-example/em.html#source-44) - [Probit Model](https://m-clark.github.io/models-by-example/em.html#probit-model) - [Data Setup](https://m-clark.github.io/models-by-example/em.html#data-setup-43) - [Probit via Maximum Likelihood](https://m-clark.github.io/models-by-example/em.html#probit-via-maximum-likelihood) - [EM for Latent Variable Approach to Probit](https://m-clark.github.io/models-by-example/em.html#em-for-latent-variable-approach-to-probit) - [Source](https://m-clark.github.io/models-by-example/em.html#source-45) - [PCA](https://m-clark.github.io/models-by-example/em.html#pca) - [Data Setup](https://m-clark.github.io/models-by-example/em.html#data-setup-44) - [Function](https://m-clark.github.io/models-by-example/em.html#function-27) - [Estimation](https://m-clark.github.io/models-by-example/em.html#estimation-43) - [Comparison](https://m-clark.github.io/models-by-example/em.html#comparison-50) - [Visualize](https://m-clark.github.io/models-by-example/em.html#visualize-4) - [Source](https://m-clark.github.io/models-by-example/em.html#source-46) - [Probabilistic PCA](https://m-clark.github.io/models-by-example/em.html#probabilistic-pca) - [Data Setup](https://m-clark.github.io/models-by-example/em.html#data-setup-45) - [Function](https://m-clark.github.io/models-by-example/em.html#function-28) - [Estimation](https://m-clark.github.io/models-by-example/em.html#estimation-44) - [Comparison](https://m-clark.github.io/models-by-example/em.html#comparison-51) - [Visualize](https://m-clark.github.io/models-by-example/em.html#visualize-5) - [Missing Data Example](https://m-clark.github.io/models-by-example/em.html#missing-data-example) - [Source](https://m-clark.github.io/models-by-example/em.html#source-47) - [State Space Model](https://m-clark.github.io/models-by-example/em.html#state-space-model) - [Data Setup](https://m-clark.github.io/models-by-example/em.html#data-setup-47) - [Function](https://m-clark.github.io/models-by-example/em.html#function-30) - [Estimation](https://m-clark.github.io/models-by-example/em.html#estimation-46) - [Visualization](https://m-clark.github.io/models-by-example/em.html#visualization-6) - [Source](https://m-clark.github.io/models-by-example/em.html#source-48) - [Gradient Descent](https://m-clark.github.io/models-by-example/gradient-descent.html) - [Data Setup](https://m-clark.github.io/models-by-example/gradient-descent.html#data-setup-48) - [Function](https://m-clark.github.io/models-by-example/gradient-descent.html#function-31) - [Estimation](https://m-clark.github.io/models-by-example/gradient-descent.html#estimation-47) - [Comparison](https://m-clark.github.io/models-by-example/gradient-descent.html#comparison-53) - [Source](https://m-clark.github.io/models-by-example/gradient-descent.html#source-49) - [Stochastic Gradient Descent](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html) - [Data Setup](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#data-setup-49) - [Function](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#function-32) - [Estimation](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#estimation-48) - [Comparison](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#comparison-54) - [Visualize Estimates](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#visualize-estimates) - [Data Set Shift](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#data-set-shift) - [Data Setup](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#data-setup-50) - [Estimation](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#estimation-49) - [Comparison](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#comparison-55) - [Visualize Estimates](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#visualize-estimates-1) - [SGD Variants](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#sgd-variants) - [Data Setup](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#data-setup-51) - [Function](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#function-33) - [Estimation](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#estimation-50) - [Comparison](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#comparison-56) - [Visualize Estimates](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#visualize-estimates-2) - [Source](https://m-clark.github.io/models-by-example/stochastic-gradient-descent.html#source-50) - [Metropolis Hastings](https://m-clark.github.io/models-by-example/metropolis-hastings.html) - [Data Setup](https://m-clark.github.io/models-by-example/metropolis-hastings.html#data-setup-52) - [Functions](https://m-clark.github.io/models-by-example/metropolis-hastings.html#functions-11) - [Estimation](https://m-clark.github.io/models-by-example/metropolis-hastings.html#estimation-51) - [Comparison](https://m-clark.github.io/models-by-example/metropolis-hastings.html#comparison-57) - [Source](https://m-clark.github.io/models-by-example/metropolis-hastings.html#source-51) - [Hamiltonian Monte Carlo](https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html) - [Data Setup](https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html#data-setup-53) - [Functions](https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html#functions-12) - [Estimation](https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html#estimation-52) - [Comparison](https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html#comparison-58) - [Source](https://m-clark.github.io/models-by-example/hamiltonian-monte-carlo.html#source-52) - Other - [Supplemental](https://m-clark.github.io/models-by-example/supplemental.html) - [Other Languages](https://m-clark.github.io/models-by-example/supplemental.html#other-languages) - [Python Demos](https://m-clark.github.io/models-by-example/supplemental.html#python-demos) - [Linear Regression](https://m-clark.github.io/models-by-example/supplemental.html#python-linreg) - [Logistic Regression](https://m-clark.github.io/models-by-example/supplemental.html#python-logreg) - [Quantile Regression](https://m-clark.github.io/models-by-example/supplemental.html#python-qreg) - [Stochastic Gradient Descent](https://m-clark.github.io/models-by-example/supplemental.html#python-sgd) - [Nelder-Mead](https://m-clark.github.io/models-by-example/supplemental.html#python-nelder) - [HMM](https://m-clark.github.io/models-by-example/supplemental.html#python-hmm) - [Julia Demos](https://m-clark.github.io/models-by-example/supplemental.html#julia-demos) - [Mixed Models](https://m-clark.github.io/models-by-example/supplemental.html#mixed-models-1) - [Matlab Demos](https://m-clark.github.io/models-by-example/supplemental.html#matlab-demos) - [Mixed Models](https://m-clark.github.io/models-by-example/supplemental.html#mixed-models-2) - [Gaussian Processes](https://m-clark.github.io/models-by-example/supplemental.html#matlab-gp) - [![MC logo](https://m-clark.github.io/models-by-example/img/mc_logo.png)](https://m-clark.github.io/) # [Model Estimation by Example](https://m-clark.github.io/models-by-example/) # Zero-Inflated Model Log likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples and comparison to pscl package output. Also includes approach based on Hilbe GLM text. Zero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory. ## Poisson ### Data Setup Get the data. ``` library(tidyverse) fish = haven::read_dta("http://www.stata-press.com/data/r11/fish.dta") ``` ### Function The log likelihood function. ``` zip_ll <- function(y, X, par) { # arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois' # Extract parameters logitpars = par[grep('logit', names(par))] poispars = par[grep('pois', names(par))] # Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example Xlogit = X LPlogit = Xlogit %% logitpars logi0 = plogis(LPlogit) # alternative 1/(1+exp(-LPlogit)) # Poisson part Xpois = X mupois = exp(Xpois %% poispars) # LLs logliklogit = log( logi0 + exp(log(1 - logi0) - mupois) ) loglikpois = log(1 - logi0) + dpois(y, lambda = mupois, log = TRUE) # Hilbe formulation # logliklogit = log(logi0 + (1 - logi0)exp(- mupois) ) # loglikpois = log(1-logi0) -mupois + log(mupois)y #not necessary: - log(gamma(y+1)) y0 = y == 0 # 0 values yc = y > 0 # Count part loglik = sum(logliklogit[y0]) + sum(loglikpois[yc]) -loglik } ``` ### Estimation Get starting values or simply do zeros. ``` # for zip: need 'logit', 'pois' initial_model = glm( count ~ persons + livebait, data = fish, x = TRUE, y = TRUE, "poisson" ) # starts = c(logit = coef(initial_model), pois = coef(initial_model)) starts = c(rep(0, 3), rep(0, 3)) names(starts) = c(paste0('pois.', names(coef(initial_model))), paste0('logit.', names(coef(initial_model)))) ``` Estimate with optim. ``` fit = optim( par = starts , fn = zip_ll, X = initial_model$x, y = initial_model$y, method = "BFGS", control = list(maxit = 5000, reltol = 1e-12), hessian = TRUE ) # fit ``` ### Comparison Extract for clean display. ``` B = fit$par se = sqrt(diag(solve((fit$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)2 ``` Results from pscl. ``` library(pscl) fit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = "poisson") ``` Compare. \| \| coef \| B \| se \| Z \| p \| \|---\|---\|---\|---\|---\|---\| \| pscl \| X.Intercept. \| \-2.006 \| 0\.324 \| \-6.196 \| 0\.000 \| \| pscl \| persons \| 0\.747 \| 0\.043 \| 17\.516 \| 0\.000 \| \| pscl \| livebait \| 1\.809 \| 0\.292 \| 6\.195 \| 0\.000 \| \| pscl \| X.Intercept..1 \| 0\.303 \| 0\.674 \| 0\.449 \| 0\.654 \| \| pscl \| persons.1 \| \-0.069 \| 0\.129 \| \-0.537 \| 0\.591 \| \| pscl \| livebait.1 \| \-0.031 \| 0\.558 \| \-0.056 \| 0\.956 \| \| zip\_poisson\_ll \| pois.(Intercept) \| \-2.006 \| 0\.324 \| \-6.196 \| 0\.000 \| \| zip\_poisson\_ll \| pois.persons \| 0\.747 \| 0\.043 \| 17\.516 \| 0\.000 \| \| zip\_poisson\_ll \| pois.livebait \| 1\.809 \| 0\.292 \| 6\.195 \| 0\.000 \| \| zip\_poisson\_ll \| logit.(Intercept) \| 0\.303 \| 0\.674 \| 0\.449 \| 0\.654 \| \| zip\_poisson\_ll \| logit.persons \| \-0.069 \| 0\.129 \| \-0.537 \| 0\.591 \| \| zip\_poisson\_ll \| logit.livebait \| \-0.031 \| 0\.558 \| \-0.056 \| 0\.956 \| ## Negative Binomial ### Function ``` zinb_ll <- function(y, X, par) { # arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta' # Extract parameters logitpars = par[grep('logit', names(par))] negbinpars = par[grep('negbin', names(par))] theta = exp(par[grep('theta', names(par))]) # Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example Xlogit = X LPlogit = Xlogit %% logitpars logi0 = plogis(LPlogit) # Negbin part Xnegbin = X munb = exp(Xnegbin %% negbinpars) # LLs logliklogit = log( logi0 + exp( log(1 - logi0) + suppressWarnings(dnbinom(0, size = theta, mu = munb, log = TRUE)) ) ) logliknegbin = log(1 - logi0) + suppressWarnings(dnbinom(y, size = theta, mu = munb, log = TRUE)) # Hilbe formulation # theta part # alpha = 1/theta # m = 1/alpha # p = 1/(1 + alphamunb) # logliklogit = log( logi0 + (1 - logi0)(p^m) ) # logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + mlog(p) + ylog(1-p) # gamma(y+1) not needed y0 = y == 0 # 0 values yc = y > 0 # Count part loglik = sum(logliklogit[y0]) + sum(logliknegbin[yc]) -loglik } ``` ### Estimation Get starting values or simply do zeros. ``` # for zinb: 'logit', 'negbin', 'theta' initial_model = model.matrix(count ~ persons + livebait, data = fish) # to get X matrix startlogi = glm(count == 0 ~ persons + livebait, data = fish, family = "binomial") startcount = glm(count ~ persons + livebait, data = fish, family = "poisson") starts = c( negbin = coef(startcount), logit = coef(startlogi), theta = 1 ) # starts = c(negbin = rep(0, 3), # logit = rep(0, 3), # theta = log(1)) ``` Estimate with optim. ``` fit_nb = optim( par = starts , fn = zinb_ll, X = initial_model, y = fish$count, method = "BFGS", control = list(maxit = 5000, reltol = 1e-12), hessian = TRUE ) # fit_nb ``` ### Comparison Extract for clean display. ``` B = fit_nb$par se = sqrt(diag(solve((fit_nb$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)2 ``` Results from pscl. ``` # pscl results library(pscl) fit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = "negbin") ``` Compare. \| \| coef \| B \| se \| Z \| p \| \|---\|---\|---\|---\|---\|---\| \| pscl \| X.Intercept. \| \-2.803 \| 0\.558 \| \-5.026 \| 0\.000 \| \| pscl \| persons \| 0\.849 \| 0\.124 \| 6\.833 \| 0\.000 \| \| pscl \| livebait \| 1\.791 \| 0\.511 \| 3\.504 \| 0\.000 \| \| pscl \| Log.theta. \| \-0.969 \| 0\.302 \| \-3.206 \| 0\.001 \| \| pscl \| X.Intercept..1 \| \-4.276 \| 4\.278 \| \-1.000 \| 0\.318 \| \| pscl \| persons.1 \| 0\.560 \| 0\.517 \| 1\.084 \| 0\.279 \| \| pscl \| livebait.1 \| 1\.168 \| 3\.661 \| 0\.319 \| 0\.750 \| \| zinb\_ll \| negbin.(Intercept) \| \-2.803 \| 0\.558 \| \-5.026 \| 0\.000 \| \| zinb\_ll \| negbin.persons \| 0\.849 \| 0\.124 \| 6\.834 \| 0\.000 \| \| zinb\_ll \| negbin.livebait \| 1\.791 \| 0\.511 \| 3\.504 \| 0\.000 \| \| zinb\_ll \| logit.(Intercept) \| \-4.276 \| 4\.278 \| \-1.000 \| 0\.318 \| \| zinb\_ll \| logit.persons \| 0\.560 \| 0\.517 \| 1\.084 \| 0\.279 \| \| zinb\_ll \| logit.livebait \| 1\.168 \| 3\.661 \| 0\.319 \| 0\.750 \| \| zinb\_ll \| theta \| \-0.969 \| 0\.302 \| \-3.206 \| 0\.001 \| ## Supplemental Example This supplemental example uses the bioChemists data. It contains a sample of 915 biochemistry graduate students with the following information: - art: count of articles produced during last 3 years of Ph.D. - fem: factor indicating gender of student, with levels Men and Women - mar: factor indicating marital status of student, with levels Single and Married - kid5: number of children aged 5 or younger - phd: prestige of Ph.D. department - ment: count of articles produced by Ph.D. mentor during last 3 years ``` data("bioChemists", package = "pscl") initial_model = model.matrix(art ~ fem + mar + kid5 + phd + ment, data = bioChemists) # to get X matrix startlogi = glm(art==0 ~ fem + mar + kid5 + phd + ment, data = bioChemists, family = "binomial") startcount = glm(art ~ fem + mar + kid5 + phd + ment, data = bioChemists, family = "quasipoisson") starts = c( negbin = coef(startcount), logit = coef(startlogi), theta = summary(startcount)$dispersion ) # starts = c(negbin = rep(0, 6), # logit = rep(0, 6), # theta = 1) fit_nb_pub = optim( par = starts , fn = zinb_ll, X = initial_model, y = bioChemists$art, method = "BFGS", control = list(maxit = 5000, reltol = 1e-12), hessian = TRUE ) # fit_nb_pub B = fit_nb_pub$par se = sqrt(diag(solve((fit_nb_pub$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)*2 library(pscl) fit_pscl = zeroinfl(art ~ . \| ., data = bioChemists, dist = "negbin") summary(fit_pscl)$coefficients ``` ``` $count Estimate Std. Error z value Pr(>\|z\|) (Intercept) 0.4167465901 0.143596450 2.90220678 3.705439e-03 femWomen -0.1955076374 0.075592558 -2.58633447 9.700275e-03 marMarried 0.0975826042 0.084451953 1.15548073 2.478936e-01 kid5 -0.1517320709 0.054206071 -2.79917119 5.123397e-03 phd -0.0006997593 0.036269674 -0.01929323 9.846072e-01 ment 0.0247861500 0.003492672 7.09661548 1.278491e-12 Log(theta) 0.9763577454 0.135469554 7.20721163 5.710921e-13 $zero Estimate Std. Error z value Pr(>\|z\|) (Intercept) -0.19160645 1.3227962 -0.1448496 0.884829645 femWomen 0.63587048 0.8488959 0.7490559 0.453823498 marMarried -1.49943716 0.9386562 -1.5974296 0.110169987 kid5 0.62840922 0.4427746 1.4192531 0.155825245 phd -0.03773288 0.3080059 -0.1225070 0.902497523 ment -0.88227364 0.3162186 -2.7900755 0.005269575 ``` ``` round(data.frame(B, se, Z, p), 4) ``` ``` B se Z p negbin.(Intercept) 0.4167 0.1436 2.9021 0.0037 negbin.femWomen -0.1955 0.0756 -2.5863 0.0097 negbin.marMarried 0.0976 0.0845 1.1555 0.2479 negbin.kid5 -0.1517 0.0542 -2.7992 0.0051 negbin.phd -0.0007 0.0363 -0.0195 0.9845 negbin.ment 0.0248 0.0035 7.0967 0.0000 logit.(Intercept) -0.1916 1.3229 -0.1449 0.8848 logit.femWomen 0.6359 0.8490 0.7491 0.4538 logit.marMarried -1.4995 0.9387 -1.5974 0.1102 logit.kid5 0.6284 0.4428 1.4192 0.1558 logit.phd -0.0377 0.3080 -0.1225 0.9025 logit.ment -0.8823 0.3162 -2.7901 0.0053 theta 0.9763 0.1355 7.2071 0.0000 ``` ## Source Original code for zip\_ll found at <https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/poiszeroinfl.R> Original code for zinb\_ll found at <https://github.com/m-clark/Miscellaneous-R-Code/blob/master/ModelFitting/NBzeroinfl.R>
Readable Markdown	Log likelihood function to estimate parameters for a Zero-inflated Poisson model. With examples and comparison to pscl package output. Also includes approach based on Hilbe GLM text. Zero-inflated models are applied to situations in which target data has relatively many of one value, usually zero, to go along with the other observed values. They are two-part models, a logistic model for whether an observation is zero or not, and a count model for the other part. The key distinction from hurdle count models is that the count distribution contributes to the excess zeros. While the typical application is count data, the approach can be applied to any distribution in theory. ## Poisson ### Data Setup Get the data. ``` library(tidyverse) fish = haven::read_dta("http://www.stata-press.com/data/r11/fish.dta") ``` ### Function The log likelihood function. ``` zip_ll <- function(y, X, par) { # arguments are response y, predictor matrix X, and parameter named starting points of 'logit' and 'pois' # Extract parameters logitpars = par[grep('logit', names(par))] poispars = par[grep('pois', names(par))] # Logit part; in this function Xlogit = Xpois but one could split X argument into Xlogi and Xpois for example Xlogit = X LPlogit = Xlogit %% logitpars logi0 = plogis(LPlogit) # alternative 1/(1+exp(-LPlogit)) # Poisson part Xpois = X mupois = exp(Xpois %% poispars) # LLs logliklogit = log( logi0 + exp(log(1 - logi0) - mupois) ) loglikpois = log(1 - logi0) + dpois(y, lambda = mupois, log = TRUE) # Hilbe formulation # logliklogit = log(logi0 + (1 - logi0)exp(- mupois) ) # loglikpois = log(1-logi0) -mupois + log(mupois)y #not necessary: - log(gamma(y+1)) y0 = y == 0 # 0 values yc = y > 0 # Count part loglik = sum(logliklogit[y0]) + sum(loglikpois[yc]) -loglik } ``` ### Estimation Get starting values or simply do zeros. ``` # for zip: need 'logit', 'pois' initial_model = glm( count ~ persons + livebait, data = fish, x = TRUE, y = TRUE, "poisson" ) # starts = c(logit = coef(initial_model), pois = coef(initial_model)) starts = c(rep(0, 3), rep(0, 3)) names(starts) = c(paste0('pois.', names(coef(initial_model))), paste0('logit.', names(coef(initial_model)))) ``` Estimate with optim. ``` fit = optim( par = starts , fn = zip_ll, X = initial_model$x, y = initial_model$y, method = "BFGS", control = list(maxit = 5000, reltol = 1e-12), hessian = TRUE ) # fit ``` ### Comparison Extract for clean display. ``` B = fit$par se = sqrt(diag(solve((fit$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)2 ``` Results from pscl. ``` library(pscl) fit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = "poisson") ``` Compare. \| \| coef \| B \| se \| Z \| p \| \|---\|---\|---\|---\|---\|---\| \| pscl \| X.Intercept. \| \-2.006 \| 0\.324 \| \-6.196 \| 0\.000 \| \| pscl \| persons \| 0\.747 \| 0\.043 \| 17\.516 \| 0\.000 \| \| pscl \| livebait \| 1\.809 \| 0\.292 \| 6\.195 \| 0\.000 \| \| pscl \| X.Intercept..1 \| 0\.303 \| 0\.674 \| 0\.449 \| 0\.654 \| \| pscl \| persons.1 \| \-0.069 \| 0\.129 \| \-0.537 \| 0\.591 \| \| pscl \| livebait.1 \| \-0.031 \| 0\.558 \| \-0.056 \| 0\.956 \| \| zip\_poisson\_ll \| pois.(Intercept) \| \-2.006 \| 0\.324 \| \-6.196 \| 0\.000 \| \| zip\_poisson\_ll \| pois.persons \| 0\.747 \| 0\.043 \| 17\.516 \| 0\.000 \| \| zip\_poisson\_ll \| pois.livebait \| 1\.809 \| 0\.292 \| 6\.195 \| 0\.000 \| \| zip\_poisson\_ll \| logit.(Intercept) \| 0\.303 \| 0\.674 \| 0\.449 \| 0\.654 \| \| zip\_poisson\_ll \| logit.persons \| \-0.069 \| 0\.129 \| \-0.537 \| 0\.591 \| \| zip\_poisson\_ll \| logit.livebait \| \-0.031 \| 0\.558 \| \-0.056 \| 0\.956 \| ## Negative Binomial ### Function ``` zinb_ll <- function(y, X, par) { # arguments are response y, predictor matrix X, and parameter named starting points of 'logit', 'negbin', and 'theta' # Extract parameters logitpars = par[grep('logit', names(par))] negbinpars = par[grep('negbin', names(par))] theta = exp(par[grep('theta', names(par))]) # Logit part; in this function Xlogit = Xnegbin but one could split X argument into Xlogit and Xnegbin for example Xlogit = X LPlogit = Xlogit %% logitpars logi0 = plogis(LPlogit) # Negbin part Xnegbin = X munb = exp(Xnegbin %% negbinpars) # LLs logliklogit = log( logi0 + exp( log(1 - logi0) + suppressWarnings(dnbinom(0, size = theta, mu = munb, log = TRUE)) ) ) logliknegbin = log(1 - logi0) + suppressWarnings(dnbinom(y, size = theta, mu = munb, log = TRUE)) # Hilbe formulation # theta part # alpha = 1/theta # m = 1/alpha # p = 1/(1 + alphamunb) # logliklogit = log( logi0 + (1 - logi0)(p^m) ) # logliknegbin = log(1-logi0) + log(gamma(m+y)) - log(gamma(m)) + mlog(p) + ylog(1-p) # gamma(y+1) not needed y0 = y == 0 # 0 values yc = y > 0 # Count part loglik = sum(logliklogit[y0]) + sum(logliknegbin[yc]) -loglik } ``` ### Estimation Get starting values or simply do zeros. ``` # for zinb: 'logit', 'negbin', 'theta' initial_model = model.matrix(count ~ persons + livebait, data = fish) # to get X matrix startlogi = glm(count == 0 ~ persons + livebait, data = fish, family = "binomial") startcount = glm(count ~ persons + livebait, data = fish, family = "poisson") starts = c( negbin = coef(startcount), logit = coef(startlogi), theta = 1 ) # starts = c(negbin = rep(0, 3), # logit = rep(0, 3), # theta = log(1)) ``` Estimate with optim. ``` fit_nb = optim( par = starts , fn = zinb_ll, X = initial_model, y = fish$count, method = "BFGS", control = list(maxit = 5000, reltol = 1e-12), hessian = TRUE ) # fit_nb ``` ### Comparison Extract for clean display. ``` B = fit_nb$par se = sqrt(diag(solve((fit_nb$hessian)))) Z = B/se p = pnorm(abs(Z), lower = FALSE)2 ``` Results from pscl. ``` # pscl results library(pscl) fit_pscl = zeroinfl(count ~ persons + livebait, data = fish, dist = "negbin") ``` Compare. \| \| coef \| B \| se \| Z \| p \| \|---\|---\|---\|---\|---\|---\| \| pscl \| X.Intercept. \| \-2.803 \| 0\.558 \| \-5.026 \| 0\.000 \| \| pscl \| persons \| 0\.849 \| 0\.124 \| 6\.833 \| 0\.000 \| \| pscl \| livebait \| 1\.791 \| 0\.511 \| 3\.504 \| 0\.000 \| \| pscl \| Log.theta. \| \-0.969 \| 0\.302 \| \-3.206 \| 0\.001 \| \| pscl \| X.Intercept..1 \| \-4.276 \| 4\.278 \| \-1.000 \| 0\.318 \| \| pscl \| persons.1 \| 0\.560 \| 0\.517 \| 1\.084 \| 0\.279 \| \| pscl \| livebait.1 \| 1\.168 \| 3\.661 \| 0\.319 \| 0\.750 \| \| zinb\_ll \| negbin.(Intercept) \| \-2.803 \| 0\.558 \| \-5.026 \| 0\.000 \| \| zinb\_ll \| negbin.persons \| 0\.849 \| 0\.124 \| 6\.834 \| 0\.000 \| \| zinb\_ll \| negbin.livebait \| 1\.791 \| 0\.511 \| 3\.504 \| 0\.000 \| \| zinb\_ll \| logit.(Intercept) \| \-4.276 \| 4\.278 \| \-1.000 \| 0\.318 \| \| zinb\_ll \| logit.persons \| 0\.560 \| 0\.517 \| 1\.084 \| 0\.279 \| \| zinb\_ll \| logit.livebait \| 1\.168 \| 3\.661 \| 0\.319 \| 0\.750 \| \| zinb\_ll \| theta \| \-0.969 \| 0\.302 \| \-3.206 \| 0\.001 \|
Shard	143 (laksa)
Root Hash	2566890010099092343
Unparsed URL	io,github!m-clark,/models-by-example/zi.html s443