STA442程序语言代写、代写R编程

日期：2021-10-08 10:14

STA442 Methods of Applied Statistics

Due 8 Oct 2021

1 Question 1: School

The dataset described at

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,

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)