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

日期：2021-10-07 10:55

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