Homework 1, Generalized linear models
STA442 Methods of Applied Statistics
Due 8 Oct 2021
1 Question 1: School
The dataset described at http://www.bristol.ac.uk/cmm/learning/mmsoftware/data-rev.html#chem97
contains chemistry test scores for 31,000 individuals from 2,200 schools in the UK. The outcome of interest is
the “Average GCSE sore” variable, a student’s final high school grade coded as ‘grade’ below. The maximum
grade is 8 (in this dataset, the original grades have been scaled somehow), and it appears that the number of
grade points below the maximum (variable ‘y’ below) looks like a Gamma distribution (don’t take my word
for it, check yourself). Each row is an individual, the first column is an id variable for region and the second
column indicates the school the subject belongs to.
The question of interest is to identify the important components of variation in GCSE scores. How important
are differences between schools, differences between regions, sex differences, or how old the student is? All
the students completed high school at age 18, it is hypothesized that those born near the end of a calendar
year (and therefore starting school a few months older than those born at the start of a calendar year) have
a slight advantage and should achieve higher grades.
Consider the code below, some of which is intended to be helpful and some of which is deliberately misleading.
xFile = Pmisc::downloadIfOld("http://www.bristol.ac.uk/cmm/media/migrated/datasets.zip")
x = read.table(grep("chem97", xFile, value = TRUE), col.names = c("region",
"school", "indiv", "chem", "sexNum", "ageMonthC", "grade"))
x$sex = factor(x$sexNum, levels = c(0, 1), labels = c("M",
"F"))
x$age = (222 + x$ageMonthC)/12
x$y = pmax(0.05, 8 - x$grade)
library("INLA")
xres = inla(y ~ 0 + age + f(region), data = x, family = "gamma",
control.family = list(hyper = list(prec = list(prior = "loggamma",
param = c(1e-04, 1e-04)))))
Pmisc::priorPostSd(xres, group = "random")$summary
mean sd 0.025quant 0.5quant 0.975quant mode
SD for region 0.1284088 0.009897604 0.1097051 0.1281781 0.1485965 0.1279285
Pmisc::priorPostSd(xres, group = "family")$summary
mean sd 0.025quant
SD parameter for the Gamma observations 0.5586705 0.002150167 0.5544532
0.5quant 0.975quant mode
SD parameter for the Gamma observations 0.5586649 0.5629152 0.5586649
1. Define a statistical model suitable for answering the question of interest, justifying your choice of
model and explaining all the various terms in your model.
1
2. Give prior distributions to all unknown model parameters and justify your choice. A plausible assump-
tion is differences between regions or schools or sexes are unlikely to cause doubling of GCSE grades,
a 50% increase or decrease of grades is rare but could well happen, changes of 20% are more likely.
3. Write a one-paragraph summary of the main results of your analysis, accompanied by a nicely formatted
table of results.
2 Question 2: Smoke
Data from the 2014 American National Youth Tobacco Survey is available on http://pbrown.ca/teaching/appliedstats/data,
where there is an R version of the 2014 dataset smoke2014.RData, a pdf documentation file
NYTS2014-Codebook-p.pdf, and the code used to create the R version of the data smokingData2014.R.
You can obtain the data with:
smokeFile = "smokeDownload2014.RData"
if (!file.exists(smokeFile)) {
download.file("http://pbrown.ca/teaching/appliedstats/data/smoke2014.RData",
smokeFile)
}
(load(smokeFile))
[1] "smoke" "smokeFormats"
The smoke object is a data.frame containing the data, the smokeFormats gives some explanation of the vari-
ables. The colName and label columns of smokeFormats contain variable names in smoke and descriptions
respectively.
smokeFormats[smokeFormats[, "colName"] == "chewing_tobacco_snuff_or",
c("colName", "label")]
colName
150 chewing_tobacco_snuff_or
label
150 RECODE: Used chewing tobacco, snuff, or dip on 1 or more days in the past 30 days
Consider the following model and set of results
# get rid of 9-11 and 19 year olds and missing age and
# race
smokeSub = smoke[which(smoke$Age >= 12 & smoke$Age <= 18 &
!is.na(smoke$Race) & !is.na(smoke$chewing_tobacco_snuff_or) &
(!is.na(smoke$Sex))), ]
smokeSub$ageFac = relevel(factor(smokeSub$Age), "15")
smokeSub$y = as.numeric(smokeSub$chewing_tobacco_snuff_or)
lincombDf = do.call(expand.grid, lapply(smokeSub[, c("ageFac",
"Sex", "Race", "RuralUrban")], levels))
lincombDf$y = -99
lincombList = inla.make.lincombs(as.data.frame(model.matrix(y ~
ageFac * Sex * RuralUrban * Race, lincombDf)))
library("INLA", quietly = TRUE)
smokeModel = inla(y ~ ageFac * Sex * RuralUrban * Race +
f(state) + f(school), lincomb = lincombList, data = smokeSub,
family = "binomial")
Warning in writeChar(lines[i], fp, nchars = nchar(lines[i]), eos = NULL):
problem writing to connection
2
knitr::kable(1/sqrt(smokeModel$summary.hyper[, c(4, 5, 3)]),
digits = 3)
0.5quant 0.975quant 0.025quant
Precision for state 0.008 0.003 0.02
Precision for school 0.803 0.637 1.00
An interesting plot is Figure 1.
smokePred = smokeModel$summary.lincomb.derived[,
paste0(c(0.5, 0.025, 0.975), 'quant')]
smokePred = exp(smokePred)/(1+exp(smokePred))
smokePred$diff = smokePred$'0.975quant' - smokePred$'0.025quant'
lincombDf$Age = as.numeric(as.character(lincombDf$ageFac))
lincombDf$ageShift = lincombDf$Age + 0.06*(as.numeric(lincombDf$Race)-2) +
0.3*(lincombDf$RuralUrban == 'Urban')
Spch = c('Rural' = 15, 'Urban' = 1)
Scol = c(black = 'black', white = 'red', hispanic='blue')
toPlot = (lincombDf$Race %in% names(Scol)) & (smokePred$diff < 0.9) &
lincombDf$Sex == 'M'
lincombDfHere = lincombDf[toPlot,]
smokePredHere = smokePred[toPlot,]
plot(
lincombDfHere$ageShift,
smokePredHere$'0.5quant',
pch = Spch[as.character(lincombDfHere$RuralUrban)],
col = Scol[as.character(lincombDfHere$Race)],
# log='y',
ylim = c(0,max(smokePredHere)),
xlab='age', ylab='prob', #yaxt='n',
yaxs='i', bty='l')
#forY = 1/c(4,10,25,100,500)
#axis(2, at=forY, mapply(format, forY), las=1)
segments(lincombDfHere$ageShift, smokePredHere$'0.025quant',
lincombDfHere$ageShift, smokePredHere$'0.975quant',
col = Scol[as.character(lincombDfHere$Race)])
legend('topleft', bty='n',
ncol = 2,
pch=c(rep(NA, length(Scol)), Spch),
lty = rep(c(1,NA), c(length(Scol), length(Spch))),
col = c(Scol, rep('black', length(Spch))),
legend=c(names(Scol), names(Spch)))
Consider the following two hypotheses:
1. State-level differences in chewing tobacco usage amongst high school students are much larger than
differences between schools within a state. Tobacco chewing is believed to be strongly regional, very
popular in some states and rare in others.
2. If American TV is an accurate reflection of reality, tobacco chewing is mostly done by Cowboys.
Cowboys are Male, white and live in rural areas. Tobacco chewing is fairly common amongst White
rural American high school students of every age, and virtually unknown for other ethnic groups and
3
12 13 14 15 16 17 180
.0
0
0.
10
0.
20
0.
30
age
pr
ob
black
white
hispanic
Rural
Urban
Figure 1: Plot of predicted probabilities
in urban areas.
Write a short consulting report addressing these hypotheses. This should include the following:
• a one-paragraph summary stating your conclusions, which could be understood by a child health and
welfare professional or an executive in the marketing department of a large tobacco firm;
• a writeup of roughly one page of text (not including figures and tables) containing
– an introduction restating the problem as you’ve interpreted it in relation to this dataset,
– a methods section giving the statistical models used (in mathematical notation, not R syntax)
and justifying their use, and
– a results section where the results are described and interpreted; and
• an appendix containing your code.
The report will be assessed in terms of:
• clarity of presentation,
• the use of an appropriate model and providing jusification for it,
• demonstration of an understanding of the statistical models used, and
• drawing conclusions which are consistent with the analysis.
Some words of advice
• Write in sentences and paragraphs.
• Provide captions for ALL figures and tables
• Don’t use default axis labels on plots and ensure text on plots is large enough to read comfortably
• Round numbers to 2 or 3 decimal places so tables look tidy.
4
• Don’t show raw R output. Put things in Latex or Markdown tables (using knitr::kable or
Hmisc::latex)
• Give parameter estimates and confidence intervals on the ‘natural’ scale where possible (probabilities
or odds rather than log-odds ratios)
5

欢迎咨询51作业君