{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "# Beetle Mortality Data\n", "Author(s): Paul Miles | Date Created: June 18, 2018\n", "\n", "This is a binomial logistic regression example with dose-response data.\n", "\n", "A. Dobson, \"An Introduction to Generalized Linear Models\", Chapman & Hall/CRC, 2002." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9.0\n" ] } ], "source": [ "# import required packages\n", "import numpy as np\n", "import math\n", "from pymcmcstat.MCMC import MCMC\n", "import matplotlib.pyplot as plt\n", "import pymcmcstat\n", "print(pymcmcstat.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define mortality data" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Beetle mortality data\n", "dose = np.array([1.6907, 1.7242, 1.7552, 1.7842,\n", " 1.8113, 1.8369, 1.8610, 1.8839])\n", "number_of_beetles = np.array([59, 60, 62, 56, 63, 59, 62, 60])\n", "number_of_beetles_killed = np.array([6, 13, 18, 28, 52, 53, 61, 60])\n", "\n", "x = np.array([dose, number_of_beetles]).T\n", "y = number_of_beetles_killed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Define model options\n", "The user is given several modeling options. The user chooses the desired options by defining the appropriate keywork to `beetle_link`.\n", "\n", "- `logitmodelfun`: \n", "$$y(d,\\theta) = \\frac{1}{1 + \\exp(\\theta_0 + \\theta_1d)}$$\n", "- `loglogmodelfun`:\n", "$$y(d,\\theta) = 1 - \\exp(-\\exp(\\theta_0 + \\theta_1d))$$\n", "- `probitmodelfun`:\n", "$$y(d,\\theta) = \\frac{1}{2}\\Big[1 + \\text{erf}\\Big(\\frac{\\theta_0 + \\theta_1d -\\mu}{\\sigma\\sqrt{2}}\\Big)\\Big], \\quad \\mu = 0, \\quad \\sigma^2 = 1$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def logitmodelfun(d, t):\n", " return 1/(1+np.exp(t[0]+t[1]*d))\n", "def loglogmodelfun(d, t):\n", " return 1 - np.exp(-np.exp(t[0] + t[1]*d))\n", "def nordf(x, mu = 0, sigma2 = 1):\n", " # NORDF the standard normal (Gaussian) cumulative distribution.\n", " # NORPF(x,mu,sigma2) x quantile, mu mean, sigma2 variance\n", " # Marko Laine