Skip to content

Commit dc3cf44

Browse files
committed
new notebook demonstrating overdispersion
1 parent aafe914 commit dc3cf44

File tree

1 file changed

+379
-0
lines changed

1 file changed

+379
-0
lines changed

Lecture-7-DataPreprocessing3.ipynb

Lines changed: 379 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,379 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"## Introduction\n",
8+
"\n",
9+
"In the last notebook, we introduced generalize linear models (GLMs) that can be used to analyze non-normal data, such as count data or percentages. In general, GLMs tend to be much more powerful than linear models. However, there are a few statistical assumptions that you should be aware of...\n",
10+
"\n",
11+
"For this notebook, we will continue to use the glucosinolate data that we downloaded earlier.\n",
12+
"\n",
13+
"### Load the glucosinolate data"
14+
]
15+
},
16+
{
17+
"cell_type": "code",
18+
"execution_count": null,
19+
"metadata": {},
20+
"outputs": [],
21+
"source": [
22+
"# rm(list=ls());\n",
23+
"# open the glucosinolate file in R\n",
24+
"# same file as before...\n",
25+
"glucosinolateFileName <- \"data/cmeyer_glucs2015/bmeyer_etal.txt\"; \n",
26+
"glucs <- read.table(glucosinolateFileName, header=T, sep=\"\\t\", as.is=T, stringsAsFactors=FALSE); \n",
27+
"glucs <- glucs[order(glucs[,\"accession_id\"]),];\n",
28+
"glucs[,3:ncol(glucs)] <- round(glucs[,3:ncol(glucs)]);"
29+
]
30+
},
31+
{
32+
"cell_type": "code",
33+
"execution_count": null,
34+
"metadata": {},
35+
"outputs": [],
36+
"source": [
37+
"# what's in the working environment?\n",
38+
"ls();"
39+
]
40+
},
41+
{
42+
"cell_type": "markdown",
43+
"metadata": {},
44+
"source": [
45+
"### What is the statistical model?\n",
46+
"\n",
47+
"With this notebook, we will continue to discuss *Poisson* family analyses. Poisson models (and their derivatives) are used to analyze count data, which - as you will remember - are very common.\n",
48+
"\n",
49+
"Examples:\n",
50+
" - The number of insects within several $1m^2$ quadrats\n",
51+
" - The number of ions that hit a detector in a time interval\n",
52+
" - The number of meteorites that hit the Earth each year\n",
53+
"\n",
54+
" \n",
55+
"The glucosinolate data are also count data (i.e. the number of glucosinolate ions that hit the QQQ detector in a given time frame)."
56+
]
57+
},
58+
{
59+
"cell_type": "code",
60+
"execution_count": null,
61+
"metadata": {},
62+
"outputs": [],
63+
"source": [
64+
"str(glucs)"
65+
]
66+
},
67+
{
68+
"cell_type": "code",
69+
"execution_count": null,
70+
"metadata": {},
71+
"outputs": [],
72+
"source": [
73+
"head(glucs);"
74+
]
75+
},
76+
{
77+
"cell_type": "code",
78+
"execution_count": null,
79+
"metadata": {},
80+
"outputs": [],
81+
"source": [
82+
"# if lme4, ggplot2, and gridExtra aren't installed, install them...\n",
83+
"if( !require(\"lme4\" )){ \n",
84+
" install.packages(\"lme4\"); \n",
85+
"}\n",
86+
"\n",
87+
"if( !require(\"ggplot2\" )){ \n",
88+
" install.packages(\"ggplot2\"); \n",
89+
"}"
90+
]
91+
},
92+
{
93+
"cell_type": "markdown",
94+
"metadata": {},
95+
"source": [
96+
"\n",
97+
"Recall that we used linear and generalized linear mixed models (i.e. fitted and random effects) to model BLUPs. As an example with a Poisson GLM:"
98+
]
99+
},
100+
{
101+
"cell_type": "code",
102+
"execution_count": null,
103+
"metadata": {},
104+
"outputs": [],
105+
"source": [
106+
"# previously, we used the following code to determine the total 'effort' or offset\n",
107+
"totalIonCounts <- rowSums( glucs[,-c(1,2)] );\n",
108+
"cat(\"The range of ion counts is now:\", range(totalIonCounts), \"\\n\");"
109+
]
110+
},
111+
{
112+
"cell_type": "code",
113+
"execution_count": null,
114+
"metadata": {},
115+
"outputs": [],
116+
"source": [
117+
"# the 0 indicates that there are missing data; get rid of those samples:\n",
118+
"dropouts <- which( totalIonCounts == 0 );\n",
119+
"glucs2 <- glucs[-dropouts,];\n",
120+
"dim(glucs2);\n",
121+
"\n",
122+
"totalIonCounts <- rowSums( glucs2[,-c(1,2)] );\n",
123+
"cat(\"The range of ion counts is now:\", range(totalIonCounts), \"\\n\");"
124+
]
125+
},
126+
{
127+
"cell_type": "code",
128+
"execution_count": null,
129+
"metadata": {},
130+
"outputs": [],
131+
"source": [
132+
"# we can now fit the mixed model:\n",
133+
"glmer0 <- glmer( G2P ~ 1 + offset(log(totalIonCounts)) + (1|accession_id), data=glucs2, control=glmerControl(optimizer=\"bobyqa\", optCtrl=list(maxfun=2e5)), family=\"poisson\" );\n",
134+
"glmBlups <- ranef(glmer0)$accession_id; # these are the best-linear unbiased predictors:\n",
135+
"glmBlups <- data.frame( accession_id=rownames(glmBlups), pois_blup=glmBlups[,1], stringsAsFactors=FALSE );\n",
136+
"tail(glmBlups);"
137+
]
138+
},
139+
{
140+
"cell_type": "markdown",
141+
"metadata": {},
142+
"source": [
143+
"### Model assumptions\n",
144+
"\n",
145+
"The Poisson distribution is used with non-negative integers that have no natural upper bound. However, there are other assumptions, including independence among events. For example, the number of people buying lunch each minute is not really Poisson distributed, because (a) people go to lunch together and (b) people eat at certain times (among other factors).\n",
146+
"\n",
147+
"Furthermore, the mean and variance of a Poisson random variable are expected to be equal and are described by the single parameter $\\lambda$. \n",
148+
"\n",
149+
"When the mean and variance are *not* equal, then the data are either underdispersed or overdispersed. Underdispersion (mean > variance) is rare compared to overdispersion (mean < variance) and is often ignored.\n",
150+
"\n",
151+
"<!--# however, these are count data (see the output from head above), which suggests \n",
152+
"# that a Poisson-family model should be used; let's try a simple quasi-Poisson GLM:\n",
153+
"glm0 <- glm( G2P ~ 1, data=glucs, family=\"quasipoisson\" );\n",
154+
"glm0.res <- residuals( glm0 );\n",
155+
"glm0.res[1:10];\n",
156+
"\n",
157+
"# note: the names are no longer the accession Ids, but the row names from the glucs data frame. Do you know why?\n",
158+
"# let's use a workaround:\n",
159+
"################################################################################\n",
160+
"## my version of stack, which avoids factor generation\n",
161+
"################################################################################\n",
162+
"mstack <- function(arg, newHeaders, setRowNames=T, sorted=TRUE, decreasing=F){\n",
163+
" values <- data.frame(names=I(names(arg)), values=as.numeric(arg));\n",
164+
"\n",
165+
" if( setRowNames ){\n",
166+
" rownames(values) <- values[,\"names\"];\n",
167+
" }\n",
168+
" \n",
169+
" if( sorted ) {\n",
170+
" values <- values[order(values[,\"values\"], decreasing=decreasing),];\t\n",
171+
" }\n",
172+
" \n",
173+
" colnames(values) <- newHeaders;\n",
174+
" return(values);\n",
175+
"}\n",
176+
"\n",
177+
"glm0.res <- mstack( glm0.res, newHeaders=c(\"row_id\", \"residual\"), sorted=FALSE );\n",
178+
"head(glm0.res);-->\n",
179+
"\n",
180+
"We can test for overdispersion using the [following code from Ben Bolker](http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#testing-for-overdispersioncomputing-overdispersion-factor):"
181+
]
182+
},
183+
{
184+
"cell_type": "code",
185+
"execution_count": null,
186+
"metadata": {},
187+
"outputs": [],
188+
"source": [
189+
"# a test for overdispersion, from Ben Bolker's GLMM website:\n",
190+
"overdisp_fun <- function(model) {\n",
191+
" rdf <- df.residual(model)\n",
192+
" rp <- residuals(model,type=\"pearson\")\n",
193+
" Pearson.chisq <- sum(rp^2)\n",
194+
" prat <- Pearson.chisq/rdf\n",
195+
" pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)\n",
196+
" c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)\n",
197+
"}"
198+
]
199+
},
200+
{
201+
"cell_type": "code",
202+
"execution_count": null,
203+
"metadata": {},
204+
"outputs": [],
205+
"source": [
206+
"# determine whether the data are overdispersed:\n",
207+
"glmer0 <- glmer( G2P ~ 1 + offset(log(totalIonCounts)) + (1|accession_id), data=glucs2, family=\"poisson\" );\n",
208+
"overdisp_fun( glmer0 );"
209+
]
210+
},
211+
{
212+
"cell_type": "markdown",
213+
"metadata": {},
214+
"source": [
215+
"The data are highly overdispersed.\n"
216+
]
217+
},
218+
{
219+
"cell_type": "code",
220+
"execution_count": null,
221+
"metadata": {},
222+
"outputs": [],
223+
"source": [
224+
"# one workaround is to simply update the coefficient table by multiplying the\n",
225+
"# standard error by the square root of the dispersion factor (phi) and then \n",
226+
"# updating the Z-values and the corresponding P-values such as:\n",
227+
"\n",
228+
"modelToTest <- coef(summary(glmer0));\n",
229+
"phi <- overdisp_fun(glmer0)[\"ratio\"];\n",
230+
"\n",
231+
"modelToTest <- within(as.data.frame(modelToTest),\n",
232+
"{ `Std. Error` <- `Std. Error`*sqrt(phi)\n",
233+
" `z value` <- Estimate/`Std. Error`\n",
234+
" `Pr(>|z|)` <- 2*pnorm(abs(`z value`), lower.tail=FALSE)\n",
235+
"})\n",
236+
"\n",
237+
"cat(\"Before:\\n\");\n",
238+
"coef(summary(glmer0));\n",
239+
"\n",
240+
"cat(\"\\nAnd after:\\n\");\n",
241+
"printCoefmat(modelToTest, digits=3);"
242+
]
243+
},
244+
{
245+
"cell_type": "code",
246+
"execution_count": null,
247+
"metadata": {},
248+
"outputs": [],
249+
"source": [
250+
"# another workaround: simply use glmer.nb as an alternative\n",
251+
"# however glmer.nb is \"slow and fragile\" (see the Bolker wiki page noted above):\n",
252+
"glmer.nb0 <- glmer.nb( G2P ~ 1 + offset(log(totalIonCounts)) + (1|accession_id), data=glucs2);"
253+
]
254+
},
255+
{
256+
"cell_type": "code",
257+
"execution_count": null,
258+
"metadata": {},
259+
"outputs": [],
260+
"source": [
261+
"# can we use quasi models with mixed models? Not with lme4. One alternative is to fit an olre\n",
262+
"# or an 'observation level random effect'. This is somewhat conservative, but will aid you in accounting for overdispersion\n",
263+
"glucs2 <- data.frame(olre=1:nrow(glucs2), glucs2, stringsAsFactors=FALSE);\n",
264+
"glucs2[1:3,];"
265+
]
266+
},
267+
{
268+
"cell_type": "code",
269+
"execution_count": null,
270+
"metadata": {},
271+
"outputs": [],
272+
"source": [
273+
"# with the olre as a random effect:\n",
274+
"glmer.olre <- glmer( G2P ~ 1 + offset(log(totalIonCounts)) + (1|olre) + (1|accession_id), data=glucs2, family=\"poisson\");"
275+
]
276+
},
277+
{
278+
"cell_type": "markdown",
279+
"metadata": {},
280+
"source": [
281+
"\n",
282+
"### Last but not least, in a pinch, you could use residuals instead of BLUPs"
283+
]
284+
},
285+
{
286+
"cell_type": "code",
287+
"execution_count": null,
288+
"metadata": {},
289+
"outputs": [],
290+
"source": [
291+
"# another possible workaround\n",
292+
"glm0 <- glm( G2P ~ 1, offset=log(totalIonCounts), data=glucs2, family=\"quasipoisson\");\n",
293+
"summary(glm0);"
294+
]
295+
},
296+
{
297+
"cell_type": "code",
298+
"execution_count": null,
299+
"metadata": {},
300+
"outputs": [],
301+
"source": [
302+
"# like BLUPs, residuals are regularly used in GWAS\n",
303+
"glm0.res <- residuals( glm0 );\n",
304+
"glm0.res[1:10];"
305+
]
306+
},
307+
{
308+
"cell_type": "code",
309+
"execution_count": null,
310+
"metadata": {},
311+
"outputs": [],
312+
"source": [
313+
"# note: the names are no longer the accession Ids, but the row names from the glucs data frame. Do you know why?\n",
314+
"# let's use a workaround:\n",
315+
"################################################################################\n",
316+
"## my version of stack, which avoids factor generation\n",
317+
"################################################################################\n",
318+
"mstack <- function(arg, newHeaders, setRowNames=T, sorted=TRUE, decreasing=F){\n",
319+
" values <- data.frame(names=I(names(arg)), values=as.numeric(arg));\n",
320+
"\n",
321+
" if( setRowNames ){\n",
322+
" rownames(values) <- values[,\"names\"];\n",
323+
" }\n",
324+
" \n",
325+
" if( sorted ) {\n",
326+
" values <- values[order(values[,\"values\"], decreasing=decreasing),];\t\n",
327+
" }\n",
328+
" \n",
329+
" colnames(values) <- newHeaders;\n",
330+
" return(values);\n",
331+
"}\n",
332+
"\n",
333+
"glm0.res <- mstack( glm0.res, newHeaders=c(\"row_id\", \"residual\"), sorted=FALSE );\n",
334+
"glucs2 <- data.frame( row_id=rownames(glucs2), glucs2, stringsAsFactors=FALSE);\n",
335+
"\n",
336+
"# the residuals here only correspond to the metabolite we've been working with, so \n",
337+
"# to protect yourself from making downstream mistakes, subset to the metadata and that column\n",
338+
"output <- merge(glm0.res, glucs2[,c(\"row_id\", \"accession_id\", \"sample_weight\", \"G2P\")], by=\"row_id\" );"
339+
]
340+
},
341+
{
342+
"cell_type": "code",
343+
"execution_count": null,
344+
"metadata": {
345+
"scrolled": true
346+
},
347+
"outputs": [],
348+
"source": [
349+
"head(output);"
350+
]
351+
},
352+
{
353+
"cell_type": "markdown",
354+
"metadata": {},
355+
"source": [
356+
"<br><br>\n",
357+
"\n",
358+
"One could then take the averages of each residual for each accession id, using an approach like the tapply approach described earlier."
359+
]
360+
}
361+
],
362+
"metadata": {
363+
"kernelspec": {
364+
"display_name": "R",
365+
"language": "R",
366+
"name": "ir"
367+
},
368+
"language_info": {
369+
"codemirror_mode": "r",
370+
"file_extension": ".r",
371+
"mimetype": "text/x-r-source",
372+
"name": "R",
373+
"pygments_lexer": "r",
374+
"version": "3.4.3"
375+
}
376+
},
377+
"nbformat": 4,
378+
"nbformat_minor": 4
379+
}

0 commit comments

Comments
 (0)