diff --git a/.quarto/idx/index.qmd.json b/.quarto/idx/index.qmd.json index 215fdbf..0a143e7 100644 --- a/.quarto/idx/index.qmd.json +++ b/.quarto/idx/index.qmd.json @@ -1 +1 @@ -{"title":"Course overview","markdown":{"yaml":{"title":"Course overview","number-sections":false},"headingText":"Core aims","containsRefs":false,"markdown":"\n\nWelcome to the wonderful world of generalised linear models!\n\nThese sessions are intended to enable you to construct and use generalised linear models confidently.\n\nAs with all of our statistics courses our focus is not on mathematical derivations, but on developing an intuitive understanding of the underlying statistical concepts.\n\nAt the same time this is also *not* a \"how to mindlessly use a stats program\" course! We hope that at the end of this course you feel like you have a better grasp on what it is we're trying to do, and gained sufficient confidence in your coding skills to implement these statistical concepts in your own research!\n\n\nTo introduce sufficient understanding and coding experience for analysing data with non-continuous response variables.\n\n::: callout-note\n## Course aims\n\nTo know what to do when presented with an arbitrary data set e.g.\n\n1. Construct\n a. a logistic model for binary response variables\n b. a logistic model for proportion response variables\n c. a Poisson model for count response variables\n d. ~~a Negative Binomial model for count response variables~~ (to be added later)\n2. Plot the data and the fitted curve in each case for both continuous and categorical predictors\n3. Assess the significance of fit\n4. Assess assumption of the model\n:::\n","srcMarkdownNoYaml":"\n\nWelcome to the wonderful world of generalised linear models!\n\nThese sessions are intended to enable you to construct and use generalised linear models confidently.\n\nAs with all of our statistics courses our focus is not on mathematical derivations, but on developing an intuitive understanding of the underlying statistical concepts.\n\nAt the same time this is also *not* a \"how to mindlessly use a stats program\" course! We hope that at the end of this course you feel like you have a better grasp on what it is we're trying to do, and gained sufficient confidence in your coding skills to implement these statistical concepts in your own research!\n\n## Core aims\n\nTo introduce sufficient understanding and coding experience for analysing data with non-continuous response variables.\n\n::: callout-note\n## Course aims\n\nTo know what to do when presented with an arbitrary data set e.g.\n\n1. Construct\n a. a logistic model for binary response variables\n b. a logistic model for proportion response variables\n c. a Poisson model for count response variables\n d. ~~a Negative Binomial model for count response variables~~ (to be added later)\n2. Plot the data and the fitted curve in each case for both continuous and categorical predictors\n3. Assess the significance of fit\n4. Assess assumption of the model\n:::\n"},"formats":{"courseformat-html":{"identifier":{"display-name":"HTML","target-format":"courseformat-html","base-format":"html","extension-name":"courseformat"},"execute":{"fig-width":7,"fig-height":5,"fig-format":"retina","fig-dpi":96,"df-print":"default","error":false,"eval":true,"cache":null,"freeze":"auto","echo":true,"output":true,"warning":true,"include":true,"keep-md":false,"keep-ipynb":false,"ipynb":null,"enabled":null,"daemon":null,"daemon-restart":false,"debug":false,"ipynb-filters":[],"ipynb-shell-interactivity":null,"plotly-connected":true,"engine":"markdown"},"render":{"keep-tex":false,"keep-typ":false,"keep-source":false,"keep-hidden":false,"prefer-html":false,"output-divs":true,"output-ext":"html","fig-align":"default","fig-pos":null,"fig-env":null,"code-fold":"none","code-overflow":"scroll","code-link":true,"code-line-numbers":false,"code-tools":false,"tbl-colwidths":"auto","merge-includes":true,"inline-includes":false,"preserve-yaml":false,"latex-auto-mk":true,"latex-auto-install":true,"latex-clean":true,"latex-min-runs":1,"latex-max-runs":10,"latex-makeindex":"makeindex","latex-makeindex-opts":[],"latex-tlmgr-opts":[],"latex-input-paths":[],"latex-output-dir":null,"link-external-icon":false,"link-external-newwindow":false,"self-contained-math":false,"format-resources":[],"notebook-links":true,"shortcodes":[]},"pandoc":{"standalone":true,"wrap":"none","default-image-extension":"png","to":"html","toc":true,"number-sections":false,"filters":["courseformat"],"output-file":"index.html"},"language":{"toc-title-document":"Table of contents","toc-title-website":"On this page","related-formats-title":"Other Formats","related-notebooks-title":"Notebooks","source-notebooks-prefix":"Source","other-links-title":"Other Links","code-links-title":"Code Links","launch-dev-container-title":"Launch Dev Container","launch-binder-title":"Launch Binder","article-notebook-label":"Article Notebook","notebook-preview-download":"Download Notebook","notebook-preview-download-src":"Download Source","notebook-preview-back":"Back to Article","manuscript-meca-bundle":"MECA Bundle","section-title-abstract":"Abstract","section-title-appendices":"Appendices","section-title-footnotes":"Footnotes","section-title-references":"References","section-title-reuse":"Reuse","section-title-copyright":"Copyright","section-title-citation":"Citation","appendix-attribution-cite-as":"For attribution, please cite this work as:","appendix-attribution-bibtex":"BibTeX citation:","title-block-author-single":"Author","title-block-author-plural":"Authors","title-block-affiliation-single":"Affiliation","title-block-affiliation-plural":"Affiliations","title-block-published":"Published","title-block-modified":"Modified","title-block-keywords":"Keywords","callout-tip-title":"Tip","callout-note-title":"Note","callout-warning-title":"Warning","callout-important-title":"Important","callout-caution-title":"Caution","code-summary":"Code","code-tools-menu-caption":"Code","code-tools-show-all-code":"Show All Code","code-tools-hide-all-code":"Hide All Code","code-tools-view-source":"View Source","code-tools-source-code":"Source Code","tools-share":"Share","tools-download":"Download","code-line":"Line","code-lines":"Lines","copy-button-tooltip":"Copy to Clipboard","copy-button-tooltip-success":"Copied!","repo-action-links-edit":"Edit this page","repo-action-links-source":"View source","repo-action-links-issue":"Report an issue","back-to-top":"Back to top","search-no-results-text":"No results","search-matching-documents-text":"matching documents","search-copy-link-title":"Copy link to search","search-hide-matches-text":"Hide additional matches","search-more-match-text":"more match in this document","search-more-matches-text":"more matches in this document","search-clear-button-title":"Clear","search-text-placeholder":"","search-detached-cancel-button-title":"Cancel","search-submit-button-title":"Submit","search-label":"Search","toggle-section":"Toggle section","toggle-sidebar":"Toggle sidebar navigation","toggle-dark-mode":"Toggle dark mode","toggle-reader-mode":"Toggle reader mode","toggle-navigation":"Toggle navigation","crossref-fig-title":"Figure","crossref-tbl-title":"Table","crossref-lst-title":"Listing","crossref-thm-title":"Theorem","crossref-lem-title":"Lemma","crossref-cor-title":"Corollary","crossref-prp-title":"Proposition","crossref-cnj-title":"Conjecture","crossref-def-title":"Definition","crossref-exm-title":"Example","crossref-exr-title":"Exercise","crossref-ch-prefix":"Chapter","crossref-apx-prefix":"Appendix","crossref-sec-prefix":"Section","crossref-eq-prefix":"Equation","crossref-lof-title":"List of Figures","crossref-lot-title":"List of Tables","crossref-lol-title":"List of Listings","environment-proof-title":"Proof","environment-remark-title":"Remark","environment-solution-title":"Solution","listing-page-order-by":"Order By","listing-page-order-by-default":"Default","listing-page-order-by-date-asc":"Oldest","listing-page-order-by-date-desc":"Newest","listing-page-order-by-number-desc":"High to Low","listing-page-order-by-number-asc":"Low to High","listing-page-field-date":"Date","listing-page-field-title":"Title","listing-page-field-description":"Description","listing-page-field-author":"Author","listing-page-field-filename":"File Name","listing-page-field-filemodified":"Modified","listing-page-field-subtitle":"Subtitle","listing-page-field-readingtime":"Reading Time","listing-page-field-wordcount":"Word Count","listing-page-field-categories":"Categories","listing-page-minutes-compact":"{0} min","listing-page-category-all":"All","listing-page-no-matches":"No matching items","listing-page-words":"{0} words"},"metadata":{"lang":"en","fig-responsive":true,"quarto-version":"1.4.546","theme":["default","_extensions/cambiotraining/courseformat/theme.scss"],"number-depth":3,"code-copy":true,"revealjs-plugins":[],"bibliography":["references.bib"],"knitr":{"opts_knit":{"cache.path":".knitr_cache"}},"title":"Course overview"},"extensions":{"book":{"multiFile":true}}}},"projectFormats":["courseformat-html"]} \ No newline at end of file +{"title":"Course overview","markdown":{"yaml":{"title":"Course overview","author":["Vicki Hodgson, Matt Castle, Martin van Rongen"],"number-sections":false},"headingText":"Core aims","containsRefs":false,"markdown":"\n\nWelcome to the wonderful world of generalised linear models!\n\nThese sessions are intended to enable you to construct and use generalised linear models confidently.\n\nAs with all of our statistics courses our focus is not on mathematical derivations, but on developing an intuitive understanding of the underlying statistical concepts.\n\nAt the same time this is also *not* a \"how to mindlessly use a stats program\" course! We hope that at the end of this course you feel like you have a better grasp on what it is we're trying to do, and gained sufficient confidence in your coding skills to implement these statistical concepts in your own research!\n\n\nTo introduce sufficient understanding and coding experience for analysing data with non-continuous response variables.\n\n::: callout-note\n## Course aims\n\nTo know what to do when presented with an arbitrary data set e.g.\n\n1. Construct\n a. a logistic model for binary response variables\n b. a logistic model for proportion response variables\n c. a Poisson model for count response variables\n d. ~~a Negative Binomial model for count response variables~~ (to be added later)\n2. Plot the data and the fitted curve in each case for both continuous and categorical predictors\n3. Assess the significance of fit\n4. Assess assumption of the model\n:::\n\n## Authors\n\nAbout the author(s):\n\n- **Vicki Hodgson**\n \n \n _Affiliation_: Bioinformatics Training Facility, University of Cambridge \n _Roles_: writing - review & editing; conceptualisation; coding\n- **Martin van Rongen**\n \n \n _Affiliation_: Bioinformatics Training Facility, University of Cambridge \n _Roles_: writing - review & editing; conceptualisation; coding\n- **Matt Castle**\n \n _Affiliation_: Bioinformatics Training Facility, University of Cambridge \n _Roles_: conceptualisation; writing\n\n","srcMarkdownNoYaml":"\n\nWelcome to the wonderful world of generalised linear models!\n\nThese sessions are intended to enable you to construct and use generalised linear models confidently.\n\nAs with all of our statistics courses our focus is not on mathematical derivations, but on developing an intuitive understanding of the underlying statistical concepts.\n\nAt the same time this is also *not* a \"how to mindlessly use a stats program\" course! We hope that at the end of this course you feel like you have a better grasp on what it is we're trying to do, and gained sufficient confidence in your coding skills to implement these statistical concepts in your own research!\n\n## Core aims\n\nTo introduce sufficient understanding and coding experience for analysing data with non-continuous response variables.\n\n::: callout-note\n## Course aims\n\nTo know what to do when presented with an arbitrary data set e.g.\n\n1. Construct\n a. a logistic model for binary response variables\n b. a logistic model for proportion response variables\n c. a Poisson model for count response variables\n d. ~~a Negative Binomial model for count response variables~~ (to be added later)\n2. Plot the data and the fitted curve in each case for both continuous and categorical predictors\n3. Assess the significance of fit\n4. Assess assumption of the model\n:::\n\n## Authors\n\nAbout the author(s):\n\n- **Vicki Hodgson**\n \n \n _Affiliation_: Bioinformatics Training Facility, University of Cambridge \n _Roles_: writing - review & editing; conceptualisation; coding\n- **Martin van Rongen**\n \n \n _Affiliation_: Bioinformatics Training Facility, University of Cambridge \n _Roles_: writing - review & editing; conceptualisation; coding\n- **Matt Castle**\n \n _Affiliation_: Bioinformatics Training Facility, University of Cambridge \n _Roles_: conceptualisation; writing\n\n"},"formats":{"courseformat-html":{"identifier":{"display-name":"HTML","target-format":"courseformat-html","base-format":"html","extension-name":"courseformat"},"execute":{"fig-width":7,"fig-height":5,"fig-format":"retina","fig-dpi":96,"df-print":"default","error":false,"eval":true,"cache":null,"freeze":"auto","echo":true,"output":true,"warning":true,"include":true,"keep-md":false,"keep-ipynb":false,"ipynb":null,"enabled":null,"daemon":null,"daemon-restart":false,"debug":false,"ipynb-filters":[],"ipynb-shell-interactivity":null,"plotly-connected":true,"engine":"markdown"},"render":{"keep-tex":false,"keep-typ":false,"keep-source":false,"keep-hidden":false,"prefer-html":false,"output-divs":true,"output-ext":"html","fig-align":"default","fig-pos":null,"fig-env":null,"code-fold":"none","code-overflow":"scroll","code-link":true,"code-line-numbers":false,"code-tools":false,"tbl-colwidths":"auto","merge-includes":true,"inline-includes":false,"preserve-yaml":false,"latex-auto-mk":true,"latex-auto-install":true,"latex-clean":true,"latex-min-runs":1,"latex-max-runs":10,"latex-makeindex":"makeindex","latex-makeindex-opts":[],"latex-tlmgr-opts":[],"latex-input-paths":[],"latex-output-dir":null,"link-external-icon":false,"link-external-newwindow":false,"self-contained-math":false,"format-resources":[],"notebook-links":true,"shortcodes":[]},"pandoc":{"standalone":true,"wrap":"none","default-image-extension":"png","to":"html","toc":true,"number-sections":false,"filters":["courseformat"],"output-file":"index.html"},"language":{"toc-title-document":"Table of contents","toc-title-website":"On this page","related-formats-title":"Other Formats","related-notebooks-title":"Notebooks","source-notebooks-prefix":"Source","other-links-title":"Other Links","code-links-title":"Code Links","launch-dev-container-title":"Launch Dev Container","launch-binder-title":"Launch Binder","article-notebook-label":"Article Notebook","notebook-preview-download":"Download Notebook","notebook-preview-download-src":"Download Source","notebook-preview-back":"Back to Article","manuscript-meca-bundle":"MECA Bundle","section-title-abstract":"Abstract","section-title-appendices":"Appendices","section-title-footnotes":"Footnotes","section-title-references":"References","section-title-reuse":"Reuse","section-title-copyright":"Copyright","section-title-citation":"Citation","appendix-attribution-cite-as":"For attribution, please cite this work as:","appendix-attribution-bibtex":"BibTeX citation:","title-block-author-single":"Author","title-block-author-plural":"Authors","title-block-affiliation-single":"Affiliation","title-block-affiliation-plural":"Affiliations","title-block-published":"Published","title-block-modified":"Modified","title-block-keywords":"Keywords","callout-tip-title":"Tip","callout-note-title":"Note","callout-warning-title":"Warning","callout-important-title":"Important","callout-caution-title":"Caution","code-summary":"Code","code-tools-menu-caption":"Code","code-tools-show-all-code":"Show All Code","code-tools-hide-all-code":"Hide All Code","code-tools-view-source":"View Source","code-tools-source-code":"Source Code","tools-share":"Share","tools-download":"Download","code-line":"Line","code-lines":"Lines","copy-button-tooltip":"Copy to Clipboard","copy-button-tooltip-success":"Copied!","repo-action-links-edit":"Edit this page","repo-action-links-source":"View source","repo-action-links-issue":"Report an issue","back-to-top":"Back to top","search-no-results-text":"No results","search-matching-documents-text":"matching documents","search-copy-link-title":"Copy link to search","search-hide-matches-text":"Hide additional matches","search-more-match-text":"more match in this document","search-more-matches-text":"more matches in this document","search-clear-button-title":"Clear","search-text-placeholder":"","search-detached-cancel-button-title":"Cancel","search-submit-button-title":"Submit","search-label":"Search","toggle-section":"Toggle section","toggle-sidebar":"Toggle sidebar navigation","toggle-dark-mode":"Toggle dark mode","toggle-reader-mode":"Toggle reader mode","toggle-navigation":"Toggle navigation","crossref-fig-title":"Figure","crossref-tbl-title":"Table","crossref-lst-title":"Listing","crossref-thm-title":"Theorem","crossref-lem-title":"Lemma","crossref-cor-title":"Corollary","crossref-prp-title":"Proposition","crossref-cnj-title":"Conjecture","crossref-def-title":"Definition","crossref-exm-title":"Example","crossref-exr-title":"Exercise","crossref-ch-prefix":"Chapter","crossref-apx-prefix":"Appendix","crossref-sec-prefix":"Section","crossref-eq-prefix":"Equation","crossref-lof-title":"List of Figures","crossref-lot-title":"List of Tables","crossref-lol-title":"List of Listings","environment-proof-title":"Proof","environment-remark-title":"Remark","environment-solution-title":"Solution","listing-page-order-by":"Order By","listing-page-order-by-default":"Default","listing-page-order-by-date-asc":"Oldest","listing-page-order-by-date-desc":"Newest","listing-page-order-by-number-desc":"High to Low","listing-page-order-by-number-asc":"Low to High","listing-page-field-date":"Date","listing-page-field-title":"Title","listing-page-field-description":"Description","listing-page-field-author":"Author","listing-page-field-filename":"File Name","listing-page-field-filemodified":"Modified","listing-page-field-subtitle":"Subtitle","listing-page-field-readingtime":"Reading Time","listing-page-field-wordcount":"Word Count","listing-page-field-categories":"Categories","listing-page-minutes-compact":"{0} min","listing-page-category-all":"All","listing-page-no-matches":"No matching items","listing-page-words":"{0} words"},"metadata":{"lang":"en","fig-responsive":true,"quarto-version":"1.4.546","theme":["default","_extensions/cambiotraining/courseformat/theme.scss"],"number-depth":3,"code-copy":true,"revealjs-plugins":[],"bibliography":["references.bib"],"knitr":{"opts_knit":{"cache.path":".knitr_cache"}},"title":"Course overview","author":["Vicki Hodgson, Matt Castle, Martin van Rongen"]},"extensions":{"book":{"multiFile":true}}}},"projectFormats":["courseformat-html"]} \ No newline at end of file diff --git a/_freeze/materials/checking-assumptions/execute-results/html.json b/_freeze/materials/checking-assumptions/execute-results/html.json new file mode 100644 index 0000000..c947048 --- /dev/null +++ b/_freeze/materials/checking-assumptions/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "e60a019a08a121704741c223e9cc90c0", + "result": { + "engine": "knitr", + "markdown": "---\ntitle: \"Checking assumptions\"\noutput: html_document\n---\n\n\n\n::: {.cell}\n\n:::\n\n::: {.cell}\n\n:::\n\n\nAlthough generalised linear models do allow us to relax certain assumptions compared to standard linear models (linearity, equality of variance of residuals, and normality of residuals).\n\nHowever, we cannot relax all of them. This section of the materials will talk through the important assumptions for GLMs, and how to assess them.\n\n## Libraries and functions\n\n::: {.callout-note collapse=\"true\"}\n## Click to expand\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\nlibrary(ggResidpanel)\n```\n:::\n\n## Python\n\n::: {.cell}\n\n```{.python .cell-code}\nfrom scipy.stats import *\n```\n:::\n\n:::\n\n:::\n\n## Assumption 1: Distribution of response variable\n\nAlthough we don't expect our response variable $y$ to be continuous and normally distributed (as we did in linear modelling), we do still expect its distribution to come from the \"exponential family\" of distributions.\n\nThe exponential family contains the following distributions, among others:\n\n- normal\n- exponential\n- Poisson \n- Bernoulli\n- binomial (for fixed number of trials)\n- chi-squared\n\nYou can use a histogram to visualise the distribution of your response variable, but it is typically most useful just to think about the nature of your response variable. For instance, binary variables will follow a Bernoulli distribution, proportional variables follow a binomial distribution, and most count variables will follow a Poisson distribution.\n\nIf you have a very unusual variable that doesn't follow one of these exponential family distributions, however, then a GLM will not be an appropriate choice. In other words, a GLM is not necessarily a magic fix!\n\n## Assumption 2: Correct link function\n\nA closely-related assumption to assumption 1 above, is that we have chosen the correct link function for our model.\n\nIf we have done so, then there should be a linear relationship between our *transformed* model and our response variable; in other words, if we have chosen the right link function, then we have correctly \"linearised\" our model.\n\n## Assumption 3: Independence\n\nWe expect that the each observation or datapoint in our sample is independent of all the others. Specifically, we expect that our set of $y$ response variables are independent of one another.\n\nFor this to be true, we have to make sure:\n\n- that we aren't treating technical replicates as true/biological replicates;\n- that we don't have observations/datapoints in our sample that are artificially similar to each other (compared to other datapoints);\n- that we don't have any nuisance/confounding variables that create \"clusters\" or hierarchy in our dataset;\n- that we haven't got repeated measures, i.e., multiple measurements/rows per individual in our sample\n\nThere is no diagnostic plot for assessing this assumption. To determine whether your data are independent, you need to understand your experimental design.\n\nYou might find [this page](https://cambiotraining.github.io/experimental-design/materials/04-replication.html#criteria-for-true-independent-replication) useful if you're looking for more information on what counts as truly independent data.\n\n## Good science: No influential observations\n\nAs with linear models, though this isn't always considered a \"formal\" assumption, we do want to ensure that there aren't any datapoints that are overly influencing our model.\n\nA datapoint is overly influential, i.e., has high leverage, if removing that point from the dataset would cause large changes in the model coefficients. Datapoints with high leverage are typically those that don't follow the same general \"trend\" as the rest of the data.\n\nThe easiest way to check for overly influential points is to construct a Cook's distance plot.\n\nLet's try that out, using the `diabetes` example dataset.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndiabetes <- read_csv(\"data/diabetes.csv\")\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\nRows: 728 Columns: 3\n── Column specification ────────────────────────────────────────────────────────\nDelimiter: \",\"\ndbl (3): glucose, diastolic, test_result\n\nℹ Use `spec()` to retrieve the full column specification for this data.\nℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.\n```\n\n\n:::\n\n```{.r .cell-code}\nglm_dia <- glm(test_result ~ glucose * diastolic,\n family = \"binomial\",\n data = diabetes)\n```\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\ndiabetes_py = pd.read_csv(\"data/diabetes.csv\")\n\nmodel = smf.glm(formula = \"test_result ~ glucose * diastolic\", \n family = sm.families.Binomial(), \n data = diabetes_py)\n \nglm_dia_py = model.fit()\n```\n:::\n\n:::\n\nOnce our model is fitted, we can fit a Cook's distance plot:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\nresid_panel(glm_dia, plots = \"cookd\")\n```\n\n::: {.cell-output-display}\n![](checking-assumptions_files/figure-html/unnamed-chunk-7-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n::: {.cell}\n\n:::\n\n:::\n\nGood news - there don't appear to be any overly influential points!\n\n## Dispersion\n\nAnother thing that we want to check, primarily in Poisson regression, is whether our dispersion parameter is correct.\n\n::: {.callout-note collapse=\"true}\n\n#### First, let's unpack what dispersion is!\n\nDispersion, in statistics, is a general term to describe the variability, scatter, or spread of a distribution. Variance is a common measure of dispersion that hopefully you are familiar with.\n\nIn a normal distribution, the mean (average) and the variance (dispersion) are independent of each other; we need both numbers, or parameters, to understand the shape of the distribution.\n\nOther distributions, however, require different parameters to describe them in full. For a Poisson distribution, we need just one parameter $lambda$, which captures the expected rate of occurrences/expected count. The mean and variance of a Poisson distribution are actually expected to be the same.\n\nIn the context of a model, you can think about the dispersion as the degree to which the data are spread out around the model curve. A dispersion parameter of 1 means the data are spread out exactly as we expect; <1 is called underdispersion; and >1 is called overdispersion.\n:::\n\n### A \"hidden assumption\"\n\nWhen we fit a linear model, because we're assuming a normal distribution, we take the time to estimate the dispersion - by measuring the variance. \n\nWhen performing Poisson regression, however, we make an extra \"hidden\" assumption, in setting the dispersion parameter to 1. In other words, we expect the errors to have a certain spread to them that matches our theoretical distribution/model. This means we don't have to waste time and statistical power in estimating the dispersion.\n\nHowever, if our data are underdispersed or overdispersed, then we might be violating this assumption we've made. \n\nUnderdispersion is quite rare. It's far more likely that you'll encounter overdispersion; in Poisson regression, this is usually caused by the presence of lots of zeroes in your response variable (known as zero-inflation).\n\nIn these situations, you may wish to fit a different GLM to the data. Negative binomial regression, for instance, is a common alternative for zero-inflated count data.\n\n### Checking the dispersion parameter\n\nThe easiest way to check dispersion in a model is to calculate the ratio of the residual deviance to the residual degrees of freedom.\n\nLet's practice doing this using a Poisson regression fitted to the `islands` dataset that you saw earlier in the course.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\nislands <- read_csv(\"data/islands.csv\")\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\nRows: 35 Columns: 2\n── Column specification ────────────────────────────────────────────────────────\nDelimiter: \",\"\ndbl (2): species, area\n\nℹ Use `spec()` to retrieve the full column specification for this data.\nℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.\n```\n\n\n:::\n\n```{.r .cell-code}\nglm_isl <- glm(species ~ area,\n data = islands, family = \"poisson\")\n```\n:::\n\n\n## Python\n\n::: {.cell}\n\n```{.python .cell-code}\nislands_py = pd.read_csv(\"data/islands.csv\")\n\nmodel = smf.glm(formula = \"species ~ area\",\n family = sm.families.Poisson(),\n data = islands_py)\n\nglm_isl_py = model.fit()\n```\n:::\n\n:::\n\nIf we take a look at the model output, we can see the two quantities we care about - residual deviance and residual degrees of freedom:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_isl)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = species ~ area, family = \"poisson\", data = islands)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) 4.241129 0.041322 102.64 <2e-16 ***\narea 0.035613 0.001247 28.55 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for poisson family taken to be 1)\n\n Null deviance: 856.899 on 34 degrees of freedom\nResidual deviance: 30.437 on 33 degrees of freedom\nAIC: 282.66\n\nNumber of Fisher Scoring iterations: 3\n```\n\n\n:::\n:::\n\n\n## Python\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_isl_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: species No. Observations: 35\nModel: GLM Df Residuals: 33\nModel Family: Poisson Df Model: 1\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -139.33\nDate: Fri, 17 May 2024 Deviance: 30.437\nTime: 12:40:44 Pearson chi2: 30.3\nNo. Iterations: 4 Pseudo R-squ. (CS): 1.000\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 4.2411 0.041 102.636 0.000 4.160 4.322\narea 0.0356 0.001 28.551 0.000 0.033 0.038\n==============================================================================\n```\n\n\n:::\n:::\n\n:::\n\nThe residual deviance is 30.437, on 33 residual degrees of freedom. All we need to do is divide one by the other to get our dispersion parameter.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_isl$deviance/glm_isl$df.residual\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0.922334\n```\n\n\n:::\n:::\n\n\n## Python\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_isl_py.deviance/glm_isl_py.df_resid)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n0.9223340414458532\n```\n\n\n:::\n:::\n\n:::\n\nThe dispersion parameter here is 0.922. That's pretty good - not far off 1 at all.\n\nBut how can we check whether it is *significantly* different from 1? \n\nWell, you've actually already got the knowledge you need to do this, from the previous course section on significance testing. Specifically, the chi-squared goodness-of-fit test can be used to check whether the dispersion is within sensible limits.\n\nYou may have noticed that the two values we're using for the dispersion parameter are the same two numbers that we used in those chi-squared tests. For this Poisson regression fitted to the `islands` dataset, that goodness-of-fit test would look like this:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(glm_isl$deviance, glm_isl$df.residual)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0.595347\n```\n\n\n:::\n:::\n\n\n## Python\n\n::: {.cell}\n\n```{.python .cell-code}\npvalue = chi2.sf(glm_isl_py.deviance, glm_isl_py.df_resid)\n\nprint(pvalue)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n0.5953470127463187\n```\n\n\n:::\n:::\n\n:::\n\nIf our chi-squared goodness-of-fit test returns a large (insignificant) p-value, as it does here, that tells us that we don't need to worry about the dispersion.\n\nIf our chi-squared goodness-of-fit test returned a small, significant p-value, this would tell us our model doesn't fit the data well. And, since dispersion is all about the spread of points around the model, it makes sense that these two things are so closely related!\n\n## Summary\n\nWhile generalised linear models make fewer assumptions than standard linear models, we do still expect certain things to be true about the model and our variables for GLMs to be valid. Checking most of these assumptions requires understanding your dataset, and diagnostic plots play a less heavy role.\n\n::: {.callout-tip}\n#### Key points\n\n- For a generalised linear model, we assume that we have chosen the correct link function, that our response variable follows a distribution from the exponential family, and that our data are independent\n- To assess these assumptions, we need to understand our dataset and variables\n- We can also use visualisation to determine whether we have overly influential (high leverage) datapoints\n- For Poisson regression, we should also investigate the dispersion parameter of our model, which we expect to be close to 1\n:::\n", + "supporting": [ + "checking-assumptions_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_freeze/materials/checking-assumptions/figure-html/unnamed-chunk-7-1.png b/_freeze/materials/checking-assumptions/figure-html/unnamed-chunk-7-1.png new file mode 100644 index 0000000..b1a0a52 Binary files /dev/null and b/_freeze/materials/checking-assumptions/figure-html/unnamed-chunk-7-1.png differ diff --git a/_freeze/materials/glm-practical-poisson/execute-results/html.json b/_freeze/materials/glm-practical-poisson/execute-results/html.json index 5a4ad40..d6865f9 100644 --- a/_freeze/materials/glm-practical-poisson/execute-results/html.json +++ b/_freeze/materials/glm-practical-poisson/execute-results/html.json @@ -1,8 +1,8 @@ { - "hash": "a7944edaa2f705428f83c869d555dfd3", + "hash": "bb8efe99ebfffabd6145784d8da531a1", "result": { "engine": "knitr", - "markdown": "---\ntitle: \"Count data\"\n---\n\n::: {.cell}\n\n:::\n\n::: {.cell}\n\n:::\n\n\n::: {.callout-tip}\n## Learning outcomes\n\n**Questions**\n\n- How do we analyse count data?\n\n**Objectives**\n\n- Be able to perform a poisson regression on count data\n:::\n\n## Libraries and functions\n\n::: {.callout-note collapse=\"true\"}\n## Click to expand\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n### Libraries\n### Functions\n\n## Python\n\n### Libraries\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# A maths library\nimport math\n# A Python data analysis and manipulation tool\nimport pandas as pd\n\n# Python equivalent of `ggplot2`\nfrom plotnine import *\n\n# Statistical models, conducting tests and statistical data exploration\nimport statsmodels.api as sm\n\n# Convenience interface for specifying models using formula strings and DataFrames\nimport statsmodels.formula.api as smf\n\n# Needed for additional probability functionality\nfrom scipy.stats import *\n```\n:::\n\n\n### Functions\n\n:::\n:::\n\nThe examples in this section use the following data sets:\n\n`data/islands.csv`\n\nThis is a data set comprising 35 observations of two variables (one dependent and one predictor). This records the number of species recorded on different small islands along with the area (km2) of the islands. The variables are `species` and `area`.\n\nThe second data set is on seat belts.\n\nThe `seatbelts` data set is a multiple time-series data set that was commissioned by the Department of Transport in 1984 to measure differences in deaths before and after front seat belt legislation was introduced on 31st January 1983. It provides monthly total numerical data on a number of incidents including those related to death and injury in Road Traffic Accidents (RTA's). The data set starts in January 1969 and observations run until December 1984.\n\nYou can find the file in `data/seatbelts.csv`\n\n## Load and visualise the data\n\nFirst we load the data, then we visualise it.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nislands <- read_csv(\"data/islands.csv\")\n```\n:::\n\n\nLet's have a glimpse at the data:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nislands\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 35 × 2\n species area\n \n 1 114 12.1\n 2 130 13.4\n 3 113 13.7\n 4 109 14.5\n 5 118 16.8\n 6 136 19.0\n 7 149 19.6\n 8 162 20.6\n 9 145 20.9\n10 148 21.0\n# ℹ 25 more rows\n```\n\n\n:::\n:::\n\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nislands_py = pd.read_csv(\"data/islands.csv\")\n```\n:::\n\n\nLet's have a glimpse at the data:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nislands_py.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n species area\n0 114 12.076133\n1 130 13.405439\n2 113 13.723525\n3 109 14.540359\n4 118 16.792122\n```\n\n\n:::\n:::\n\n\n:::\n\nLooking at the data, we can see that there are two columns: `species`, which contains the number of species recorded on each island and `area`, which contains the surface area of the island in square kilometers.\n\nWe can plot the data:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(islands, aes(x = area, y = species)) +\n geom_point()\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-8-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(islands_py, aes(x = \"area\", y = \"species\")) +\n geom_point())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-9-1.png){width=614}\n:::\n:::\n\n:::\n\nIt looks as though `area` may have an effect on the number of species that we observe on each island. We note that the response variable is count data and so we try to construct a Poisson regression.\n\n## Constructing a model\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_isl <- glm(species ~ area,\n data = islands, family = \"poisson\")\n```\n:::\n\n\nand we look at the model summary:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_isl)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = species ~ area, family = \"poisson\", data = islands)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) 4.241129 0.041322 102.64 <2e-16 ***\narea 0.035613 0.001247 28.55 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for poisson family taken to be 1)\n\n Null deviance: 856.899 on 34 degrees of freedom\nResidual deviance: 30.437 on 33 degrees of freedom\nAIC: 282.66\n\nNumber of Fisher Scoring iterations: 3\n```\n\n\n:::\n:::\n\n\nThe output is strikingly similar to the logistic regression models (who’d have guessed, eh?) and the main numbers to extract from the output are the two numbers underneath `Estimate.Std` in the `Coefficients` table:\n\n```\n(Intercept) 4.241129\narea 0.035613\n```\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# create a generalised linear model\nmodel = smf.glm(formula = \"species ~ area\",\n family = sm.families.Poisson(),\n data = islands_py)\n# and get the fitted parameters of the model\nglm_isl_py = model.fit()\n```\n:::\n\n\nLet's look at the model output:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_isl_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: species No. Observations: 35\nModel: GLM Df Residuals: 33\nModel Family: Poisson Df Model: 1\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -139.33\nDate: Tue, 06 Feb 2024 Deviance: 30.437\nTime: 16:16:33 Pearson chi2: 30.3\nNo. Iterations: 4 Pseudo R-squ. (CS): 1.000\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 4.2411 0.041 102.636 0.000 4.160 4.322\narea 0.0356 0.001 28.551 0.000 0.033 0.038\n==============================================================================\n```\n\n\n:::\n:::\n\n\n:::\n\nThese are the coefficients of the Poisson model equation and need to be placed in the following formula in order to estimate the expected number of species as a function of island size:\n\n$$ E(species) = \\exp(4.24 + 0.036 \\times area) $$\n\nInterpreting this requires a bit of thought (not much, but a bit).\nThe intercept coefficient, `4.24`, is related to the number of species we would expect on an island of zero area (this is statistics, not real life. You’d do well to remember that before you worry too much about what that even means). But in order to turn this number into something meaningful we have to exponentiate it. Since `exp(4.24) ≈ 70`, we can say that the baseline number of species the model expects on any island is 70. This isn’t actually the interesting bit though.\n\nThe coefficient of `area` is the fun bit. For starters we can see that it is a positive number which does mean that increasing `area` leads to increasing numbers of `species`. Good so far.\n\nBut what does the value `0.036` actually mean? Well, if we exponentiate it as well, we get `exp(0.036) ≈ 1.04`. This means that for every increase in `area` of 1 km^2 (the original units of the area variable), the number of species on the island is multiplied by `1.04`. So, an island of area 1 km^2 will have `1.04 x 70 ≈ 72` species.\n\nSo, in order to interpret Poisson coefficients, you have to exponentiate them.\n\n## Plotting the Poisson regression\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(islands, aes(area, species)) +\n geom_point() +\n geom_smooth(method = \"glm\", se = FALSE, fullrange = TRUE, \n method.args = list(family = poisson)) +\n xlim(10,50)\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\n`geom_smooth()` using formula = 'y ~ x'\n```\n\n\n:::\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-14-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = pd.DataFrame({'area': list(range(10, 50))})\n\nmodel[\"pred\"] = glm_isl_py.predict(model)\n\nmodel.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n area pred\n0 10 99.212463\n1 11 102.809432\n2 12 106.536811\n3 13 110.399326\n4 14 114.401877\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(islands_py,\n aes(x = \"area\",\n y = \"species\")) +\n geom_point() +\n geom_line(model, aes(x = \"area\", y = \"pred\"), colour = \"blue\", size = 1))\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-16-1.png){width=614}\n:::\n:::\n\n:::\n\n## Assumptions\n\nAs we mentioned earlier, Poisson regressions require that the variance of the data at any point is the same as the mean of the data at that point. We checked that earlier by looking at the residual deviance values.\n\nWe can look for influential points using the Cook’s distance plot:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nplot(glm_isl , which=4)\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-17-3.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# extract the Cook's distances\nglm_isl_py_resid = pd.DataFrame(glm_isl_py.\n get_influence().\n summary_frame()[\"cooks_d\"])\n\n# add row index \nglm_isl_py_resid['obs'] = glm_isl_py_resid.reset_index().index\n```\n:::\n\n\nWe can use these to create the plot:\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(glm_isl_py_resid,\n aes(x = \"obs\",\n y = \"cooks_d\")) +\n geom_segment(aes(x = \"obs\", y = \"cooks_d\", xend = \"obs\", yend = 0)) +\n geom_point())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-19-1.png){width=614}\n:::\n:::\n\n\n:::\n\nNone of our points have particularly large Cook’s distances and so life is rosy.\n\n## Assessing significance\n\nWe can ask the same three questions we asked before.\n\n1. Is the model well-specified?\n2. Is the overall model better than the null model?\n3. Are any of the individual predictors significant?\n\nAgain, in this case, questions 2 and 3 are effectively asking the same thing because we still only have a single predictor variable.\n\nTo assess if the model is any good we’ll again use the residual deviance and the residual degrees of freedom.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(30.437, 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0.5953482\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(30.437, 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n0.5953481872979622\n```\n\n\n:::\n:::\n\n\n:::\n\nThis gives a probability of `0.60`. This suggests that this model is actually a reasonably decent one and that the data are pretty well supported by the model. For Poisson models this has an extra interpretation. This can be used to assess whether we have significant over-dispersion in our data.\n\nFor a Poisson model to be appropriate we need that the variance of the data to be exactly the same as the mean of the data. Visually, this would correspond to the data spreading out more for higher predicted values of `species.` However, we don’t want the data to spread out too much. If that happens then a Poisson model wouldn’t be appropriate.\n\nThe easy way to check this is to look at the ratio of the residual deviance to the residual degrees of freedom (in this case `0.922`). For a Poisson model to be valid, this ratio should be about 1. If the ratio is significantly bigger than 1 then we say that we have over-dispersion in the model and we wouldn’t be able to trust any of the significance testing that we are about to do using a Poisson regression.\n\nThankfully the probability we have just created (`0.60`) is exactly the right one we need to look at to assess whether we have significant over-dispersion in our model.\n\nSecondly, to assess whether the overall model, with all of the terms, is better than the null model we’ll look at the difference in deviances and the difference in degrees of freedom:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(856.899 - 30.437, 34 - 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(856.899 - 30.437, 34 - 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n9.524927068555617e-182\n```\n\n\n:::\n:::\n\n:::\n\nThis gives a reported p-value of pretty much zero, which is pretty damn small. So, yes, this model is better than nothing at all and species does appear to change with some of our predictors\n\nFinally, we’ll construct an analysis of deviance table to look at the individual terms:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nanova(glm_isl , test = \"Chisq\")\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nAnalysis of Deviance Table\n\nModel: poisson, link: log\n\nResponse: species\n\nTerms added sequentially (first to last)\n\n Df Deviance Resid. Df Resid. Dev Pr(>Chi) \nNULL 34 856.90 \narea 1 826.46 33 30.44 < 2.2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n```\n\n\n:::\n:::\n\n\nThe p-value in this table is just as small as we’d expect given our previous result (`<2.2e-16` is pretty close to 0), and we have the nice consistent result that `area` definitely has an effect on `species`.\n\n## Python\n\nAs mentioned before, this is not quite possible in Python.\n:::\n\n## Exercises\n\n### Seat belts {#sec-exr_seatbelts}\n\n:::{.callout-exercise}\n\n\n{{< level 2 >}}\n\n\n\nFor this exercise we'll be using the data from `data/seatbelts.csv`.\n\nI'd like you to do the following:\n\n1. Load the data\n2. Visualise the data and create a poisson regression model\n3. Plot the regression model on top of the data\n4. Assess if the model is a decent predictor for the number of fatalities\n\n::: {.callout-answer collapse=\"true\"}\n\n#### Load and visualise the data\n\nFirst we load the data, then we visualise it.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nseatbelts <- read_csv(\"data/seatbelts.csv\")\n```\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nseatbelts_py = pd.read_csv(\"data/seatbelts.csv\")\n```\n:::\n\n\nLet's have a glimpse at the data:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nseatbelts_py.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n casualties drivers front rear ... van_killed law year month\n0 107 1687 867 269 ... 12 0 1969 Jan\n1 97 1508 825 265 ... 6 0 1969 Feb\n2 102 1507 806 319 ... 12 0 1969 Mar\n3 87 1385 814 407 ... 8 0 1969 Apr\n4 119 1632 991 454 ... 10 0 1969 May\n\n[5 rows x 10 columns]\n```\n\n\n:::\n:::\n\n:::\n\nThe data tracks the number of drivers killed in road traffic accidents, before and after the seat belt law was introduced. The information on whether the law was in place is encoded in the `law` column as `0` (law not in place) or `1` (law in place).\n\nThere are many more observations when the law was *not* in place, so we need to keep this in mind when we're interpreting the data.\n\nFirst we have a look at the data comparing no law vs law:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\nWe have to convert the `law` column to a factor, otherwise R will see it as numerical.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nseatbelts %>% \n ggplot(aes(as_factor(law), casualties)) +\n geom_boxplot()\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-28-1.png){width=672}\n:::\n:::\n\n\nThe data are recorded by month and year, so we can also display the number of drivers killed by year:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nseatbelts %>% \n ggplot(aes(year, casualties)) +\n geom_point()\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-29-1.png){width=672}\n:::\n:::\n\n\n## Python\n\nWe have to convert the `law` column to a factor, otherwise R will see it as numerical.\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(seatbelts_py,\n aes(x = seatbelts_py.law.astype(object),\n y = \"casualties\")) +\n geom_boxplot())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-30-1.png){width=614}\n:::\n:::\n\n\nThe data are recorded by month and year, so we can also display the number of casualties by year:\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(seatbelts_py,\n aes(x = \"year\",\n y = \"casualties\")) +\n geom_point())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-31-3.png){width=614}\n:::\n:::\n\n\n:::\n\nThe data look a bit weird. There is quite some variation within years (keeping in mind that the data are aggregated monthly). The data also seems to wave around a bit... with some vague peaks (e.g. 1972 - 1973) and some troughs (e.g. around 1976).\n\nSo my initial thought is that these data are going to be a bit tricky to interpret. But that's OK.\n\n#### Constructing a model\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_stb <- glm(casualties ~ year,\n data = seatbelts, family = \"poisson\")\n```\n:::\n\n\nand we look at the model summary:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_stb)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = casualties ~ year, family = \"poisson\", data = seatbelts)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) 37.168958 2.796636 13.29 <2e-16 ***\nyear -0.016373 0.001415 -11.57 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for poisson family taken to be 1)\n\n Null deviance: 984.50 on 191 degrees of freedom\nResidual deviance: 850.41 on 190 degrees of freedom\nAIC: 2127.2\n\nNumber of Fisher Scoring iterations: 4\n```\n\n\n:::\n:::\n\n\n```\n(Intercept) 37.168958\nyear -0.016373\n```\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# create a linear model\nmodel = smf.glm(formula = \"casualties ~ year\",\n family = sm.families.Poisson(),\n data = seatbelts_py)\n# and get the fitted parameters of the model\nglm_stb_py = model.fit()\n```\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_stb_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: casualties No. Observations: 192\nModel: GLM Df Residuals: 190\nModel Family: Poisson Df Model: 1\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -1061.6\nDate: Tue, 06 Feb 2024 Deviance: 850.41\nTime: 16:16:38 Pearson chi2: 862.\nNo. Iterations: 4 Pseudo R-squ. (CS): 0.5026\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 37.1690 2.797 13.291 0.000 31.688 42.650\nyear -0.0164 0.001 -11.569 0.000 -0.019 -0.014\n==============================================================================\n```\n\n\n:::\n:::\n\n\n```\n======================\n coef \n----------------------\nIntercept 37.1690 \nyear -0.0164 \n======================\n```\n:::\n\nThese are the coefficients of the Poisson model equation and need to be placed in the following formula in order to estimate the expected number of species as a function of island size:\n\n$$ E(casualties) = \\exp(37.17 - 0.164 \\times year) $$\n\n#### Assessing significance\n\nIs the model well-specified?\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(850.41, 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(850.41, 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n3.1319689119997022e-84\n```\n\n\n:::\n:::\n\n:::\n\nThis value indicates that the model is actually pretty good. Remember, it is between $[0, 1]$ and the closer to zero, the better the model.\n\nHow about the overall fit?\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(984.50 - 850.41, 191 - 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0\n```\n\n\n:::\n:::\n\n\n## Python\n\nFirst we need to define the null model:\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# create a linear model\nmodel = smf.glm(formula = \"casualties ~ 1\",\n family = sm.families.Poisson(),\n data = seatbelts_py)\n# and get the fitted parameters of the model\nglm_stb_null_py = model.fit()\n\nprint(glm_stb_null_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: casualties No. Observations: 192\nModel: GLM Df Residuals: 191\nModel Family: Poisson Df Model: 0\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -1128.6\nDate: Tue, 06 Feb 2024 Deviance: 984.50\nTime: 16:16:38 Pearson chi2: 1.00e+03\nNo. Iterations: 4 Pseudo R-squ. (CS): 1.942e-13\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 4.8106 0.007 738.670 0.000 4.798 4.823\n==============================================================================\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(984.50 - 850.41, 191 - 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n5.2214097202831414e-31\n```\n\n\n:::\n:::\n\n:::\n\nAgain, this indicates that the model is markedly better than the null model.\n\n#### Plotting the regression\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(seatbelts, aes(year, casualties)) +\n geom_point() +\n geom_smooth(method = \"glm\", se = FALSE, fullrange = TRUE, \n method.args = list(family = poisson)) +\n xlim(1970,1985)\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-41-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = pd.DataFrame({'year': list(range(1968, 1985))})\n\nmodel[\"pred\"] = glm_stb_py.predict(model)\n\nmodel.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n year pred\n0 1968 140.737690\n1 1969 138.452153\n2 1970 136.203733\n3 1971 133.991827\n4 1972 131.815842\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(seatbelts_py,\n aes(x = \"year\",\n y = \"casualties\")) +\n geom_point() +\n geom_line(model, aes(x = \"year\", y = \"pred\"), colour = \"blue\", size = 1))\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-43-1.png){width=614}\n:::\n:::\n\n:::\n\n\n#### Conclusions\n\nThe model we constructed appears to be a decent predictor for the number of fatalities.\n\n:::\n:::\n\n## Summary\n\n::: {.callout-tip}\n#### Key points\n\n- Poisson regression is useful when dealing with count data\n:::\n", + "markdown": "---\ntitle: \"Count data\"\n---\n\n::: {.cell}\n\n:::\n\n::: {.cell}\n\n:::\n\n\n::: {.callout-tip}\n## Learning outcomes\n\n**Questions**\n\n- How do we analyse count data?\n\n**Objectives**\n\n- Be able to perform a poisson regression on count data\n:::\n\n## Libraries and functions\n\n::: {.callout-note collapse=\"true\"}\n## Click to expand\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n### Libraries\n### Functions\n\n## Python\n\n### Libraries\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# A maths library\nimport math\n# A Python data analysis and manipulation tool\nimport pandas as pd\n\n# Python equivalent of `ggplot2`\nfrom plotnine import *\n\n# Statistical models, conducting tests and statistical data exploration\nimport statsmodels.api as sm\n\n# Convenience interface for specifying models using formula strings and DataFrames\nimport statsmodels.formula.api as smf\n\n# Needed for additional probability functionality\nfrom scipy.stats import *\n```\n:::\n\n\n### Functions\n\n:::\n:::\n\nThe examples in this section use the following data sets:\n\n`data/islands.csv`\n\nThis is a data set comprising 35 observations of two variables (one dependent and one predictor). This records the number of species recorded on different small islands along with the area (km2) of the islands. The variables are `species` and `area`.\n\nThe second data set is on seat belts.\n\nThe `seatbelts` data set is a multiple time-series data set that was commissioned by the Department of Transport in 1984 to measure differences in deaths before and after front seat belt legislation was introduced on 31st January 1983. It provides monthly total numerical data on a number of incidents including those related to death and injury in Road Traffic Accidents (RTA's). The data set starts in January 1969 and observations run until December 1984.\n\nYou can find the file in `data/seatbelts.csv`\n\n## Load and visualise the data\n\nFirst we load the data, then we visualise it.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nislands <- read_csv(\"data/islands.csv\")\n```\n:::\n\n\nLet's have a glimpse at the data:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nislands\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n# A tibble: 35 × 2\n species area\n \n 1 114 12.1\n 2 130 13.4\n 3 113 13.7\n 4 109 14.5\n 5 118 16.8\n 6 136 19.0\n 7 149 19.6\n 8 162 20.6\n 9 145 20.9\n10 148 21.0\n# ℹ 25 more rows\n```\n\n\n:::\n:::\n\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nislands_py = pd.read_csv(\"data/islands.csv\")\n```\n:::\n\n\nLet's have a glimpse at the data:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nislands_py.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n species area\n0 114 12.076133\n1 130 13.405439\n2 113 13.723525\n3 109 14.540359\n4 118 16.792122\n```\n\n\n:::\n:::\n\n\n:::\n\nLooking at the data, we can see that there are two columns: `species`, which contains the number of species recorded on each island and `area`, which contains the surface area of the island in square kilometers.\n\nWe can plot the data:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(islands, aes(x = area, y = species)) +\n geom_point()\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-8-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(islands_py, aes(x = \"area\", y = \"species\")) +\n geom_point())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-9-1.png){width=614}\n:::\n:::\n\n:::\n\nIt looks as though `area` may have an effect on the number of species that we observe on each island. We note that the response variable is count data and so we try to construct a Poisson regression.\n\n## Constructing a model\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_isl <- glm(species ~ area,\n data = islands, family = \"poisson\")\n```\n:::\n\n\nand we look at the model summary:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_isl)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = species ~ area, family = \"poisson\", data = islands)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) 4.241129 0.041322 102.64 <2e-16 ***\narea 0.035613 0.001247 28.55 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for poisson family taken to be 1)\n\n Null deviance: 856.899 on 34 degrees of freedom\nResidual deviance: 30.437 on 33 degrees of freedom\nAIC: 282.66\n\nNumber of Fisher Scoring iterations: 3\n```\n\n\n:::\n:::\n\n\nThe output is strikingly similar to the logistic regression models (who’d have guessed, eh?) and the main numbers to extract from the output are the two numbers underneath `Estimate.Std` in the `Coefficients` table:\n\n```\n(Intercept) 4.241129\narea 0.035613\n```\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# create a generalised linear model\nmodel = smf.glm(formula = \"species ~ area\",\n family = sm.families.Poisson(),\n data = islands_py)\n# and get the fitted parameters of the model\nglm_isl_py = model.fit()\n```\n:::\n\n\nLet's look at the model output:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_isl_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: species No. Observations: 35\nModel: GLM Df Residuals: 33\nModel Family: Poisson Df Model: 1\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -139.33\nDate: Fri, 17 May 2024 Deviance: 30.437\nTime: 12:30:01 Pearson chi2: 30.3\nNo. Iterations: 4 Pseudo R-squ. (CS): 1.000\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 4.2411 0.041 102.636 0.000 4.160 4.322\narea 0.0356 0.001 28.551 0.000 0.033 0.038\n==============================================================================\n```\n\n\n:::\n:::\n\n\n:::\n\nThese are the coefficients of the Poisson model equation and need to be placed in the following formula in order to estimate the expected number of species as a function of island size:\n\n$$ E(species) = \\exp(4.24 + 0.036 \\times area) $$\n\nInterpreting this requires a bit of thought (not much, but a bit).\nThe intercept coefficient, `4.24`, is related to the number of species we would expect on an island of zero area (this is statistics, not real life. You’d do well to remember that before you worry too much about what that even means). But in order to turn this number into something meaningful we have to exponentiate it. Since `exp(4.24) ≈ 70`, we can say that the baseline number of species the model expects on any island is 70. This isn’t actually the interesting bit though.\n\nThe coefficient of `area` is the fun bit. For starters we can see that it is a positive number which does mean that increasing `area` leads to increasing numbers of `species`. Good so far.\n\nBut what does the value `0.036` actually mean? Well, if we exponentiate it as well, we get `exp(0.036) ≈ 1.04`. This means that for every increase in `area` of 1 km^2 (the original units of the area variable), the number of species on the island is multiplied by `1.04`. So, an island of area 1 km^2 will have `1.04 x 70 ≈ 72` species.\n\nSo, in order to interpret Poisson coefficients, you have to exponentiate them.\n\n## Plotting the Poisson regression\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(islands, aes(area, species)) +\n geom_point() +\n geom_smooth(method = \"glm\", se = FALSE, fullrange = TRUE, \n method.args = list(family = poisson)) +\n xlim(10,50)\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\n`geom_smooth()` using formula = 'y ~ x'\n```\n\n\n:::\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-14-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = pd.DataFrame({'area': list(range(10, 50))})\n\nmodel[\"pred\"] = glm_isl_py.predict(model)\n\nmodel.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n area pred\n0 10 99.212463\n1 11 102.809432\n2 12 106.536811\n3 13 110.399326\n4 14 114.401877\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(islands_py,\n aes(x = \"area\",\n y = \"species\")) +\n geom_point() +\n geom_line(model, aes(x = \"area\", y = \"pred\"), colour = \"blue\", size = 1))\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-16-1.png){width=614}\n:::\n:::\n\n:::\n\n## Assumptions\n\nAs we mentioned earlier, Poisson regressions require that the variance of the data at any point is the same as the mean of the data at that point. We checked that earlier by looking at the residual deviance values.\n\nWe can look for influential points using the Cook’s distance plot:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nplot(glm_isl , which = 4)\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-17-3.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# extract the Cook's distances\nglm_isl_py_resid = pd.DataFrame(glm_isl_py.\n get_influence().\n summary_frame()[\"cooks_d\"])\n\n# add row index \nglm_isl_py_resid['obs'] = glm_isl_py_resid.reset_index().index\n```\n:::\n\n\nWe can use these to create the plot:\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(glm_isl_py_resid,\n aes(x = \"obs\",\n y = \"cooks_d\")) +\n geom_segment(aes(x = \"obs\", y = \"cooks_d\", xend = \"obs\", yend = 0)) +\n geom_point())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-19-1.png){width=614}\n:::\n:::\n\n\n:::\n\nNone of our points have particularly large Cook’s distances and so life is rosy.\n\n## Assessing significance\n\nWe can ask the same three questions we asked before.\n\n1. Is the model well-specified?\n2. Is the overall model better than the null model?\n3. Are any of the individual predictors significant?\n\nAgain, in this case, questions 2 and 3 are effectively asking the same thing because we still only have a single predictor variable.\n\nTo assess if the model is any good we’ll again use the residual deviance and the residual degrees of freedom.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(30.437, 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0.5953482\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(30.437, 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n0.5953481872979622\n```\n\n\n:::\n:::\n\n\n:::\n\nThis gives a probability of `0.60`. This suggests that this model is actually a reasonably decent one and that the data are pretty well supported by the model. For Poisson models this has an extra interpretation. This can be used to assess whether we have significant over-dispersion in our data.\n\nFor a Poisson model to be appropriate we need that the variance of the data to be exactly the same as the mean of the data. Visually, this would correspond to the data spreading out more for higher predicted values of `species.` However, we don’t want the data to spread out too much. If that happens then a Poisson model wouldn’t be appropriate.\n\nThe easy way to check this is to look at the ratio of the residual deviance to the residual degrees of freedom (in this case `0.922`). For a Poisson model to be valid, this ratio should be about 1. If the ratio is significantly bigger than 1 then we say that we have over-dispersion in the model and we wouldn’t be able to trust any of the significance testing that we are about to do using a Poisson regression.\n\nThankfully the probability we have just created (`0.60`) is exactly the right one we need to look at to assess whether we have significant over-dispersion in our model.\n\nSecondly, to assess whether the overall model, with all of the terms, is better than the null model we’ll look at the difference in deviances and the difference in degrees of freedom:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(856.899 - 30.437, 34 - 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(856.899 - 30.437, 34 - 33)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n9.524927068555617e-182\n```\n\n\n:::\n:::\n\n:::\n\nThis gives a reported p-value of pretty much zero, which is pretty damn small. So, yes, this model is better than nothing at all and species does appear to change with some of our predictors\n\nFinally, we’ll construct an analysis of deviance table to look at the individual terms:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nanova(glm_isl , test = \"Chisq\")\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nAnalysis of Deviance Table\n\nModel: poisson, link: log\n\nResponse: species\n\nTerms added sequentially (first to last)\n\n Df Deviance Resid. Df Resid. Dev Pr(>Chi) \nNULL 34 856.90 \narea 1 826.46 33 30.44 < 2.2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n```\n\n\n:::\n:::\n\n\nThe p-value in this table is just as small as we’d expect given our previous result (`<2.2e-16` is pretty close to 0), and we have the nice consistent result that `area` definitely has an effect on `species`.\n\n## Python\n\nAs mentioned before, this is not quite possible in Python.\n:::\n\n## Exercises\n\n### Seat belts {#sec-exr_seatbelts}\n\n:::{.callout-exercise}\n\n\n{{< level 2 >}}\n\n\n\nFor this exercise we'll be using the data from `data/seatbelts.csv`.\n\nI'd like you to do the following:\n\n1. Load the data\n2. Visualise the data and create a poisson regression model\n3. Plot the regression model on top of the data\n4. Assess if the model is a decent predictor for the number of fatalities\n\n::: {.callout-answer collapse=\"true\"}\n\n#### Load and visualise the data\n\nFirst we load the data, then we visualise it.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nseatbelts <- read_csv(\"data/seatbelts.csv\")\n```\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nseatbelts_py = pd.read_csv(\"data/seatbelts.csv\")\n```\n:::\n\n\nLet's have a glimpse at the data:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nseatbelts_py.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n casualties drivers front rear ... van_killed law year month\n0 107 1687 867 269 ... 12 0 1969 Jan\n1 97 1508 825 265 ... 6 0 1969 Feb\n2 102 1507 806 319 ... 12 0 1969 Mar\n3 87 1385 814 407 ... 8 0 1969 Apr\n4 119 1632 991 454 ... 10 0 1969 May\n\n[5 rows x 10 columns]\n```\n\n\n:::\n:::\n\n:::\n\nThe data tracks the number of drivers killed in road traffic accidents, before and after the seat belt law was introduced. The information on whether the law was in place is encoded in the `law` column as `0` (law not in place) or `1` (law in place).\n\nThere are many more observations when the law was *not* in place, so we need to keep this in mind when we're interpreting the data.\n\nFirst we have a look at the data comparing no law vs law:\n\n::: {.panel-tabset group=\"language\"}\n## R\n\nWe have to convert the `law` column to a factor, otherwise R will see it as numerical.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nseatbelts %>% \n ggplot(aes(as_factor(law), casualties)) +\n geom_boxplot()\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-28-1.png){width=672}\n:::\n:::\n\n\nThe data are recorded by month and year, so we can also display the number of drivers killed by year:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nseatbelts %>% \n ggplot(aes(year, casualties)) +\n geom_point()\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-29-1.png){width=672}\n:::\n:::\n\n\n## Python\n\nWe have to convert the `law` column to a factor, otherwise R will see it as numerical.\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(seatbelts_py,\n aes(x = seatbelts_py.law.astype(object),\n y = \"casualties\")) +\n geom_boxplot())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-30-1.png){width=614}\n:::\n:::\n\n\nThe data are recorded by month and year, so we can also display the number of casualties by year:\n\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(seatbelts_py,\n aes(x = \"year\",\n y = \"casualties\")) +\n geom_point())\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-31-3.png){width=614}\n:::\n:::\n\n\n:::\n\nThe data look a bit weird. There is quite some variation within years (keeping in mind that the data are aggregated monthly). The data also seems to wave around a bit... with some vague peaks (e.g. 1972 - 1973) and some troughs (e.g. around 1976).\n\nSo my initial thought is that these data are going to be a bit tricky to interpret. But that's OK.\n\n#### Constructing a model\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_stb <- glm(casualties ~ year,\n data = seatbelts, family = \"poisson\")\n```\n:::\n\n\nand we look at the model summary:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_stb)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = casualties ~ year, family = \"poisson\", data = seatbelts)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) 37.168958 2.796636 13.29 <2e-16 ***\nyear -0.016373 0.001415 -11.57 <2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for poisson family taken to be 1)\n\n Null deviance: 984.50 on 191 degrees of freedom\nResidual deviance: 850.41 on 190 degrees of freedom\nAIC: 2127.2\n\nNumber of Fisher Scoring iterations: 4\n```\n\n\n:::\n:::\n\n\n```\n(Intercept) 37.168958\nyear -0.016373\n```\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# create a linear model\nmodel = smf.glm(formula = \"casualties ~ year\",\n family = sm.families.Poisson(),\n data = seatbelts_py)\n# and get the fitted parameters of the model\nglm_stb_py = model.fit()\n```\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_stb_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: casualties No. Observations: 192\nModel: GLM Df Residuals: 190\nModel Family: Poisson Df Model: 1\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -1061.6\nDate: Fri, 17 May 2024 Deviance: 850.41\nTime: 12:30:06 Pearson chi2: 862.\nNo. Iterations: 4 Pseudo R-squ. (CS): 0.5026\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 37.1690 2.797 13.291 0.000 31.688 42.650\nyear -0.0164 0.001 -11.569 0.000 -0.019 -0.014\n==============================================================================\n```\n\n\n:::\n:::\n\n\n```\n======================\n coef \n----------------------\nIntercept 37.1690 \nyear -0.0164 \n======================\n```\n:::\n\nThese are the coefficients of the Poisson model equation and need to be placed in the following formula in order to estimate the expected number of species as a function of island size:\n\n$$ E(casualties) = \\exp(37.17 - 0.164 \\times year) $$\n\n#### Assessing significance\n\nIs the model well-specified?\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(850.41, 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(850.41, 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n3.1319689119997022e-84\n```\n\n\n:::\n:::\n\n:::\n\nThis value indicates that the model is actually pretty good. Remember, it is between $[0, 1]$ and the closer to zero, the better the model.\n\nHow about the overall fit?\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(984.50 - 850.41, 191 - 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0\n```\n\n\n:::\n:::\n\n\n## Python\n\nFirst we need to define the null model:\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# create a linear model\nmodel = smf.glm(formula = \"casualties ~ 1\",\n family = sm.families.Poisson(),\n data = seatbelts_py)\n# and get the fitted parameters of the model\nglm_stb_null_py = model.fit()\n\nprint(glm_stb_null_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: casualties No. Observations: 192\nModel: GLM Df Residuals: 191\nModel Family: Poisson Df Model: 0\nLink Function: Log Scale: 1.0000\nMethod: IRLS Log-Likelihood: -1128.6\nDate: Fri, 17 May 2024 Deviance: 984.50\nTime: 12:30:07 Pearson chi2: 1.00e+03\nNo. Iterations: 4 Pseudo R-squ. (CS): 1.942e-13\nCovariance Type: nonrobust \n==============================================================================\n coef std err z P>|z| [0.025 0.975]\n------------------------------------------------------------------------------\nIntercept 4.8106 0.007 738.670 0.000 4.798 4.823\n==============================================================================\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\nchi2.sf(984.50 - 850.41, 191 - 190)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n5.2214097202831414e-31\n```\n\n\n:::\n:::\n\n:::\n\nAgain, this indicates that the model is markedly better than the null model.\n\n#### Plotting the regression\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nggplot(seatbelts, aes(year, casualties)) +\n geom_point() +\n geom_smooth(method = \"glm\", se = FALSE, fullrange = TRUE, \n method.args = list(family = poisson)) +\n xlim(1970, 1985)\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-41-1.png){width=672}\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = pd.DataFrame({'year': list(range(1968, 1985))})\n\nmodel[\"pred\"] = glm_stb_py.predict(model)\n\nmodel.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n year pred\n0 1968 140.737690\n1 1969 138.452153\n2 1970 136.203733\n3 1971 133.991827\n4 1972 131.815842\n```\n\n\n:::\n:::\n\n::: {.cell}\n\n```{.python .cell-code}\n(ggplot(seatbelts_py,\n aes(x = \"year\",\n y = \"casualties\")) +\n geom_point() +\n geom_line(model, aes(x = \"year\", y = \"pred\"), colour = \"blue\", size = 1))\n```\n\n::: {.cell-output-display}\n![](glm-practical-poisson_files/figure-html/unnamed-chunk-43-1.png){width=614}\n:::\n:::\n\n:::\n\n\n#### Conclusions\n\nThe model we constructed appears to be a decent predictor for the number of fatalities.\n\n:::\n:::\n\n## Summary\n\n::: {.callout-tip}\n#### Key points\n\n- Poisson regression is useful when dealing with count data\n:::\n", "supporting": [ "glm-practical-poisson_files" ], diff --git a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-16-1.png b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-16-1.png index 1c0b16a..416b5d3 100644 Binary files a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-16-1.png and b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-16-1.png differ diff --git a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-19-1.png b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-19-1.png index 829445c..f466d73 100644 Binary files a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-19-1.png and b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-19-1.png differ diff --git a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-30-1.png b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-30-1.png index 6ea8e6f..437228b 100644 Binary files a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-30-1.png and b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-30-1.png differ diff --git a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-31-3.png b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-31-3.png index b7ce4a8..d55a9ec 100644 Binary files a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-31-3.png and b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-31-3.png differ diff --git a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-43-1.png b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-43-1.png index a1b6b49..52baf93 100644 Binary files a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-43-1.png and b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-43-1.png differ diff --git a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-9-1.png b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-9-1.png index 72005bc..c9bcee1 100644 Binary files a/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-9-1.png and b/_freeze/materials/glm-practical-poisson/figure-html/unnamed-chunk-9-1.png differ diff --git a/_freeze/materials/significance-testing/execute-results/html.json b/_freeze/materials/significance-testing/execute-results/html.json new file mode 100644 index 0000000..d77988d --- /dev/null +++ b/_freeze/materials/significance-testing/execute-results/html.json @@ -0,0 +1,17 @@ +{ + "hash": "9b1997b4e655bbd2efb84094264a9b0c", + "result": { + "engine": "knitr", + "markdown": "---\ntitle: \"Significance & goodness-of-fit\"\noutput: html_document\n---\n\n\n\n::: {.cell}\n\n:::\n\n::: {.cell}\n\n:::\n\n\nGeneralised linear models are fitted a little differently to standard linear models - namely, using maximum likelihood estimation instead of ordinary least squares for estimating the model coefficients.\n\nAs a result, we can no longer use F-tests for significance, or interpret $R^2$ values in quite the same way. This section will discuss new techniques for significance and goodness-of-fit testing, specifically for use with GLMs.\n\n## Libraries and functions\n\n::: {.callout-note collapse=\"true\"}\n## Click to expand\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n::: {.cell}\n\n```{.r .cell-code}\ninstall.packages(\"lmtest\")\nlibrary(lmtest)\n```\n:::\n\n## Python\n\n::: {.cell}\n\n```{.python .cell-code}\nfrom scipy.stats import *\n```\n:::\n\n:::\n\n:::\n\n## Deviance\n\nSeveral of the tests and metrics we'll discuss below are based heavily on deviance. So, what is deviance, and where does it come from?\n\nFitting a model using maximum likelihood estimation - the method that we use for GLMs, among other models - is all about finding the parameters that maximise the **likelihood**, or joint probability, of the sample. In other words, how likely is it that you would sample a set of data points like these, if they were being drawn from an underlying population where your model is true? Each model that you fit has its own likelihood.\n\nNow, for each dataset, there is a \"saturated\", or perfect, model. This model has the same number of parameters in it as there are data points, meaning the data are fitted exactly - as if connecting the dots between them. The **saturated model** has the largest possible likelihood of any model fitted to the dataset.\n\nOf course, we don't actually use the saturated model for drawing real conclusions, but we can use it as a baseline for comparison. We compare each model that we fit to this saturated model, to calculate the **deviance**. Deviance is defined as the difference between the log-likelihood of your fitted model and the log-likelihood of the saturated model (multiplied by 2). \n\nBecause deviance is all about capturing the discrepancy between fitted and actual values, it's performing a similar function to the residual sum of squares (RSS) in a standard linear model. In fact, the RSS is really just a specific type of deviance.\n\n![Different models and their deviances](images/deviance.png){width=70%}\n\n## Significance testing\n\nThere are a few different potential sources of p-values for a generalised linear model. \n\nHere, we'll briefly discuss the p-values that are reported \"as standard\" in a typical GLM model output.\n\nThen, we'll spend most of our time focusing on likelihood ratio tests, perhaps the most effective way to assess significance in a GLM.\n\n### Revisiting the diabetes dataset\n\nAs a worked example, we'll use a logistic regression fitted to the `diabetes` dataset that we saw in a previous section.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\ndiabetes <- read_csv(\"data/diabetes.csv\")\n```\n\n::: {.cell-output .cell-output-stderr}\n\n```\nRows: 728 Columns: 3\n── Column specification ────────────────────────────────────────────────────────\nDelimiter: \",\"\ndbl (3): glucose, diastolic, test_result\n\nℹ Use `spec()` to retrieve the full column specification for this data.\nℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\ndiabetes_py = pd.read_csv(\"data/diabetes.csv\")\n\ndiabetes_py.head()\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n glucose diastolic test_result\n0 148 72 1\n1 85 66 0\n2 183 64 1\n3 89 66 0\n4 137 40 1\n```\n\n\n:::\n:::\n\n:::\n\nAs a reminder, this dataset contains three variables:\n\n- `test_result`, binary results of a diabetes test result (1 for positive, 0 for negative)\n- `glucose`, the results of a glucose tolerance test\n- `diastolic` blood pressure\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_dia <- glm(test_result ~ glucose * diastolic,\n family = \"binomial\",\n data = diabetes)\n```\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = smf.glm(formula = \"test_result ~ glucose * diastolic\", \n family = sm.families.Binomial(), \n data = diabetes_py)\n \nglm_dia_py = model.fit()\n```\n:::\n\n:::\n\n### Wald tests\n\nLet's use the `summary` function to see the model we've just fitted.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_dia)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = test_result ~ glucose * diastolic, family = \"binomial\", \n data = diabetes)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) -8.5710565 2.7032318 -3.171 0.00152 **\nglucose 0.0547050 0.0209256 2.614 0.00894 **\ndiastolic 0.0423651 0.0363681 1.165 0.24406 \nglucose:diastolic -0.0002221 0.0002790 -0.796 0.42590 \n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for binomial family taken to be 1)\n\n Null deviance: 936.60 on 727 degrees of freedom\nResidual deviance: 748.01 on 724 degrees of freedom\nAIC: 756.01\n\nNumber of Fisher Scoring iterations: 4\n```\n\n\n:::\n:::\n\n\n## Python\n\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_dia_py.summary())\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n Generalized Linear Model Regression Results \n==============================================================================\nDep. Variable: test_result No. Observations: 728\nModel: GLM Df Residuals: 724\nModel Family: Binomial Df Model: 3\nLink Function: Logit Scale: 1.0000\nMethod: IRLS Log-Likelihood: -374.00\nDate: Fri, 17 May 2024 Deviance: 748.01\nTime: 12:43:46 Pearson chi2: 720.\nNo. Iterations: 5 Pseudo R-squ. (CS): 0.2282\nCovariance Type: nonrobust \n=====================================================================================\n coef std err z P>|z| [0.025 0.975]\n-------------------------------------------------------------------------------------\nIntercept -8.5711 2.703 -3.171 0.002 -13.869 -3.273\nglucose 0.0547 0.021 2.614 0.009 0.014 0.096\ndiastolic 0.0424 0.036 1.165 0.244 -0.029 0.114\nglucose:diastolic -0.0002 0.000 -0.796 0.426 -0.001 0.000\n=====================================================================================\n```\n\n\n:::\n:::\n\n:::\n\nWhichever language you're using, you may have spotted some p-values being reported directly here in the model summaries. Specifically, each individual parameter, or coefficient, has its own z-value and associated p-value.\n\nA hypothesis test has automatically been performed for each of the parameters in your model, including the intercept and interaction. In each case, something called a **Wald test** has been performed.\n\nThe null hypothesis for these Wald tests is that the value of the coefficient = 0. The idea is that if a coefficient isn't significantly different from 0, then that parameter isn't useful and could be dropped from the model. These tests are the equivalent of the t-tests that are calculated as part of the `summary` output for standard linear models.\n\nImportantly, these Wald tests *don't* tell you about the significance of the overall model. For that, we're going to need something else: a likelihood ratio test.\n\n### Likelihood ratio tests (LRTs)\n\nWhen we were assessing the significance of standard linear models, we were able to use the F-statistic to determine:\n\n- the significance of the model versus a null model, and\n- the significance of individual predictors.\n\nWe can't use these F-tests for GLMs, but we can use LRTs in a really similar way, to calculate p-values for both the model as a whole, and for individual variables.\n\nThese tests are all built on the idea of deviance, or the likelihood ratio, as discussed above on this page. We can compare any two models fitted to the same dataset by looking at the difference in their deviances, also known as the difference in their log-likelihoods, or more simply as a likelihood ratio.\n\nHelpfully, this likelihood ratio approximately follows a chi-square distribution, which we can capitalise on that to calculate a p-value. All we need is the number of degrees of freedom, which is equal to the difference in the number of parameters of the two models you're comparing.\n\n::: {.callout-warning}\nImportantly, we are only able to use this sort of test when one of the two models that we are comparing is a \"simpler\" version of the other, i.e., one model has a subset of the parameters of the other model. \n\nSo while we could perform an LRT just fine between these two models: `Y ~ A + B + C` and `Y ~ A + B + C + D`, or between any model and the null (`Y ~ 1`), we would not be able to use this test to compare `Y ~ A + B + C` and `Y ~ A + B + D`.\n:::\n\n#### Testing the model versus the null\n\nSince LRTs involve making a comparison between two models, we must first decide which models we're comparing, and check that one model is a \"subset\" of the other.\n\nLet's use an example from a previous section of the course, where we fitted a logistic regression to the `diabetes` dataset. \n\n::: {.panel-tabset group=\"language\"}\n## R\n\nThe first step is to create the two models that we want to compare: our original model, and the null model (with and without predictors, respectively).\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_dia <- glm(test_result ~ glucose * diastolic,\n family = \"binomial\",\n data = diabetes)\n\nglm_null <- glm(test_result ~ 1, \n family = binomial, \n data = diabetes)\n```\n:::\n\n\nThen, we use the `lrtest` function from the `lmtest` package to perform the test itself; we include both the models that we want to compare, listing them one after another.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlrtest(glm_dia, glm_null)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nLikelihood ratio test\n\nModel 1: test_result ~ glucose * diastolic\nModel 2: test_result ~ 1\n #Df LogLik Df Chisq Pr(>Chisq) \n1 4 -374.0 \n2 1 -468.3 -3 188.59 < 2.2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n```\n\n\n:::\n:::\n\n\nWe can see from the output that our chi-square statistic is significant, with a really small p-value. This tells us that, for the difference in degrees of freedom (here, that's 3), the change in deviance is actually quite big. (In this case, you can use `summary(glm_dia)` to see those deviances - 936 versus 748!)\n\nIn other words, our model is better than the null.\n\n## Python\n\nThe first step is to create the two models that we want to compare: our original model, and the null model (with and without our predictor, respectively).\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = smf.glm(formula = \"test_result ~ glucose * diastolic\", \n family = sm.families.Binomial(), \n data = diabetes_py)\n \nglm_dia_py = model.fit()\n\nmodel = smf.glm(formula = \"test_result ~ 1\",\n family = sm.families.Binomial(),\n data = diabetes_py)\n\nglm_null_py = model.fit()\n```\n:::\n\n\nUnlike in R, there isn't a nice neat function for extracting the $\\chi^2$ value, so we have to do a little bit of work by hand.\n\n\n::: {.cell}\n\n```{.python .cell-code}\n# calculate the likelihood ratio (i.e. the chi-square value)\nlrstat = -2*(glm_null_py.llf - glm_dia_py.llf)\n\n# calculate the associated p-value\npvalue = chi2.sf(lrstat, glm_dia_py.df_model - glm_null_py.df_model)\n\nprint(lrstat, pvalue)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n188.59314837444526 1.2288700360045209e-40\n```\n\n\n:::\n:::\n\n\nThis gives us the likelihood ratio, based on the log-likelihoods that we've extracted directly from the models, which approximates a chi-square distribution. \n\nWe've also calculated the associated p-value, by providing the difference in degrees of freedom between the two models (in this case, that's simply 1, but for more complicated models it's easier to extract the degrees of freedom directly from the model as we've done here).\n\nHere, we have a large chi-square statistic and a small p-value. This tells us that, for the difference in degrees of freedom (here, that's 1), the change in deviance is actually quite big. (In this case, you can use `summary(glm_dia)` to see those deviances - 936 versus 748!)\n\nIn other words, our model is better than the null.\n:::\n\n### Testing individual predictors\n\nAs well as testing the overall model versus the null, we might want to test particular predictors to determine whether they are individually significant.\n\nThe way to achieve this is essentially to perform a series of \"targeted\" likelihood ratio tests. In each LRT, we'll compare two models that are almost identical - one with, and one without, our variable of interest in each case.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\nThe first step is to construct a new model that doesn't contain our predictor of interest. Let's test the `glucose:diastolic` interaction.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nglm_dia_add <- glm(test_result ~ glucose + diastolic,\n family = \"binomial\",\n data = diabetes)\n```\n:::\n\n\nNow, we use the `lrtest` function to compare the models with and without the interaction:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nlrtest(glm_dia, glm_dia_add)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nLikelihood ratio test\n\nModel 1: test_result ~ glucose * diastolic\nModel 2: test_result ~ glucose + diastolic\n #Df LogLik Df Chisq Pr(>Chisq)\n1 4 -374.00 \n2 3 -374.32 -1 0.6288 0.4278\n```\n\n\n:::\n:::\n\n\nThis tells us that our interaction `glucose:diastolic` isn't significant - our more complex model doesn't have a meaningful reduction in deviance.\n\nThis might, however, seem like a slightly clunky way to test each individual predictor. Luckily, we can also use our trusty `anova` function with an extra argument to tell us about individual predictors. \n\nBy specifying that we want to use a chi-squared test, we are able to construct an analysis of deviance table (as opposed to an analysis of variance table) that will perform the likelihood ratio tests for us for each predictor:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nanova(glm_dia, test=\"Chisq\")\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nAnalysis of Deviance Table\n\nModel: binomial, link: logit\n\nResponse: test_result\n\nTerms added sequentially (first to last)\n\n Df Deviance Resid. Df Resid. Dev Pr(>Chi) \nNULL 727 936.60 \nglucose 1 184.401 726 752.20 < 2e-16 ***\ndiastolic 1 3.564 725 748.64 0.05905 . \nglucose:diastolic 1 0.629 724 748.01 0.42779 \n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n```\n\n\n:::\n:::\n\n\nYou'll spot that the p-values we get from the analysis of deviance table match the p-values you could calculate yourself using `lrtest`; this is just more efficient when you have a complex model!\n\n## Python\n\nThe first step is to construct a new model that doesn't contain our predictor of interest. Let's test the `glucose:diastolic` interaction.\n\n\n::: {.cell}\n\n```{.python .cell-code}\nmodel = smf.glm(formula = \"test_result ~ glucose + diastolic\", \n family = sm.families.Binomial(), \n data = diabetes_py)\n \nglm_dia_add_py = model.fit()\n```\n:::\n\n\nWe'll then use the same code we used above, to compare the models with and without the interaction:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nlrstat = -2*(glm_dia_add_py.llf - glm_dia_py.llf)\n\npvalue = chi2.sf(lrstat, glm_dia_py.df_model - glm_dia_add_py.df_model)\n\nprint(lrstat, pvalue)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n0.6288201373599804 0.42778842576800746\n```\n\n\n:::\n:::\n\n\nThis tells us that our interaction `glucose:diastolic` isn't significant - our more complex model doesn't have a meaningful reduction in deviance.\n:::\n\n## Goodness-of-fit\n\nGoodness-of-fit is all about how well a model fits the data, and typically involves summarising the discrepancy between the actual data points, and the fitted/predicted values that the model produces.\n\nThough closely linked, it's important to realise that goodness-of-fit and significance don't come hand-in-hand automatically: we might find a model that is significantly better than the null, but is still overall pretty rubbish at matching the data. So, to understand the quality of our model better, we should ideally perform both types of test. \n\n### Chi-square tests\n\nOnce again, we can make use of deviance and chi-square tests, this time to assess goodness-of-fit.\n\nAbove, we used likelihood ratio tests to assess the null hypothesis that our candidate fitted model and the null model had the same deviance.\n\nNow, however, we will test the null hypothesis that the fitted model and the saturated (perfect) model have the same deviance, i.e., that they both fit the data equally well. In most hypothesis tests, we want to reject the null hypothesis, but in this case, we'd like it to be true.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\nRunning a goodness-of-fit chi-square test in R can be done using the `pchisq` function. We need to include two arguments: 1) the residual deviance, and 2) the residual degrees of freedom. Both of these can be found in the `summary` output, but you can use the `$` syntax to call these properties directly like so:\n\n\n::: {.cell}\n\n```{.r .cell-code}\n1 - pchisq(glm_dia$deviance, glm_dia$df.residual)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 0.2605931\n```\n\n\n:::\n:::\n\n\n## Python\n\nThe syntax is very similar to the LRT we ran above, but now instead of including information about both our candidate model and the null, we instead just need 1) the residual deviance, and 2) the residual degrees of freedom:\n\n\n::: {.cell}\n\n```{.python .cell-code}\npvalue = chi2.sf(glm_dia_py.deviance, glm_dia_py.df_resid)\n\nprint(pvalue)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n0.26059314630406843\n```\n\n\n:::\n:::\n\n:::\n\nYou can think about this p-value, roughly, as \"the probability that this model is good\". We're not below our significance threshold, which means that we're not rejecting our null hypothesis (which is a good thing) - but it's also not a huge probability. This suggests that there's probably other variables we could measure and include in a future experiment, to give a better overall model.\n\n### AIC values\n\nYou might remember AIC values from standard linear modelling. AIC values are useful, because they tell us about overall model quality, factoring in both goodness-of-fit and model complexity.\n\nOne of the best things about the Akaike information criterion (AIC) is that it isn't specific to linear models - it works for models fitted with maximum likelihood estimation.\n\nIn fact, if you look at the formula for AIC, you'll see why:\n\n$$\nAIC = 2k - 2ln(\\hat{L})\n$$\n\nwhere $k$ represents the number of parameters in the model, and $\\hat{L}$ is the maximised likelihood function. In other words, the two parts of the equation represent the complexity of the model, versus the log-likelihood.\n\nThis means that AIC can be used for model comparison for GLMs in precisely the same way as it's used for linear models: lower AIC indicates a better-quality model.\n\n::: {.panel-tabset group=\"language\"}\n## R\n\nThe AIC value is given as standard, near the bottom of the `summary` output (just below the deviance values). You can also print it directly using the `$` syntax:\n\n\n::: {.cell}\n\n```{.r .cell-code}\nsummary(glm_dia)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall:\nglm(formula = test_result ~ glucose * diastolic, family = \"binomial\", \n data = diabetes)\n\nCoefficients:\n Estimate Std. Error z value Pr(>|z|) \n(Intercept) -8.5710565 2.7032318 -3.171 0.00152 **\nglucose 0.0547050 0.0209256 2.614 0.00894 **\ndiastolic 0.0423651 0.0363681 1.165 0.24406 \nglucose:diastolic -0.0002221 0.0002790 -0.796 0.42590 \n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n\n(Dispersion parameter for binomial family taken to be 1)\n\n Null deviance: 936.60 on 727 degrees of freedom\nResidual deviance: 748.01 on 724 degrees of freedom\nAIC: 756.01\n\nNumber of Fisher Scoring iterations: 4\n```\n\n\n:::\n\n```{.r .cell-code}\nglm_dia$aic\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n[1] 756.0069\n```\n\n\n:::\n:::\n\n\nIn even better news for R users, the `step` function works for GLMs just as it does for linear models, so long as you include the `test = LRT` argument.\n\n\n::: {.cell}\n\n```{.r .cell-code}\nstep(glm_dia, test = \"LRT\")\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\nStart: AIC=756.01\ntest_result ~ glucose * diastolic\n\n Df Deviance AIC LRT Pr(>Chi)\n- glucose:diastolic 1 748.64 754.64 0.62882 0.4278\n 748.01 756.01 \n\nStep: AIC=754.64\ntest_result ~ glucose + diastolic\n\n Df Deviance AIC LRT Pr(>Chi) \n 748.64 754.64 \n- diastolic 1 752.20 756.20 3.564 0.05905 . \n- glucose 1 915.52 919.52 166.884 < 2e-16 ***\n---\nSignif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n```\n\n\n:::\n\n::: {.cell-output .cell-output-stdout}\n\n```\n\nCall: glm(formula = test_result ~ glucose + diastolic, family = \"binomial\", \n data = diabetes)\n\nCoefficients:\n(Intercept) glucose diastolic \n -6.49941 0.03836 0.01407 \n\nDegrees of Freedom: 727 Total (i.e. Null); 725 Residual\nNull Deviance:\t 936.6 \nResidual Deviance: 748.6 \tAIC: 754.6\n```\n\n\n:::\n:::\n\n\n## Python\n\nThe AIC value isn't printed as standard with the model summary, but you can access it easily like so:\n\n\n::: {.cell}\n\n```{.python .cell-code}\nprint(glm_dia_py.aic)\n```\n\n::: {.cell-output .cell-output-stdout}\n\n```\n756.0068586069744\n```\n\n\n:::\n:::\n\n:::\n\n### Pseudo r-squared\n\nWe can't use $R^2$ values to represent the amount of variance explained in a GLM. This is primarily because, while linear models are fitted by minimising the squared residuals, GLMs are fitted by maximising the likelihood - an entirely different procedure.\n\nHowever, because $R^2$ values are so useful in linear modelling, statisticians have developed something called a \"pseudo $R^2$\" for GLMs.\n\n::: {.callout-note}\n#### Debate about pseudo $R^2$ values\n\nThere are two main areas of debate:\n\n1. Which version of pseudo $R^2$ to use? \n\nThere are many. Some of the most popular are McFadden's, Nagelkerke's, Cox & Snell's, and Tjur's. They all have slightly different formulae and in some cases can give quite different results. [This post](https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faq-what-are-pseudo-r-squareds/) does a nice job of discussing some of them and providing some comparisons.\n\n2. Should pseudo $R^2$ values be calculated at all? \n\nWell, it depends what you want them for. Most statisticians tend to advise that pseudo $R^2$ values are only really useful for model comparisons (i.e., comparing different GLMs fitted to the same dataset). This is in contrast to the way that we use $R^2$ values in linear models, as a measure of effect size that is generalisable across studies.\n\nSo, if you choose to use pseudo $R^2$ values, try to be thoughtful about it; and avoid the temptation to over-interpret! \n:::\n\n## Summary\n\nLikelihood and deviance are very important in generalised linear models - not just for fitting the model via maximum likelihood estimation, but for assessing significance and goodness-of-fit. To determine the quality of a model and draw conclusions from it, it's important to assess both of these things.\n\n::: {.callout-tip}\n#### Key points\n\n- Deviance is the difference between predicted and actual values, and is calculated by comparing a model's log-likelihood to that of the perfect \"saturated\" model \n- Using deviance, likelihood ratio tests can be used in lieu of F-tests for generalised linear models\n- Similarly, a chi-square goodness-of-fit test can also be performed using likelihood/deviance\n- The Akaike information criterion is also based on likelihood, and can be used to compare the quality of GLMs fitted to the same dataset\n- Other metrics that may be of use are Wald test p-values and pseudo $R^2$ values\n:::\n", + "supporting": [ + "significance-testing_files" + ], + "filters": [ + "rmarkdown/pagebreak.lua" + ], + "includes": {}, + "engineDependencies": {}, + "preserve": {}, + "postProcess": true + } +} \ No newline at end of file diff --git a/_quarto.yml b/_quarto.yml index 0c1ef0f..2cc13c2 100644 --- a/_quarto.yml +++ b/_quarto.yml @@ -33,6 +33,7 @@ metadata-files: # - "materials/_sidebar.yml" book: + bread-crumbs: false search: location: sidebar favicon: "_extensions/cambiotraining/courseformat/img/university-of-cambridge-favicon.ico" diff --git a/_site/index.html b/_site/index.html index ff094d9..c2ea896 100644 --- a/_site/index.html +++ b/_site/index.html @@ -6,6 +6,7 @@ + Course overview