readNHANES.R: Data Preprocessing

1 library(XML)

2 attach(file.path("ExpoCast-040813.RData"))

3 ## Code to update NHANES, ExpoCast.NHANES.compare, and 

4 NHANES.parent.metabolite.mapping

5 ## In NHANES:

6 ## Misspellings, etc.

7 ## 8:2-Hydroxynapthalene (should be phth)

8 ## 54: should be 1-Hydroxynaphthalene (original has pthalene)

9 ## 

10 ## 3- and 9- Hydroxybenz[a]anthracene are combined in the data tables.

11 ## 29: should be 3- and 9-Hydroxybenz[a]anthracene

12 ## 58: 9-Hydroxybenz[a]anthracene - drop

13

14 ## In NHANES.parent.metabolite.mapping$Name.1

15 ## fix the same two ...-napthalene compounds

16 ##

17 options(stringsAsFactors=FALSE)

18

19 NHANES$Name <- sub("napthalene","naphthalene",NHANES$Name)

20 NHANES.parent.metabolite.mapping$Name.1 <-

21 sub("napthalene","naphthalene", NHANES.parent.metabolite.mapping$Name.1)

22

23 indx <- grep("9-Hydroxybenz\\[a\\]anthracene", NHANES$Name)

24 NHANES <- NHANES[-indx,]

25

26 indx <- grep("9-Hydroxybenz\\[a\\]anthracene", 

27 NHANES.parent.metabolite.mapping$Name.1)

28 NHANES.parent.metabolite.mapping <- NHANES.parent.metabolite.mapping[-indx,]

29

30 NHANES$Name <- sub("3-Hydroxybenz\\[a\\]anthracene","3- and 9-

31 Hydroxybenz[a]anthracene",NHANES$Name)

32 NHANES.parent.metabolite.mapping$Name.1 <-

33 sub("3-Hydroxybenz\\[a\\]anthracene",

34 "3- and 9-Hydroxybenz[a]anthracene", NHANES.parent.metabolite.mapping$Name.1)

35

36 ## Triclosan. In NHANES, this includes the more technical chemical name, encoded

37 ## in a way I cannot decipher. This is just called Triclosan in 

38 ExpoCast.NHANES.compare,

39 ## so, make it that way in NHANES, too.

40

41 indx <- grep("Triclosan",NHANES$Name)

42 NHANES$Name[indx] <- "Triclosan"

43 ## ====================================================================

44 ## Utility function to return IUR production volume code from geometric

45 ## mean production volume

46 geometric.mean.to.IUR <-

47 function (x) 

48 {

49 if (abs(x - (exp((log(1 * 1e+06/2.204/365.25) + log(10 * 1e+06/2.204/365.25))/2))) < 

50 1e-14) {

51 string <- "1 to < 10 million lbs"

52 lower <- 1 * 1e+06/2.204/365.25

53 upper <- 10 * 1e+06/2.204/365.25

54 }

55 else if (abs(x - (exp((log(0.5 * 1e+06/2.204/365.25) + log(1 * 

56 1e+06/2.204/365.25))/2))) < 

57 1e-14) {

58 string <- "500,000 to < 1 million lbs"

59 lower <- 0.5 * 1e+06/2.204/365.25

60 upper <- 1 * 1e+06/2.204/365.25
61 }

62 else if (abs(x - exp((log(25000/2.204/365.25) + log(0.5 * 1e+06/2.204/365.25))/2)) < 

63 1e-14) {

64 string <- "25,000 to < 500,000 lbs"

65 lower <- 25000/2.204/365.25

66 upper <- 0.5 * 1e+06/2.204/365.25

67 }

68 else if (abs(x - exp((log(6250/2.204/365.25) + log(25000/2.204/365.25))/2)) < 1e-14) 

69 {

70 string <- "< 25,000 lbs"

71 lower <- 6250/2.204/365.25

72 upper <- 25000/2.204/365.25

73 } 

74 else if (abs(x - exp((log(10 * 1e+06/2.204/365.25) + log(50 *

75 1e+06/2.204/365.25))/2)) < 1e-14) {

76 string <- "10 to < 50 million lbs"

77 lower <- 10 * 1e+06/2.204/365.25

78 upper <- 50 * 1e+06/2.204/365.25

79 }

80 else if (abs(x - exp((log(50 * 1e+06/2.204/365.25) + log(100 * 

81 1e+06/2.204/365.25))/2)) < 1e-14) {

82 string <- "50 to < 100 million lbs"

83 lower <- 50 * 1e+06/2.204/365.25

84 upper <- 100 * 1e+06/2.204/365.25

85 }

86 else if (abs(x - exp((log(100 * 1e+06/2.204/365.25) + log(500 * 

87 1e+06/2.204/365.25))/2)) < 1e-14) {

88 string <- "100 to < 500 million lbs"

89 lower <- 100 * 1e+06/2.204/365.25

90 upper <- 500 * 1e+06/2.204/365.25

91 }

92 else if (abs(x - exp((log(500 * 1e+06/2.204/365.25) + log(1000 * 

93 1e+06/2.204/365.25))/2)) < 1e-14) {

94 string <- "500 million to < 1 billion lbs"

95 lower <- 500 * 1e+06/2.204/365.25

96 upper <- 1000 * 1e+06/2.204/365.25

97 }

98 else if (abs(x - 5000 * 1e+06/2.204/365.25) < 1e-14) {

99 string <- "1 billion lbs and greater"

100 lower <- 1000 * 1e+06/2.204/365.25

101 upper <- Inf

102 }

103 else {

104 string <- x

105 lower <- NA

106 upper <- NA

107 }

108 out <- data.frame(string)

109 out <- cbind(out, data.frame(lower))

110 out <- cbind(out, data.frame(upper))

111 colnames(out) <- c("IUR Code", "Lower", "Upper")

112 return(out)

113 }

114

115 ## =====================================================================

116

117 tmp <- scan("CDC-NHANES-2005.xml",character(0), sep="\n")

118

119 NHANES <- subset(NHANES, !Inferred)

120

121 get2005Table <- function(name){

122 ## escape special characters in name

123 name <- gsub("\\[","\\\\[",gsub("]","\\\\]",name))
124 name <- gsub("\\(","\\\\(",gsub(")","\\\\)",name))

125 indx <- grep(paste("Table +[[:digit:]]{1,4}\\. +",name," 

126 *.*\\(creatinine",sep=""),tmp)

127 if (length(indx) > 0) {

128 indx2 <- grep("SurveySampleyearssizeTotal, age 6 and 

129 older",tmp[indx:length(tmp)])[1]

130 dt <- sub("^ *SurveySampleyearssizeTotal, age 6 and older(.+?)Age 

131 group.*$","\\1",tmp[indx - 1 + indx2[1]])

132 ## need to pull out a unit indicator

133 dunitmult <- sub(".+\\(in (.+)g/g of creatinine\\).+$","\\1",tmp[indx - 1 + 

134 indx2[1]])

135 ## replace all < LOD with NA (NA-NA)

136 dt <- gsub("< LOD","NA(NA-NA)",dt)

137 dt <- gsub("\\*","NA(NA-NA)",dt)

138 ## split the string at ##-##

139 dt <- gsub("([[:digit:]]+) +([[:digit:]][[:digit:]]\\-

140 [[:digit:]][[:digit:]])","\\1\n\\2",dt)

141 ## dt <-gsub("([[:digit:]][[:digit:]]\\-[[:digit:]][[:digit:]])","\\1 ",dt) 

142 dt <- strsplit(dt,"\n")[[1]]

143 ## Split the rows of dt into fields.

144 year <- character(length(dt))

145 for (i in seq(along=year)) {

146 indx <- regexpr("[[:digit:]][[:digit:]]-[[:digit:]][[:digit:]] *",dt[i])

147 year[i] <- sub(" *$","",substr(dt[i],1,attr(indx,"match.length")))

148 dt[i] <- substr(dt[i],attr(indx,"match.length")+1,nchar(dt[i]))

149 }

150 mytb <- do.call("rbind",strsplit(dt,")|\\(|-"))

151 ## Values as ug/g creatinine. When unitmult == "n", mult by 0.001

152 mult <- if (dunitmult == "n") 0.001 else 1

153 data.frame(sample=year,

154 gm=as.numeric(mytb[,1]) * mult,

155 gmlcl=as.numeric(mytb[,2]) * mult,

156 gmucl=as.numeric(mytb[,3]) * mult,

157 Value50=as.numeric(mytb[,4]) * mult,

158 Value75=as.numeric(mytb[,7]) * mult,

159 Value90=as.numeric(mytb[,10]) * mult,

160 Value95=as.numeric(mytb[,13]) * mult,

161 LCL50=as.numeric(mytb[,5]) * mult,

162 UCL50=as.numeric(mytb[,6]) * mult,

163 LCL75=as.numeric(mytb[,8]) * mult,

164 UCL75=as.numeric(mytb[,9]) * mult,

165 LCL90=as.numeric(mytb[,11]) * mult,

166 UCL90=as.numeric(mytb[,12]) * mult,

167 LCL95=as.numeric(mytb[,14]) * mult,

168 UCL95=as.numeric(mytb[,15]) * mult)

169 }

170 

171 }

172

173

174

175 ## ----------------------------------------------------------------------------------

176 ## Data from 2011 report

177 tmp2 <- scan("CDC-NHANES-2011.xml", character(0), sep="\n")

178 doc <- xmlParse(tmp2)

179 ## Find the beginning and ending of each table.

180 tables <- getNodeSet(doc,"//Table")

181 ## The first row contains the chemical labels

182 labels <- sapply(getNodeSet(doc,"//Table/TR[1]/TH/P"),xmlValue)

183 creaIndex <- grep("\\(creatinine",labels)

184

185 ## extract confidence limits from a vector of strings like (129-157)

186
187 getcl <- function(v,lim=c(1,2)) {

188 v <- gsub(" ","",v)

189 v <- gsub("\\(","",gsub("\\)","",v))

190 n <- strsplit(v,"-")

191 sapply(n,function(z) as.numeric(z[lim]))

192 }

193

194

195

196 get2011Table <- function(name){

197 name <- gsub("\\[","\\\\[",gsub("]","\\\\]",name))

198 name <- gsub("\\(","\\\\(",gsub(")","\\\\)",name))

199 name <- gsub(" "," *",name)

200 indx <- grep(name, labels[creaIndex])

201 if (length(indx) > 0) {

202 ## Get the units. All ug/g creatinine. 

203 tb <- tables[creaIndex][[indx]]

204 Rows <- getNodeSet(tb,"./TR[TD]")

205 mytb <- do.call("rbind",

206 lapply(Rows,

207 function(z) {

208 sapply(getNodeSet(z,"./TD/P"),xmlValue)

209 }

210 )

211 )

212 ## There will be a warning for one of the chemicals, because the last

213 ## row in the table has only two fields, for footnotes.

214 ## We want the rows associated with Total: row 1 up to the row where the

215 ## first field contains "Age"

216 indx <- grep("Age",mytb[,1])

217 mytb <- mytb[1:indx,c(2,3:12)]

218 ## Now, drop the spacer rows

219 indx <- grep("^ +",mytb[,1])

220 mytb <- mytb[-indx,]

221 ## Return a dataframe with nrows(mytb) rows and columns for

222 ## Valuex, LCLx, UCLx, where x is 50, 75, 90, and 95.

223 mydf <-

224 data.frame(

225 sample = gsub(" +","",mytb[,1]),

226 gm = as.numeric(mytb[,2]),

227 gmlcl = getcl(mytb[,3],lim=1),

228 gmucl = getcl(mytb[,3],lim=2),

229 Value50 = as.numeric(mytb[,4]),

230 Value75 = as.numeric(mytb[,6]),

231 Value90 = as.numeric(mytb[,8]),

232 Value95 = as.numeric(mytb[,10]),

233 LCL50 = getcl(mytb[,5],lim=1),

234 UCL50 = getcl(mytb[,5],lim=2),

235 LCL75 = getcl(mytb[,7],lim=1),

236 UCL75 = getcl(mytb[,7],lim=2),

237 LCL90 = getcl(mytb[,9],lim=1),

238 UCL90 = getcl(mytb[,9],lim=2),

239 LCL95 = getcl(mytb[,11],lim=1),

240 UCL95 = getcl(mytb[,11],lim=2)

241 )

242 mydf

243 }

244 }

245

246 ## Build the new NHANES object, NHANESnew, which has the chemical info from NHANES

247 ## with the quantile data we've just collected

248 dta <- list()

249
250 for (i in 1:nrow(NHANES)) {

251 dti <- if (NHANES$Year[i] == 2005){

252 get2005Table(NHANES$Name[i])

253 } else {

254 get2011Table(NHANES$Name[i])

255 }

256 dta[[i]] <- cbind(data.frame(Name=rep(NHANES$Name[i],nrow(dti)),

257 CAS=rep(NHANES$CAS[i],nrow(dti)),

258 MW=rep(NHANES$MW[i], nrow(dti)),

259 LOD=rep(NHANES$LOD[i],nrow(dti)),

260 Year=rep(NHANES$Year[i],nrow(dti)),

261 Value=rep(NHANES$Value[i],nrow(dti))

262 ),

263 dti, stringsAsFactors=FALSE)

264 }

265 NHANES.old <- NHANES

266 NHANES <- do.call("rbind", dta)

267 ## Before conversion, print out some representative values to compare back to the pdfs

268 ## Particularly in the 2005 dataset, the pdf -> xml conversion reorganized the

269 ## table labels and data.

270 write.csv(rbind(NHANES[NHANES$Year == 

271 2005,c("CAS","Name","gm","Value95","LCL95","UCL95")],

272 NHANES[NHANES$Year == 

273 2011,c("CAS","Name","gm","Value95","LCL95","UCL95")]),

274 file="NHANES_Values-2.csv",row.names=FALSE)

275

276 ## 

277 NHANES$sample <- sub("\\*$","",NHANES$sample)

278 NHANES$sample <- factor(NHANES$sample,levels=c("99-00","01-02","03-04","05-06","07-

279 08"),ordered=TRUE)

280

281

282 parent <- ExpoCast.NHANES.compare

283 ## parent inadvertantly includes arsenic. Drop that

284 indx <- grep("rsenic", parent$Compound)

285 parent <- parent[-indx,]

286

287 rownames(parent) <- as.character(parent$CASRN)

288

289 mapping <- NHANES.parent.metabolite.mapping

290 ## Drop columns CAS.Checksum, CAS.Checksum.1, X, X.1, and X.2

291 indx <- match(c("CAS.Checksum", "CAS.Checksum.1", "X", "X.1", "X.2"), 

292 colnames(mapping))

293 mapping <- mapping[,-indx]

294

295 ## This gives three data sets:

296 ## NHANES: metabolite population quantiles for possibly multiple NHANES cycles

297 ## mapping: the mapping between parent chemicals and the metabolites measured in 

298 NHANES

299 ## parent: data on release, use, and farfield exposure for parents linked to the

300 ## chemicals measured in NHANES.

301 ## We want to eliminate every parent that either does not have a product in NHANES or 

302 is not in

303 ## NHANES itself, and eliminate every product that cannot be traced to a parent, 

304 either directly

305 ## or through mapping.

306 ## Then, based on the source of production volume information for each parent, we 

307 select

308 ## the NHANES cycle from which to keep data.

309 ## Make four lists: prodcas, mapprodcas,mapparcas,parcas

310 prodcas <- unique(as.character(NHANES$CAS))

311 mapprodcas <- unique(as.character(mapping$CAS.1))

312 mapparcas <- unique(as.character(mapping$CAS))
313 parcas <- unique(as.character(parent$CASRN))

314

315 ## 1) add chemicals from the intersection of prodcas and parcas to mapping

316

317 toadd <- intersect(prodcas, parcas)

318 indx <- which(parent$CASRN %in% toadd)

319 tmp <- parent[indx,c("Compound","CASRN")]

320 df <- data.frame(Name=tmp$Compound, CAS=tmp$CASRN,MW=rep(NA, nrow(tmp)),

321 Name.1=tmp$Compound, CAS.1=tmp$CASRN, MW.1=rep(NA, nrow(tmp)))

322 indx2 <- which(NHANES$CAS %in% toadd)

323 tmp <- unique(NHANES[indx2,c("CAS","MW")])

324 MW <- tmp$MW

325 names(MW) <- as.character(tmp$CAS)

326 df$MW <- df$MW.1 <- MW[as.character(df$CAS)]

327

328 mapping <- rbind(mapping, df)

329

330 ## 2) eliminate rows in mapping for products not in prodcas

331 ## recreate mapprodcas and mapparcas

332 mapprodcas <- unique(as.character(mapping$CAS.1))

333 mapparcas <- unique(as.character(mapping$CAS))

334

335 ## notinprod is a list of cas numbers that are products in mapping, but not in NHANES

336 notinprod <- mapprodcas[which(!(mapprodcas %in% prodcas))]

337 ## indx is the index of rows of mapping that corresponds to notinprod

338 indx <- which(mapping$CAS.1 %in% notinprod)

339 ## Drop those rows:

340 mapping <- mapping[-indx,]

341 ## rebuild mapprodcas and mapparcas

342 mapprodcas <- unique(as.character(mapping$CAS.1))

343 mapparcas <- unique(as.character(mapping$CAS))

344

345 ## 2) What are the chemicals in prodcas that are not in mapprodcas?

346 notinmapprod <- prodcas[which(!(prodcas %in% mapprodcas))]

347 notinpar <- notinmapprod[which(!(notinmapprod %in% parcas))]

348 ## > NHANES[which(NHANES$CAS == notinpar),]

349 ## Name CAS MW LOD Year Value sample

350 ## 70 9-Hydroxyfluorene 1689-64-1 182.22 5.933235e-08 2005 4.343975e-06 01-02

351 ## gm gmlcl gmucl Value50 Value75 Value90 Value95 LCL50 UCL50 LCL75 UCL75

352 ## 70 0.205 0.183 0.229 0.206 0.353 0.585 0.852 0.181 0.232 0.312 0.393

353 ## LCL90 UCL90 LCL95 UCL95

354 ## 70 0.507 0.697 0.736 1.12

355 ## This is in mapping, but the cas number is wrong. Fix it:

356 indx <- which(mapping$Name.1 == "9-Hydroxyfluorene")

357 mapping$CAS.1[indx] <- "1689-64-1"

358 ## Now rebuild mapprodcas and mapparcas again

359 mapprodcas <- unique(as.character(mapping$CAS.1))

360 mapparcas <- unique(as.character(mapping$CAS))

361 notinmapprod <- prodcas[which(!(prodcas %in% mapprodcas))]

362 ## length(notinmapprod) == 0

363

364 ## 3) Are there parents in mapping that are not in parent?

365 mapparcas <- unique(as.character(mapping$CAS))

366 parcas <- unique(as.character(parent$CASRN))

367 setdiff(mapparcas, parcas)

368 ## > setdiff(mapparcas, parcas)

369 ## [1] "7784-42-1" "42509-83-1"

370 ## These are:

371 mapping[indx <- which(as.character(mapping$CAS) %in% setdiff(mapparcas, parcas)),]

372 ## Name CAS MW Name.1 CAS.1 MW.1

373 ## 81 Arsenic 7784-42-1 77.95 Dimethylarsinic Acid 917-76-0 137.998

374 ## 94 Isazaphos-methyl 42509-83-1 313.74 Dimethylphosphate 813-78-5 126.050

375 ## 109 Isazaphos-methyl 42509-83-1 313.74 Dimethylthiophosphate 1112-38-5 142.110
376 ## Drop them from mapping

377 mapping <- mapping[-indx,]

378 ## Now, the set mapping$CAS.1 should equal the set NHANES$CAS, and

379 ## the set mapping$CAS should equal the set parent$CASRN

380 prodcas <- unique(as.character(NHANES$CAS))

381 mapprodcas <- unique(as.character(mapping$CAS.1))

382 mapparcas <- unique(as.character(mapping$CAS))

383 parcas <- unique(as.character(parent$CASRN))

384 ## Compare parents in the two structures

385 setdiff(parcas, mapparcas)

386 ## > setdiff(parcas, mapparcas)

387 ## character(0)

388 setdiff(mapparcas, parcas)

389 ## > setdiff(mapparcas, parcas)

390 ## character(0)

391 setdiff(mapprodcas, prodcas)

392 ## > setdiff(mapprodcas, prodcas)

393 ## character(0)

394 setdiff(prodcas, mapprodcas)

395 ## [1] "917-76-0"

396 NHANES[which(NHANES$CAS == setdiff(prodcas, mapprodcas)),c("Name","CAS")]

397 ## DMA, from arsenicals. We aren't dealing with arsenic, so drop this chemicals from 

398 NHANES

399 indx <- which(NHANES$CAS %in% setdiff(prodcas, mapprodcas))

400 NHANES <- NHANES[-indx,]

401

402 ## Add PUsource to mapping

403 mapping$PUsource <- as.character(parent[mapping$CAS,"P/U Source"])

404 ## ==================================================================================

405 ## Now, which sample to use? We have the 2002 CPRI, so for those chemicals, lets try 

406 using 01-02,

407 ## and the 2006 IUR, so for those chemicals, use 05-06. First, we need to match 

408 NHANES chemicals

409 ## to the production volume source of their parents. 

410 ## For each product CAS number (mapping$CAS.1) write down which year range to retain.

411 keep <- tapply(mapping$PUsource, factor(mapping$CAS.1),

412 function(x) min(c(`Not IUR HPV`="05-06",`EPA IUR HPV`="05-06",`CPRI 2002`="01-

413 02")[x]))

414

415

416 ## we work on a copy of NHANES, in case we want to use the original later

417 tst <- NHANES

418 tst$keep <- keep[tst$CAS]

419 tst$keep <- factor(tst$keep, levels=levels(tst$sample), ordered=TRUE)

420

421 ## --------------------------------------------------------------------------------

422

423 ## Now, just keep the records for each CAS value that match keep. If there are none,

424 ## keep the most recent earlier record.

425

426 tst2 <- do.call("rbind",

427 unclass(by(tst,tst$CAS,

428 function(x) {

429 indx <- which(x$sample %in% x$keep)

430 if (length(indx) > 0) {

431 ## exact match. Easy

432 z <- x[indx,,drop=FALSE]

433 z$match <- "=="

434 z

435

436 } else {

437 if (length(x$sample) == 1) {

438 ## only 1 row, so just use it.
439 x$match <- "1"

440 x

441 } else {

442 ## Find the latest value in sample earlier than keep

443 z1 <- factor(as.character(x$sample),

444 levels=levels(tst$sample), ordered=TRUE)

445 kp <- unique(x$keep)

446 if (any(z1 < kp)) {

447 indx <- which(z1 == max(z1[z1<unique(x$keep)]))

448 z <- x[indx,,drop=FALSE]

449 z$match <- "<"

450 z

451 } else {

452 z <- x[which.min(z1),,drop=FALSE]

453 z$match <- ">"

454 z

455 }

456 }

457 }

458 }

459 )))

460 Measured <- tst2

461

462 ## --------------------------------------------------------------------------

463 ## Create pred.data

464

465 pred.data <-

466 data.frame(CAS=parent$CASRN,

467 Name=parent$Compound,

468 s.U = parent$`USEtox Predicted mg/kg BW per kg/day released`,

469 s.R = parent$`RAIDAR Predicted mg/kg BW per kg/day released`,

470 V = V <- parent$`Production/Usage (kg/day)`,

471 PUsource = Z <- parent$`P/U Source`,

472 V.cpri = ifelse(Z == "CPRI 2002",V,NA))

473

474 V.cat <- character(nrow(pred.data))

475 for (i in 1:nrow(pred.data)) {

476 if (Z[i] == "CPRI 2002") {

477 V.cat[i] <- NA

478 } else {

479 V.cat[i] <- as.character(geometric.mean.to.IUR(pred.data$V[i])[,1])

480 }

481 }

482 levels <- c("< 25,000 lbs",

483 "25,000 to < 500,000 lbs",

484 "500,000 to < 1 million lbs",

485 "1 to < 10 million lbs",

486 "10 to < 50 million lbs",

487 "50 to < 100 million lbs",

488 "100 to < 500 million lbs",

489 "500 million to < 1 billion lbs",

490 "1 billion lbs and greater")

491 pred.data$V.cat <- factor(V.cat, levels=levels)

492

493 ## Create mapping

494

495 ## We need this object as it is here in a little while, so

496 ## save it

497 oldmap <- mapping

498

499 ## Simplify the mapping object to just contain Name, CAS, Name.1, and CAS.1

500 mapping <- subset(mapping,

501 select=c("Name","CAS","Name.1","CAS.1"))
502

503 ## -----------------------------------------------------------------------

504 ## Convert measurements

505

506 ## Convert all the Values and CLs using convert.urine.to.exposure()

507

508 NHANES$gm <- convert.urine.to.exposure(NHANES$gm)

509 NHANES$gmlcl <- convert.urine.to.exposure(NHANES$gmlcl)

510 NHANES$gmucl <- convert.urine.to.exposure(NHANES$gmucl)

511 NHANES$Value50 <- convert.urine.to.exposure(NHANES$Value50)

512 NHANES$Value75 <- convert.urine.to.exposure(NHANES$Value75)

513 NHANES$Value90 <- convert.urine.to.exposure(NHANES$Value90)

514 NHANES$Value95 <- convert.urine.to.exposure(NHANES$Value95)

515 NHANES$LCL50 <- convert.urine.to.exposure(NHANES$LCL50)

516 NHANES$LCL75 <- convert.urine.to.exposure(NHANES$LCL75)

517 NHANES$LCL90 <- convert.urine.to.exposure(NHANES$LCL90)

518 NHANES$LCL95 <- convert.urine.to.exposure(NHANES$LCL95)

519 NHANES$UCL50 <- convert.urine.to.exposure(NHANES$UCL50)

520 NHANES$UCL75 <- convert.urine.to.exposure(NHANES$UCL75)

521 NHANES$UCL90 <- convert.urine.to.exposure(NHANES$UCL90)

522 NHANES$UCL95 <- convert.urine.to.exposure(NHANES$UCL95)

523

524 ## Need to get the appropriate LOQs (limits of quantification).

525 ## I have extracted the text for the 2005 table and formatted

526 ## it to read as a space-delimited file. Read that in,

527 ## and match to the value of 'sample' in Measured. This value

528 ## is in ug/L; convert it to ug/g creatinine by dividing

529 ## by average g creatinine / L urine

530

531 LOQ2005 <- read.table("LOQ2005.txt",header=TRUE,

532 

533 colClasses=c("character","character","character","numeric","numeric"),

534 sep="|")

535 LOQ2011 <- read.table("LOQ2011.txt",header=TRUE,

536 colClasses=c("character","character","character",

537 "numeric","numeric","numeric","numeric","numeric"),

538 sep="|")

539 LOQ <- numeric(nrow(Measured))

540 for (i in 1:nrow(Measured)) {

541 name <- Measured$Name[i]

542 if (Measured$Year[i] == 2005) {

543 name <- gsub("\\[","\\\\[",gsub("]","\\\\]",name))

544 name <- gsub("\\(","\\\\(",gsub(")","\\\\)",name))

545 name <- paste("^",name,sep="")

546 indx <- grep(name, LOQ2005$Name, ignore.case=TRUE)

547 if (length(indx) == 0) stop(paste("Cannot find ",Measured$Name[i]," in 

548 LOQ2005",sep=""))

549 mult <- switch(LOQ2005$Units[indx],"ug/L" = 1,"ug/dL" = 10,"ng/mL" = 1, "ng/L" = 

550 0.001)

551 ## Finals in LOQ * mult are ug/L. Divide by mean creatinine in g/L to get

552 ## LOQ in ug / g creatinine

553 ## In 99-00, that was 1.364 g/L (creatinine / urine) (page 5 in NHANES 2005

554 ## document

555 ## in 01-02 that was 1.306

556 LOQ[i] <-

557 switch(as.character(Measured$sample[i]),

558 "99-00" = LOQ2005[indx,"Y1999.2000"] / 1.364,

559 "01-02"=LOQ2005[indx,"Y2001.2002"]/1.306)*mult

560 } else {

561 name <- gsub("\\[","\\\\[",gsub("]","\\\\]",name))

562 name <- gsub("\\(","\\\\(",gsub(")","\\\\)",name))

563 name <- gsub(" "," \\*",name)

564 name <- paste("^",name,sep="")
565 indx <- grep(name, LOQ2011$Name, ignore.case=TRUE)

566 if (length(indx) == 0) 

567 stop(paste("Cannot find ",Measured$Name[i]," in LOQ2011",sep=""))

568 mult <- switch(LOQ2011$Units[indx],"ug/L" = 1,"ng/mL" = 1, "pg/mL" = 0.001)

569 ## No easily found value for creatinine / urine in NHANES IV, so use the 01-02

570 ## value from NHANES III

571 LOQ[i] <-

572 switch(as.character(Measured$sample[i]),"99-00"=LOQ2011[indx,"Y1999.2000"],

573 "01-02"=LOQ2011[indx,"Y2001.2002"],

574 "03-04"=LOQ2011[indx,"Y2003.2004"],

575 "05-06"=LOQ2011[indx,"Y2005.2006"],

576 "07-08"=LOQ2011[indx,"Y2007.2008"])*mult / 1.306

577 }

578 }

579 ## LOQ is in ug/g creatinine.

580

581 Measured$LOQ <- convert.urine.to.exposure(LOQ)

582 Measured$gm <- convert.urine.to.exposure(Measured$gm)

583 Measured$gmlcl <- convert.urine.to.exposure(Measured$gmlcl)

584 Measured$gmucl <- convert.urine.to.exposure(Measured$gmucl)

585 Measured$Value50 <- convert.urine.to.exposure(Measured$Value50)

586 Measured$Value75 <- convert.urine.to.exposure(Measured$Value75)

587 Measured$Value90 <- convert.urine.to.exposure(Measured$Value90)

588 Measured$Value95 <- convert.urine.to.exposure(Measured$Value95)

589 Measured$LCL50 <- convert.urine.to.exposure(Measured$LCL50)

590 Measured$LCL75 <- convert.urine.to.exposure(Measured$LCL75)

591 Measured$LCL90 <- convert.urine.to.exposure(Measured$LCL90)

592 Measured$LCL95 <- convert.urine.to.exposure(Measured$LCL95)

593 Measured$UCL50 <- convert.urine.to.exposure(Measured$UCL50)

594 Measured$UCL75 <- convert.urine.to.exposure(Measured$UCL75)

595 Measured$UCL90 <- convert.urine.to.exposure(Measured$UCL90)

596 Measured$UCL95 <- convert.urine.to.exposure(Measured$UCL95)

597 ## Now, measurements are mg/kg bodyweight/day

598

599 ## ---------------------------------------------------------------

600 ## Build a matrix with rows for parents and columns for metabolites

601 ## that is, rows for the records in pred.data, and columns for the

602 ## records in Measured. use CAS numbers to do the linking.

603 ## The entries will be NA for values that will be estimated (that sum to 

604 ## 1 across rows),

605 ## 1 for entries representing 1:1 molar relationships between parent and metabolite,

606 ## and 0 everywhere else (the majority of entries).

607 ## Set all the entries to 0, initially.

608 ## Then, set the cells corresponding to entries in NHANES.parent.metabolite.mapping

609 ## to NA. Some of these will be reset to 1: the ones identified by Jim R, 

610 ## listed below,

611 ## and the directly observed parents. These are the CAS numbers in pred.data that

612 ## are also in NHANES.

613 Xref <- matrix(0, nrow=nrow(pred.data), ncol=nrow(Measured),

614 dimnames=list(pred.data$CAS, Measured$CAS))

615 indx <- cbind(match(mapping$CAS,rownames(Xref)),

616 match(mapping$CAS.1,colnames(Xref)))

617

618 Xref[indx] <- NA

619 ## Any row with only 1 metabolite should have the corresponding entry set to 1

620 Nmets <- apply(Xref,1,function(z) sum(is.na(z)))

621 ## This is probably crude

622 for (i in seq(along=Nmets)) {

623 if (Nmets[i] == 1) Xref[i, which(is.na(Xref[i,]))] <- 1

624 }

625 ##dirobs <- pred.data$CAS[pred.data$CAS %in% NHANES$CAS]

626 ##Xref[cbind(match(dirobs, rownames(Xref)), match(dirobs, colnames(Xref)))] <- 1

627 ## A few parent-metabolite pairs are known to be 1-1.
628 fixed <- cbind(Parent=c("2921-88-2","68359-37-5","52918-63-5","52918-63-5",

629 "333-41-5","298-00-0","52315-07-8","52645-53-1"),

630 Met=c("6515-38-4","77279-89-1","3739-38-6","53179-78-5",

631 "2814-20-2","100-02-7","3739-38-6","3739-38-6"))

632 ## Set the entries in Xref that correspond to these rows and columns

633 ## to 1

634 indx <- cbind(match(fixed[,1],rownames(Xref)),match(fixed[,2], colnames(Xref)))

635 Xref[indx] <- 1

636

637 ## We want the parents with multiple metabolites at the top of the list, keeping

638 ## the metabolites in the same order as before. First, how many metabolites for each

639 ## parent?

640 N.metabs <- apply(Xref,1, function(z) sum(is.na(z)))

641 indx <- order(N.metabs, decreasing=TRUE)

642 N.metabs <- N.metabs[indx]

643 ## Now, the names for N.metabs are used to order pred.data. First,

644 ## assign useful rownames to pred.data:

645 rownames(pred.data) <- as.character(pred.data$CAS)

646 tmp <- pred.data[names(N.metabs),]

647 pred.data <- tmp

648 Xref <- Xref[rownames(pred.data),]

649

650 ## +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

651 ## Up to here, all the mass measures are in g or mg. convert unit exposures

652 ## and metabolite (product) concentrations to moles

653 ## John says the MW values in NHANES.parent.metabolite.mapping are the most

654 ## reliable, so use those first, then go to the NHANES object for the

655 ## rest.

656

657 ## Build a reference list of molecular weights indexed by CAS number

658 ParentMW <- unique(subset(oldmap, select=c("CAS","MW")))

659

660 MetaboliteMW <- unique(subset(oldmap, select=c("CAS.1","MW.1")))

661 ## One duplicated metabolite CAS, MW difference less than 0.00005

662 MetaboliteMW <- MetaboliteMW[-which(duplicated(MetaboliteMW$CAS.1)),]

663 names(MetaboliteMW) <- c("CAS","MW")

664

665 ## Combine these two. What are missing? 

666

667 MW <- unique(rbind(ParentMW, MetaboliteMW))

668 rownames(MW) <- as.character(MW$CAS)

669

670 missingpredCAS <- as.character(pred.data$CAS[which(!(pred.data$CAS %in% 

671 rownames(MW)))])

672

673 MW <- unique(rbind(MW, subset(NHANES, CAS %in% missingpredCAS, select=c("CAS","MW"))))

674 rownames(MW) <- as.character(MW$CAS)

675

676

677

678 ## Need to convert s.U, and s.R to use molar units.

679 ## construct the MW vector to match pred.data

680 pred.data$MW <- MW[as.character(pred.data$CAS),"MW"]

681 pred.data$s.U <- with(pred.data, s.U / (MW * 1000))

682 pred.data$s.R <- with(pred.data, s.R / (MW * 1000))

683

684 ## Now, s.U and s.R are moles / kg produced / kg bodyweight

685

686 Measured$Value <- with(Measured, Value /(MW * 1000))

687 Measured$LOD <- with(Measured, LOD / (MW * 1000))

688

689 Measured$LOQ <- with(Measured, LOQ / (MW * 1000))

690 Measured$gm <- with(Measured, gm / (MW * 1000))
691 Measured$gmlcl <- with(Measured, gmlcl / (MW * 1000))

692 Measured$gmucl <- with(Measured, gmucl / (MW * 1000))

693 Measured$Value50 <- with(Measured, Value50 /(MW * 1000))

694 Measured$Value75 <- with(Measured, Value75 /(MW * 1000))

695 Measured$Value90 <- with(Measured, Value90 /(MW * 1000))

696 Measured$Value95 <- with(Measured, Value95 /(MW * 1000))

697 Measured$LCL50 <- with(Measured, LCL50 /(MW * 1000))

698 Measured$LCL75 <- with(Measured, LCL75 /(MW * 1000))

699 Measured$LCL90 <- with(Measured, LCL90 /(MW * 1000))

700 Measured$LCL95 <- with(Measured, LCL95 /(MW * 1000))

701 Measured$UCL50 <- with(Measured, UCL50 /(MW * 1000))

702 Measured$UCL75 <- with(Measured, UCL75 /(MW * 1000))

703 Measured$UCL90 <- with(Measured, UCL90 /(MW * 1000))

704 Measured$UCL95 <- with(Measured, UCL95 /(MW * 1000))

705

706 ## Measurements are moles / kg bodyweight / day

707 ## These are small numbers, so multiply by 1e9 to get

708 ## nmoles / kg bodyweight / day

709 pred.data[,c("s.U","s.R")] <- 1e9 * pred.data[,c("s.U","s.R")]

710

711 Measured[,c("Value","LOD","LOQ","gm","gmlcl","gmucl",

712 "Value50","Value75","Value90","Value95",

713 "LCL50","LCL75","LCL90","LCL95",

714 "UCL50","UCL75","UCL90","UCL95")] <-

715 1e9 * Measured[,c("Value","LOD","LOQ","gm","gmlcl","gmucl",

716 "Value50","Value75","Value90","Value95",

717 "LCL50","LCL75","LCL90","LCL95",

718 "UCL50","UCL75","UCL90","UCL95")]

719 ## Now, values in 'Measured' are nmoles / kg bodyweight / day

720

721 cutpoints <- c(6250,25000,5e5,1e6,10e6,50e6,100e6,500e6,1000e6,1e37) ## lbs / year

722 ## Convert cutpoints to kg / day

723 cutpoints <- cutpoints * 0.45359237 / 365.25

724

725 ## Standard errors for Measured. If all values are above LOD, then sexx is

726 ## (log(UCLxx) - log(LCLxx))/(2 * 1.96) else (log(UCLxx) - log(Valuexx))/1.96

727

728 ## Median

729 se2 <- with(Measured, ((log(UCL50) - log(LCL50)) / (2 * 1.96)))

730 se1 <- with(Measured, (log(UCL50) - log(Value50)) / 1.96)

731 Measured$se50 <- ifelse(is.na(se2), se1,se2)

732

733 se2 <- with(Measured, ((log(UCL75) - log(LCL75)) / (2 * 1.96)))

734 se1 <- with(Measured, (log(UCL75) - log(Value75)) / 1.96)

735 Measured$se75 <- ifelse(is.na(se2), se1,se2)

736

737 se2 <- with(Measured, ((log(UCL90) - log(LCL90)) / (2 * 1.96)))

738 se1 <- with(Measured, (log(UCL90) - log(Value90)) / 1.96)

739 Measured$se90 <- ifelse(is.na(se2), se1,se2)

740

741 se2 <- with(Measured, ((log(UCL95) - log(LCL95)) / (2 * 1.96)))

742 se1 <- with(Measured, (log(UCL95) - log(Value95)) / 1.96)

743 Measured$se95 <- ifelse(is.na(se2), se1,se2)

744

745 ## Now, what about missing values? Let's impute them

746 ## Build a data frame.

747 Z <- Measured[,c("se50","se75","se90","se95")]

748 Zd <- reshape(Z,direction="long",v.names="se", varying=names(Z),idvar="chem", 

749 timevar="Quantile")

750 out <- lm(se ~ factor(chem) + factor(Quantile), data=Zd, na.action=na.exclude)

751 out2 <- lm(se ~ factor(Quantile), data=Zd, na.action=na.exclude)

752 out3 <- lm(se ~ factor(chem), data=Zd, na.action=na.exclude)

753 ## Now, get the predictions for all the values in Zd
754 anova(out,out2)

755 anova(out,out3)

756

757 ## No real differences among quantiles, but significant differences among chems.

758 ## So, set a single se for each chemical, the sqrt of the average of the quant759 ## specific variances.

760

761 se <- apply(data.matrix(Z),1,function(x) {

762 if (!all(is.na(x))) {

763 sqrt(mean(x*x, na.rm=TRUE))

764 } else NA

765 })

766

767 ## Replace the NAs with the overall mean of the non-missing ses

768

769 se[is.na(se)] <- mean(se, na.rm=TRUE)

770

771 Measured$se <- se

772 attach("UseData-041013.RData")

773 ## Near field: Fragrance, food additive, consumer use, personal care product.

774 indx <- match(pred.data$CAS, AliciaUseTable$CASRN)

775 ##indx2 <- match(pred.data$CAS, UseData$CASRN)

776 Uses <- AliciaUseTable[indx,]

777 ##Uses2 <- UseData[indx2,]

778 Uses$NearField <- as.numeric(with(Uses, FRAGRANCE | FOOD.ADDITIVE | CONSUMER.USE | 

779 PERSONAL.CARE.PRODUCT | PHARMACEUTICAL))

780 pred.data$NearField <- Uses$NearField

781 ## One missing classification, for benzo[c]phenanthrene cas 195-19-7

782 ## This is a component of tobacco smoke, and is near field

783 pred.data$NearField[which(pred.data$CAS == "195-19-7")] <- 1

784 save(Measured, cutpoints, pred.data, N.metabs,Xref, Uses, 

785 file=file.path("..","MCMCdata.RData"))

786

787 ## Counts for uses in the NHANES data:

788

789 print(matrix(apply(data.matrix(Uses[,-c(1,2,3)]),2,sum),

790 ncol=1,

791 dimnames=list(colnames(Uses)[-c(1,2,3)],"Count")))

792

793 writeLines(c("","Breakdown for Near Field:",

794 "Codes are 0 for absence, 1 for presence of the classification",

795 "Order is: Consumer use, food additive, fragrance, personal care product, 

796 pharmaceutical",

797 "so 10100 is classified consumer use and fragrance"))

798 tmp <- table(with(Uses,

799 

800 paste(CONSUMER.USE,FOOD.ADDITIVE,FRAGRANCE,PERSONAL.CARE.PRODUCT,PHARMACEUTICAL,sep=""

801 )))

802 matrix(tmp, ncol=1, dimnames=list(names(tmp), "Count"))

803 writeLines(c("",

804 "Breakdown for the whole dataset:",

805 "Order is: consumer use, food additive, fragrance, personal care 

806 product,",

807 "antimicrobial,chemical industrial process, colorant,",

808 "fertilizer, herbicide, pesticide, pharmaceutical, toxin,inert, flame 

809 retardant"))

810 tmp <- table(with(Uses,

811 paste(CONSUMER.USE,FOOD.ADDITIVE,FRAGRANCE,PERSONAL.CARE.PRODUCT,

812 ANTIMICROBIAL,CHEMICAL.INDUSTRIAL.PROCESS,COLORANT,

813 FERTILIZER,HERBICIDE,

814 PESTICIDE,PHARMACEUTICAL,TOXIN,INERT,FLAME.RETARDANT,sep="")))

815

816 matrix(tmp, ncol=1, dimnames=list(names(tmp), "Count"))
OnlyP/model0.jag: JAGS Code to Estimate Parent Chemical Exposures 

from NHANES

1 model

2 {

3 ## +++++++++++++++++++++++++++++++++++++++++++++

4 ## The heart of the model:

5 ## Model the parent exposures

6 ## +++++++++++++++++++++++++++++++++++++++++++++

7 for (i in 1:N) {

8 ## log of product of production volume and unit exposure

9 lP[i] ~ dnorm(lPmu, tau.V)

10 P[i]<- exp(lP[i])

11 }

12 ## hyperparameter for lP

13 lPmu ~ dnorm(0, 0.001)

14

15 ## code so that prior for sd(log(P[i])) is half-Cauchy(25)

16 sd.dum ~ dt(0,1/625,1)

17 sd.V <- abs(sd.dum)

18 tau.V <- 1/pow(sd.V,2)

19 ## ++++++++++++++++++++++++++++++++++++++++++++++

20 ## Link the unobserved parent exposure to the observed

21 ## (but censored) metabolite (w/some parent) urinary

22 ## values.

23 ##

24 ## Model metabolite concentrations as (usually) unknown fractions

25 ## of the parent.

26 ## 

27 ## Phi[i,j] say what fraction of parent[i] goes to metabolite[j]

28 ## observed values are lognormally distributed around U, with known

29 ## precision, but censored. Some of the elements of Phi are known to be

30 ## 0 or 1 (thanks to Jim Rabinowitz for some of the 1's).

31 for (j in 1:M) {

32 lU[j]<- log(inprod(Phi[,j],P[]))

33 for (i in 1:4) {

34 ly[i,j] ~ dnorm(qnorm(Q[i],lU[j], tau[j]),tau.se[j])

35 is.above.LOD[i,j] ~ dinterval(ly[i,j],llod[j])

36 }

37 tau.se[j] <- 1/(se[j]*se[j])

38 ## lognormal distribution for chemical-specific population sds

39 sd[j] ~ dlnorm(lsd, tau.lsd)

40 tau[j] <- 1/(sd[j]*sd[j])

41 }

42 lsd ~ dnorm(0, 0.001)

43 lsd.dum ~ dt(0,1/625,1)

44 lsd.sd <- abs(lsd.dum)

45 tau.lsd <- 1/(lsd.sd * lsd.sd)

46 ## P is sorted so that only the first N.nontriv parents result in multiple

47 ## metabolites. For i in (N.nontriv+1:N), there is only one metabolite.

48 ## Nmets[i] gives the number of metabolites of parent i

49 ## indx[i,1:Nmets[i]] give the columns of row i that are being estimated,

50 ## for i in 1:N.nontriv. When i > N.nontriv, indx[i,1] gives the 1

51 ## metabolite matched to parent i.

52 ##

53 ## The elements of phi are essentially unknown. 

54 for (k in 1:Ndelta) {

55 delta[k] ~ dunif(0,1)

56 }

57 for (i in 1:N.nontriv) {

58 for (j in 1:Nmets[i]) {
59 phi[dstart[i] + j - 1] <- delta[dstart[i] + j - 1] / 

60 sum(delta[dstart[i]:dstop[i]])

61 Phi[i,indx[i,j]]<- phi[dstart[i] + j - 1]

62 }

63 }

64 }
OnlyP/fitmodel0.R: R Code to Estimate Parent Chemical Exposure from 

NHANES

1 library(doMC)

2 registerDoMC(cores=3)

3 library(rjags)

4 load.module("lecuyer")

5 library(random)

6 library(MASS) ## for ginv, Moore-Penrose inverse

7

8 doONEchain <- TRUE

9

10 ## JAGS model is in model0.jag

11 ## It needs the following data:

12 ## y : observed metabolite rates (NA for below LOD) (1:M),

13 ## N = number of parents,

14 ## N.nontriv = number of parents with > 1 metabolite (sorted so they come first

15 ## in R.x)

16 ## N.mets[1:N.nontriv] = number of metabolites for parent i

17 ## indx[N,M]: indx[i,1:Nmets[i]] gives the indices in U for the metabolites

18 ## i contributes. Corresponding values of Phi[i,] are non-zero.

19 ## If i > N.nontriv, indx[i,1] gives the entry in Phi[i,] that is 1

20 ## For j > N.nontriv, indx[i,j] gives column indices in Phi[i,] for 0's.

21 ## M = number of metabolites,

22 ## tau[1:M] : precision for observed metabolite rates.

23 ## lod = vector of limits of detection, tau = vector of precisions for observed

24 ## metabolite rates.

25 ## dstart[1:N.nontriv], dstop[1:N.nontriv]: indices into delta and phi

26 ## such that, e.g., delta[dstart[i]:dstop[i]] are the deltas for row i of Phi.

27 ## dstart[1] == 1, dstop[N.nontriv] == Ndelta.

28 ##

29 ## The following variables need initializing

30 ## P[i], i in 1:N

31 ## X1

32 ## X2

33 ## delta[1:Ndelta]: the non-trivial deltas corresponding to the non-trivial phis

34

35 ## Data are based on John W.'s RData file DataTables-040212.RData, and successors.

36

37 attach(file.path("..","MCMCdata.RData"))

38

39 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

40 ## JAGS inputs: nhanesdata and inits

41 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

42 if (doONEchain) {

43 nhanesdata <- list()

44 ## The Jags models need these data:

45 ## IURcuts - cutpoints for the IUR production volumes (cutpoints above)

46 nhanesdata$IURcuts <- log(cutpoints)

47 ## IURindx - vector indexing the rows of pred.data with IUR production volumes

48 nhanesdata$IURindx <- grep("IUR HPV",pred.data$PUsource)

49 ## IURvols - indicator for which interval an IUR production volume falls

50 nhanesdata$IURvols <-

51 as.numeric(pred.data$V.cat[nhanesdata$IURindx[1:length(nhanesdata$IURindx)]])

52 ## CPRIindx indexes the rows of pred.data with CPRI use volumes

53 nhanesdata$CPRIindx <- which(pred.data$PUsource == "CPRI 2002")

54 nhanesdata$lV.CPRI <- log(pred.data[nhanesdata$CPRIindx,"V"])

55 ## ls.R, ls.U are the log unit predictions from RADAR and USEtox, respectively

56 nhanesdata$ls.R <- log(pred.data$s.R)

57 nhanesdata$ls.U <- log(pred.data$s.U)

58 nhanesdata$NearField <- pred.data$NearField
59 ## Mean of ls.R and ls.U, for centering

60 nhanesdata$meanls.R <- mean(nhanesdata$ls.R)

61 nhanesdata$meanls.U <- mean(nhanesdata$ls.U)

62

63 ## ly - the log of the observed metabolite quantiles

64 nhanesdata$ly <-

65 ly <- t(log(data.matrix(Measured[,c("Value50","Value75","Value90","Value95")])))

66 ## Q, the proportions associated with the quantiles

67 nhanesdata$Q <- Q <- c(0.5,0.75,0.9,0.95)

68 ## se - precision of measured metabolite values, from reported confidence limits on 

69 quantiles

70 nhanesdata$se <- Measured$se

71 ## llod - log(Limit of Detection)

72 nhanesdata$llod <- log(Measured$LOQ)

73 ## is.above.LOD - vector with 0 for below lod, 1 for above it

74 nhanesdata$is.above.LOD <- 0 + !is.na(ly)

75 ## Some counts:

76 ## N is the number of parents

77 nhanesdata$N <- nrow(pred.data)

78 ## M is the number of metabolites

79 nhanesdata$M <- nrow(Measured)

80 ## N.IUR is the number of IUR reported production volumes (length of IURindx)

81 nhanesdata$N.IUR <- length(nhanesdata$IURindx)

82 ## N.CPRI is the number of CPRI reported use volumes (length of CPRIindx)

83 nhanesdata$N.CPRI <- length(nhanesdata$CPRIindx)

84 ## N.nontriv is the number of parents that result in >1 metabolite

85 nhanesdata$N.nontriv <- sum(N.metabs > 1)

86 ## Nmets - number of metabolites for parent[i], when > 1

87 nhanesdata$Nmets <- N.metabs[N.metabs > 1]

88 ## index - index[i,j] is the index into the metabolite list for the jth metabolite

89 ## of the ith parent

90 nhanesdata$indx <-

91 t(apply(Xref,1,function(x)c((1:length(x))[is.na(x)],

92 (1:length(x))[!is.na(x)])))

93 ## begin and end index into phi[] for each row of Phi

94 tmp <- cumsum(nhanesdata$Nmets)

95 nhanesdata$dstart <- 1+c(0,tmp[-nhanesdata$N.nontriv])

96 nhanesdata$dstop <- tmp

97

98 ## Phi

99 nhanesdata$Phi <- Xref

100 nhanesdata$Ndelta <- sum(is.na(Xref))

101

102 ## ---------------------------------------------------------------------

103 ## Initialization

104

105 inits <- list(list())

106 ## inits$lP: get estimates for lU, then work backwards to lP

107 ## assuming equal distribution of parent to metabolites when unknown

108 lU <- inits[[1]]$sd <- numeric(ncol(ly))

109 for (j in 1:ncol(ly)) {

110 ## if there are at least 2 values > LOQ, then estimate lmu and sd using least 

111 squares

112 ## against the quantiles. Otherwise, set lmu to the lLOQ - 1, sd to 1.

113 z <- ly[,j]

114 if (sum(!is.na(z)) >= 2) {

115 keep <- !is.na(z)

116 x <- Q[keep]

117 y <- z[keep]

118 pstart <- c(min(y)-1,0)

119 obfun <- function(p) {

120 sum((y - qnorm(x,p[1],exp(p[2])))^2)

121 }
122 out <- try(optim(pstart, obfun))

123 if (!is(out, "try-error")) {

124 lU[j] <- out$par[1]

125 inits[[1]]$sd[j] <- exp(out$par[2])

126 } else {

127 lU[j] <- min(y)

128 inits[[1]]$sd[j] <- 1

129 }

130 } else {

131 lU[j] <- nhanesdata$llod[j] - 1

132 inits[[1]]$sd[j] <- 1

133 }

134 }

135 XX <- Xref

136 for (i in 1:nrow(XX)) {

137 nNA <- sum(is.na(XX[i,]))

138 if (nNA > 0) {

139 XX[i,is.na(XX[i,])] <- 1/nNA

140 }

141 }

142 XXt <- t(XX)

143 meanlU <- mean(lU)

144 Usc <- exp(lU - meanlU)

145 obfun <- function(p) {

146 sum((XXt %*% exp(p) - Usc)^2)

147 }

148 pstart <- as.vector(ginv(XXt) %*% Usc)

149 pstart[pstart< 0] <- .1

150 pstart <- log(pstart)

151 out <- optim(pstart, obfun, control=list(maxit=500000))

152

153 inits[[1]]$lP <- out$par + meanlU 

154

155 inits[[1]]$sd.dum <- sd(inits[[1]]$lP)

156

157 ## for the censored values in ly

158 ily <- is.na(ly)

159 loqmx <- outer(rep(1,4),log(Measured$LOQ),"*")

160 lyinit <- ifelse(ily,loqmx-2,NA)

161 inits[[1]]$ly <- lyinit

162

163 inits[[1]]$lsd <- mean(log(inits[[1]]$sd))

164 inits[[1]]$lsd.dum <- sd(log(inits[[1]]$sd))

165

166 inits[[1]]$delta <- rep(1, nhanesdata$Ndelta)

167 inits[[1]]$lPmu <- mean(inits[[1]]$lP)

168

169 doinits <- function(chain) {inits[[chain]]}

170 assign("inits", inits, envir=environment(doinits))

171

172 ## -------------------------------------------------

173 ## Run JAGS

174

175 ## Run one chain, keeping N.Iters samples after thinning every N.Thin sample

176 ## Burnin for N.Burnin samples, which are discarded.

177 ## These samples, which are at least closed to converged, serve to provide

178 ## initial values for the three-chain run, below.

179 N.Iters <- 1000

180 N.Thin <- 500

181 N.Burnin <- 1000

182 model.jags <- jags.model("model0.jag",data=nhanesdata, inits=doinits,

183 n.chains=1,n.adapt=N.Burnin)

184
185 out.samples <- coda.samples(model.jags,

186 variable.names=c("lP","phi","lU","lPmu","sd.V",

187 "sd","lsd","lsd.sd"),

188 n.iter=N.Iters * N.Thin, thin=N.Thin)

189 

190 nms <- varnames(out.samples)

191 for (nm in nms) {

192 writeLines(c("","-----------------",nm,""))

193 print(try(heidel.diag(out.samples[,nm])))

194 print(try(summary(out.samples[,nm])))

195 }

196

197 save(out.samples, nhanesdata, inits, file="m0parms1.RData")

198 }

199 ## --------------------------------------------------------------------------------

200 ## This is the initial, hopefully converged, chain. Now, use these to start with

201 ## and run three new chains.

202

203 ## The Gelman-Rubin diagnostic works best if the initial values for the multiple

204 ## chains are hyperdispersed relative to the posterior distribution. getInits

205 ## accomplishes that by first picking ninits random directions in parameter space, 

206 ## projecting the scaled samples from the posterior from the one-chain sample passed 

207 ## as X onto those directions, and picking the Pth quantile of those projections. The 

208 ## return value is the row index of the resulting sample in the input X.

209 getInits <- function(X, ninits=3, P=0.05) {

210 ## X is an mcmc object, essentially a matrix.

211 dta <- scale(X)

212 ## Ultimately, we are looking for ninits indices into the elements of X

213 outinits <- numeric(ninits)

214 ## Pick ninits random directions. Project the data in dta on each direction

215 ## and pick the Pth quantile (since the original direction u was random

216 ## u and -u are equally likely).

217 for (i in 1:ninits) {

218 u <- rnorm(ncol(dta))

219 u <- u/sqrt(sum(u*u))

220 lngt <- dta %*% u

221 indx <- seq(along=lngt)

222 keep <- round(P*length(indx)) ## this works fine as long as length(indx) > 20

223 oind <- order(lngt)

224 outinits[i] <- oind[keep]

225 }

226 outinits

227 }

228 ## Get 3 hyperdispersed points.

229 out.coda <- out.samples[[1]]

230

231 Indx <- getInits(out.coda,3)

232

233 ## Each chemical gets a value for a lot of the parameters, like "lP". JAGS

234 ## labels those like elements of a vector, e.g., lP[1]. getsamps

235 ## returns all the variables for a given parameter.

236 getsamps <- function(nm) {

237 nm <- gsub("\\.","\\\\.",nm)

238 cindx <- grep(nm, colnames(out.coda))

239 out.coda[,cindx,drop=FALSE]

240 }

241

242 ## We run three chains.

243 inits2 <- list(list(),list(),list())

244 for (i in 1:3) {

245 inits2[[i]]$lP <- getsamps("^lP\\[")[Indx[i],]

246 inits2[[i]]$sd.dum <- getsamps("^sd.V")[Indx[i],]

247 inits2[[i]]$sd <- getsamps("^sd\\[")[Indx[i],]
248 inits2[[i]]$lsd <- getsamps("^lsd$")[Indx[i],]

249 inits2[[i]]$lsd.dum <- getsamps("lsd.sd")[Indx[i],]

250 inits2[[i]]$delta <- getsamps("phi")[Indx[i],]

251 inits2[[i]]$ly<-inits[[1]]$ly

252 inits2[[i]]$lPmu <- getsamps("lPmu")[Indx[i],]

253 inits2[[i]]$.RNG.name="lecuyer::RngStream"

254 inits2[[i]]$.RNG.seed=randomNumbers(n=1, min=1, max=1e+06, col=1)

255 }

256

257 doinits <- function(chain){inits2[[chain]]}

258 env1 <- new.env()

259 assign("inits2", inits2, pos=env1)

260 environment(doinits) <- env1

261

262 NIters <- 1000

263 Thin <- 500

264 NBurnin <- 5000

265

266 out.samps3R <- foreach(i=1:3) %dopar% {

267 model <- jags.model("model0.jag", data=nhanesdata, inits=inits2[i], n.chains=1,

268 n.adapt=NBurnin)

269 result <- coda.samples(model, variable.names=c("lP","phi","lU","lPmu","sd.V",

270 "sd","lsd","lsd.sd"),

271 n.iter=NIters*Thin, thin=Thin)

272 return(result)

273 }

274

275 out.coda3R <- do.call("mcmc.list", lapply(out.samps3R, function(z) z[[1]]))

276 save(out.samps3R, out.coda3R, file="m3parms3.RData")

277 nms <- varnames(out.coda3R)

278 for (nm in nms) {

279 writeLines(c("","-----------------",nm,""))

280 print(try(heidel.diag(out.coda3R[,nm, drop=FALSE])))

281 print(try(gelman.diag(out.coda3R[,nm,drop=FALSE])))

282 print(summary(out.coda3R[,nm,drop=FALSE]))

283 }
OnlyV/fitmodel0.R: R Code to Quantify Uncertainty in Production 

Volume

1 library(rjags)

2 library(MASS) ## for ginv, Moore-Penrose inverse

3 ## Compare RAIDAR and USEtox to NHANES urinary data

4 ## JAGS model is in nhanes2.jag

5 ## It needs the following data:

6 ## R.R : RAIDAR predicted exposure (1:N),

7 ## R.U : USEtox predicted exposure (1:N),

8 ## y : observed metabolite rates (NA for below LOD) (1:M),

9 ## N = number of parents,

10 ## N.nontriv = number of parents with > 1 metabolite (sorted so they come first

11 ## in R.x)

12 ## N.mets[1:N.nontriv] = number of metabolites for parent i

13 ## indx[N,M]: indx[i,1:Nmets[i]] gives the indices in U for the metabolites

14 ## i contributes. Corresponding values of Phi[i,] are non-zero.

15 ## If i > N.nontriv, indx[i,1] gives the entry in Phi[i,] that is 1

16 ## For j > N.nontriv, indx[i,j] gives column indices in Phi[i,] for 0's.

17 ## M = number of metabolites,

18 ## tau[1:M] : precision for observed metabolite rates.

19 ## lod = vector of limits of detection, tau = vector of precisions for observed

20 ## metabolite rates.

21 ## dstart[1:N.nontriv], dstop[1:N.nontriv]: indices into delta and phi

22 ## such that, e.g., delta[dstart[i]:dstop[i]] are the deltas for row i of Phi.

23 ## dstart[1] == 1, dstop[N.nontriv] == Ndelta.

24 ##

25 ## The following variables need initializing

26 ## P[i], i in 1:N

27 ## a

28 ## m.R

29 ## m.U

30 ## X1

31 ## X2

32 ## delta[1:Ndelta]: the non-trivial deltas corresponding to the non-trivial phis

33

34 ## Data are based on John W.'s RData file DataTables-040212.RData, and successors.

35

36 attach(file.path("..","MCMCdata.RData"))

37

38 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

39 ## JAGS inputs: nhanesdata and inits

40 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

41

42 nhanesdata <- list()

43 ## The Jags models need these data:

44 ## IURcuts - cutpoints for the IUR production volumes (cutpoints above)

45 nhanesdata$IURcuts <- log(cutpoints)

46 ## IURindx - vector indexing the rows of pred.data with IUR production volumes

47 nhanesdata$IURindx <- grep("IUR HPV",pred.data$PUsource)

48 ## IURvols - indicator for which interval an IUR production volume falls

49 nhanesdata$IURvols <-

50 as.numeric(pred.data$V.cat[nhanesdata$IURindx[1:length(nhanesdata$IURindx)]])

51 ## CPRIindx indexes the rows of pred.data with CPRI use volumes

52 nhanesdata$CPRIindx <- which(pred.data$PUsource == "CPRI 2002")

53 nhanesdata$lV.CPRI <- log(pred.data[nhanesdata$CPRIindx,"V"])

54 ## ls.R, ls.U are the log unit predictions from RADAR and USEtox, respectively

55 nhanesdata$ls.R <- log(pred.data$s.R)

56 nhanesdata$ls.U <- log(pred.data$s.U)

57 nhanesdata$NearField <- pred.data$NearField

58 ## Mean of ls.R and ls.U, for centering
59 nhanesdata$meanls.R <- mean(nhanesdata$ls.R)

60 nhanesdata$meanls.U <- mean(nhanesdata$ls.U)

61

62 ## ly - the log of the observed metabolite quantiles

63 nhanesdata$ly <-

64 ly <-

65 t(log(data.matrix(Measured[,c("Value50","Value75","Value90","Value95")])))

66 ## Q, the proportions associated with the quantiles

67 nhanesdata$Q <- Q <- c(0.5,0.75,0.9,0.95)

68 ## se - precision of measured metabolite values, from reported confidence limits on 

69 quantiles

70 nhanesdata$se <- Measured$se

71 ## llod - log(Limit of Detection)

72 nhanesdata$llod <- log(Measured$LOQ)

73 ## is.above.LOD - vector with 0 for below lod, 1 for above it

74 nhanesdata$is.above.LOD <- 0 + !is.na(ly)

75 ## Some counts:

76 ## N is the number of parents

77 nhanesdata$N <- nrow(pred.data)

78 ## M is the number of metabolites

79 nhanesdata$M <- nrow(Measured)

80 ## N.IUR is the number of IUR reported production volumes (length of IURindx)

81 nhanesdata$N.IUR <- length(nhanesdata$IURindx)

82 ## N.CPRI is the number of CPRI reported use volumes (length of CPRIindx)

83 nhanesdata$N.CPRI <- length(nhanesdata$CPRIindx)

84 ## N.nontriv is the number of parents that result in >1 metabolite

85 nhanesdata$N.nontriv <- sum(N.metabs > 1)

86 ## Nmets - number of metabolites for parent[i], when > 1

87 nhanesdata$Nmets <- N.metabs[N.metabs > 1]

88 ## index - index[i,j] is the index into the metabolite list for the jth metabolite

89 ## of the ith parent

90 nhanesdata$indx <-

91 t(apply(Xref,1,function(x)c((1:length(x))[is.na(x)],

92 (1:length(x))[!is.na(x)])))

93 ## begin and end index into phi[] for each row of Phi

94 tmp <- cumsum(nhanesdata$Nmets)

95 nhanesdata$dstart <- 1+c(0,tmp[-nhanesdata$N.nontriv])

96 nhanesdata$dstop <- tmp

97

98 ## Phi

99 nhanesdata$Phi <- Xref

100 nhanesdata$Ndelta <- sum(is.na(Xref))

101

102 ## ---------------------------------------------------------------------

103 ## Initialization

104

105 inits <- list()

106 inits$lV.IUR <- log(pred.data$V[nhanesdata$IURindx])

107 inits$sdlVmu.dum <- sd(log(pred.data$V))

108 inits$lVmu <- mean(log(pred.data$V))

109 inits$lV.CPRI.mu <- log(pred.data$V[nhanesdata$CPRIindx])

110 ## inits$lP: get estimates for lU, then work backwards to lP

111 ## assuming equal distribution of parent to metabolites when unknown

112 ## JAGS code

113

114 model <- jags.model("model0.jag",data=nhanesdata,inits=inits,n.chains=1,n.adapt=45000)

115

116

117 ## Now, run for 50000 samples, thinned every 100, and dump lP, lU, and delta

118 ## as well. We will use these to get starting values for three new chains.

119

120 out.samples <- jags.samples(model,

121 variable.names=c("lV","sd.lVmu","lVmu"),
122 n.iter=100000, thin=200)

123 out.coda <- lapply(out.samples, as.mcmc.list)

124 sapply(out.coda, function(x) try(heidel.diag(x)))

125 lapply(out.coda, function(z) summary(z))

126

127 save(out.samples, out.coda, nhanesdata, inits, file="m0parms1.RData")

128 ## --------------------------------------------------------------------------------

129 ## This is the initial, hopefully converged, chain. Now, use these to start with

130 ## and run three new chains.

131

132

133 ## getInits finds extreme values 

134 getInits <- function(X, ninits=3) {

135 ## X is a list of arrays; convert it into a matrix

136 tmp1 <- lapply(X, function(z) z[,,1])

137 ## combine the elements, transpose it (so variables are columns, replicates are 

138 rows),

139 ## and center and scale the result.

140 dta <- scale(t(do.call("rbind",tmp1)))

141 ## Ultimately, we are looking for ninits indices into the elements of X

142 outinits <- numeric(ninits)

143 ## Pick ninits random directions. Project the data in dta on each direction

144 ## and pick the 5th percentile (since the original direction u was random u and -u 

145 are equally likely).

146 for (i in 1:ninits) {

147 u <- runif(ncol(dta), -1,1)

148 u <- u/sqrt(sum(u*u))

149 lngt <- dta %*% u

150 indx <- seq(along=lngt)

151 keep <- round(0.05*length(indx)) ## this works fine as long as length(indx) > 20

152 oind <- order(lngt)

153 outinits[i] <- oind[keep]

154 }

155 outinits

156 }

157 ## Get 3 hyperdispersed points.

158 Indx <- getInits(out.samples,3)

159 

160 inits2 <- list(list(),list(),list())

161 for (i in 1:3) {

162 inits2[[i]]$lV.IUR <- out.samples[["lV"]][nhanesdata$IURindx,Indx[i],1]

163 inits2[[i]]$sdlVmu.dum <- out.samples[["sd.lVmu"]][1,Indx[i],1]

164 inits2[[i]]$lVmu <- out.samples[["lVmu"]][1,Indx[i],1]

165 inits2[[i]]$lV.CPRI.mu <- out.samples[["lV"]][nhanesdata$CPRIindx,Indx[i],1]

166 }

167

168

169 model <- jags.model("model0.jag",data=nhanesdata,inits=inits2,n.chains=3,

170 n.adapt=45000)

171

172 adaptdone <- FALSE

173 Count <- 0

174 while (!adaptdone && (Count < 10)) {

175 adaptdone <- adapt(model, 50000, end.adaptation=FALSE)

176 Count <- Count + 1

177 }

178

179 out.samps3R <- jags.samples(model,

180 variable.names=c("lV","sd.lVmu","lVmu"),

181 n.iter=500000, thin=500)

182 out.coda3R <- lapply(out.samps3R, as.mcmc.list)

183 save(out.samps3R, out.coda3R, file="m0parms3.RData")

184
185 for (z in names(out.coda3R)) {

186 writeLines(paste("Variable: ",z))

187 print(try(heidel.diag(out.coda3R[[z]])))

188 print(try(gelman.diag(out.coda3R[[z]])))

189 print(summary(out.coda3R[[z]]))

190 writeLines(c(" ","--------------------------------"," "))

191 }

192

193

194

195 
OnlyV/model0.jag: JAGS Code to Quantify Uncertainty in Production 

Volume

1 model

2 {

3 ## Production volumes are crudely estimated or reported. We have two

4 ## sources, with their own error characteristics. The CPRI data are

5 ## estimated from reported use data, as point estimates. The EPA IUR data

6 ## are reported in predetermined bins. We treat the overall distribution

7 ## of use/production as lognormal (since everything here is on a log

8 ## scale, that means lots of normal distributions). IUR data are then

9 ## interval-censored, with predefined bins reported in IURcuts,

10 ## while each CPRI estimate is sampled from a lognormal distribution

11 ## with unknown median and coefficient of variation depending on the

12 ## estimated error of the estimate. The means are distributed according

13 ## to the same distribution as the IUR data.

14 ## Model parental dose rates as a function of RAIDAR and USEtox

15 ## predictions

16 ## IUR volumes

17 for (i in 1:N.IUR) {

18 IURvols[i] ~ dinterval(lV.IUR[i], IURcuts)

19 lV.IUR[i] ~ dnorm(lVmu, tau.lVmu)

20 lV[IURindx[i]] <- lV.IUR[i]

21 }

22 lVmu ~ dnorm(0, 0.001)

23 sdlVmu.dum ~ dt(0, 1/625, 1)

24 sd.lVmu <- abs(sdlVmu.dum)

25 tau.lVmu <- 1/pow(sd.lVmu,2)

26 ## -------------------------------------------

27 ## CPRI usage

28 for (i in 1:N.CPRI) {

29 lV.CPRI[i] ~ dnorm(lV.CPRI.mu[i], 25) ## CV of 20%: 25 <- 1/(.2^2)

30 lV.CPRI.mu[i] ~ dnorm(lVmu, tau.lVmu) 

31 lV[CPRIindx[i]] <- lV.CPRI.mu[i]

32 }

33 }
analysisOmnibus-Fits.R: Evaluating Alternative Models for Predicting 

NHANES Chemical Exposures

1 library(compiler)

2 setCompilerOptions(optimize=3)

3 enableJIT(3)

4 library(nlme)

5 library(mvtnorm)

6 library(multicore)

7 source("QFunctions.R")

8 ## This takes a two-stage approach to the problem of validation and calibration

9 ## of RAIDAR and USEtox against NHANES. In the first stage, samples from the

10 ## posterior distribution of parent exposures are drawn using MCMC methods,

11 ## using code in the 'OnlyP' directory, AND samples from the production volume

12 ## posteriors are drawn using code in the 'OnlyV' directory. See those

13 ## models for the details.

14 ## The second stage is to do a conventional analysis, looking for

15 ## relationships between lP and ls.R, ls.U, and lV (and combinations thereof),

16 ## AND to look for influential observations. This analysis is done for each

17 ## sample in lV and lP computed in the first stage. We use the distributions

18 ## of the various statistics over the samples of lP and lV to make our

19 ## inferences.

20

21 ## FittingModels controls whether we actually fit models in

22 ## this run, or just print output (mainly for fixing

23 ## bugs in the output code)

24

25 FittingModels <- TRUE

26

27 ## Attach the initial data:

28

29 attach(file.path("..","MCMCdata.RData"))

30

31 ## Get the lP samples

32

33 attach(file.path("..","OnlyP","m3parms3.RData"))

34

35 tmp <- do.call("rbind",out.coda3R)

36 indx <- grep("^lP\\[",colnames(tmp))

37 lP.mx <- tmp[,indx]

38

39 detach(2)

40

41 attach(file.path("..","OnlyV", "m0parms3.RData"))

42

43 lV.mx <- do.call("rbind",out.coda3R$lV)

44 ## We use half the sample to estimate correlation and variances,

45 ## the other half to do the testing.

46 ## Get a random sample of half the rows for testing:

47

48 if (FittingModels) {

49

50 CorSampIndx <- sample(3000,1500)

51 RegSampIndx <- (1:3000)[-CorSampIndx]

52

53 ## Now compute the covariance matrix for lV

54 lP.cov <- cov(lP.mx[CorSampIndx,])

55 lP.cor <- cov2cor(lP.cov)

56

57 corparms <- lP.cor[lower.tri(lP.cor)]

58 lP.vars <- diag(lP.cov)
59

60

61 ## This will be easier if everything is in a dataframe. Everything

62 ## but lP and lV will be a constant, so just build it, and leave space

63 ## for those two variables, which will be assigned at each iteration.

64 ## Use the first row of each variable to initialize it.

65

66

67 ## Build a function that carries out the analysis we want to do for

68 ## each pair of lP and lV in lP.mx and lV.mx. The analysis is

69 ## designed to estimate parameters for each of several models,

70 ## corresponding to models fit in the full Bayesian analysis, and

71 ## compute statistics designed to help evaluate the importance of

72 ## different factors, as well as the best predictive model for lP. We

73 ## do this for two reasons: the full Bayesian model allows lP and lV

74 ## to move around, because lV is censored or measured with error, and

75 ## lP is really derived from the measurements of metabolites

76 ## compounded with unknown stoichiometry, which gives a slightly

77 ## mysterious feel to evaluating fits; computing, e.g., DIC for this

78 ## model seems to be complicated, and is certainly not immediately

79 ## obvious.

80

81 ## It seems natural to use gls() to do the fitting, but I cannot

82 ## figure out how to get it to handle the covariance structure I need

83 ## for this problem. The model is lP = f(predictions, parameters) +

84 ## epsilon + zeta epsilon is normal(0, sigma^2), and zeta is MVnorm(0,

85 ## defined covariance), where the defined covariance is estimated from

86 ## the MCMC samples. To carry this out, I make a new fitting function

87 ## mvtnormfit(), which takes as arguments:

88 ## y

89 ## x1,...,x3 (not all required)

90 ## M (model number, according to the above scheme)

91 ## lPcov: the defined covariance matrix

92 ## Return value is:

93 ## loglik: log likelihood

94 ## AIC: AIC

95 ## parameters: named appropriately

96 ## R2: defined as 1 - var(y - pred) / var(y)

97

98

99 NMs <- c("intercept_0",

100 "sigma_0",

101 "AIC_0","LL_0",

102 

103 "intercept_1","Field_1",

104 "sigma_1",

105 "AIC_1","R2_1","LL_1",

106 

107 "intercept_2","Field_2","m.lV_2",

108 "sigma_2",

109 "AIC_2","R2_2","LL_2",

110 

111 "intercept_3","Field_3","m.lR_3",

112 "sigma_3",

113 "AIC_3","R2_3","LL_3",

114 

115 "intercept_4","Field_4","m.luR_4","m.lV_4",

116 "sigma_4",

117 "AIC_4","R2_4","LL_4",

118 

119 "intercept_5","Field_5","m.lU_5",

120 "sigma_5",

121 "AIC_5","R2_5","LL_5",
122 

123 "intercept_6","Field_6","m.luU_6","m.lV_6",

124 "sigma_6",

125 "AIC_6","R2_6","LL_6",

126 

127 "intercept_7","Field_7","m.luR_7","m.luU_7","m.lV_7",

128 "sigma_7",

129 "AIC_7","R2_7","LL_7",

130 

131 "intercept_8","Field_8","m.lR_8","m.lU_8",

132 "sigma_8",

133 "AIC_8","R2_8","LL_8",

134 

135 "intercept_9","Field_9","m.luR_9","m.luU_9","m.lV_9",

136 "m.luR.N_9",

137 "sigma_9",

138 "AIC_9", "R2_9","LL_9",

139 

140 "intercept_10","Field_10","m.luR_10","m.luU_10","m.lV_10",

141 "m.luU.N_10",

142 "sigma_10",

143 "AIC_10", "R2_10","LL_10",

144 

145 "intercept_11","Field_11","m.luR_11","m.luU_11","m.lV_11",

146 "m.lV.N_11",

147 "sigma_11",

148 "AIC_11", "R2_11","LL_11",

149

150 "intercept_12","Field_12","m.luR_12","m.luU_12","m.lV_12",

151 "m.luR.N_12","m.luU.N_12","m.lV.N_12",

152 "sigma_12",

153 "AIC_12", "R2_12","LL_12",

154 

155 "intercept_13","Field_13","m.luR_13","m.luU_13","m.lV_13",

156 "m.luR.N_13","m.luU.N_13",

157 "sigma_13",

158 "AIC_13", "R2_13","LL_13",

159 

160 "intercept_14","Field_14","m.luR_14","m.luU_14","m.lV_14",

161 "m.luU.N_14","m.lV.N_14",

162 "sigma_14",

163 "AIC_14", "R2_14","LL_14",

164 

165 "intercept_15","Field_15","m.luR_15","m.luU_15","m.lV_15",

166 "m.luR.N_15","m.lV.N_15",

167 "sigma_15",

168 "AIC_15", "R2_15","LL_15",

169 

170 "intercept_16","Field_16",

171 "m.luR.N_16","m.luU.N_16","m.lV.N_16",

172 "sigma_16",

173 "AIC_16", "R2_16","LL_16",

174 

175 "intercept_17","Field_17",

176 "m.luR.N_17","m.luU.N_17",

177 "sigma_17",

178 "AIC_17", "R2_17","LL_17",

179 

180 "intercept_18","Field_18",

181 "m.luR.N_18","m.lV.N_18",

182 "sigma_18",

183 "AIC_18", "R2_18","LL_18",

184 
185 "intercept_19","Field_19",

186 "m.luU.N_19","m.lV.N_19",

187 "sigma_19",

188 "AIC_19", "R2_19","LL_19",

189 

190 "intercept_20","Field_20",

191 "m.luR.N_20",

192 "sigma_20",

193 "AIC_20", "R2_20","LL_20",

194 

195 "intercept_21", "m.lR_21",

196 "sigma_21",

197 "AIC_21","R2_21","LL_21",

198 

199 "intercept_22", "m.lU_22",

200 "sigma_22",

201 "AIC_22","R2_22","LL_22",

202 

203 "dR0<->1","dR0<->21","dR0<->22",

204 "dR1<->2","dR1<->3","dR1<->4","dR1<->5","dR1<->6","dR1<->12",

205 "dR2<->4","dR2<->6",

206 "dR3<->4",

207 "dR4<->7",

208 "dR5<->6",

209 "dR6<->7",

210 "dR7<->9","dR7<->10","dR7<->11","dR7<->12",

211 "dR8<->7",

212 "dR9<->12",

213 "dR10<->12",

214 "dR11<->12",

215 "dR16<->12",

216 "dR18<->12",

217 "dR20<->12",

218 

219 "dLL0<->1","dLL0<->21","dLL0<->22",

220 "dLL1<->2","dLL1<->3","dLL1<->4","dLL1<->5","dLL1<->6","dLL1<->12",

221 "dLL2<->4","dLL2<->6",

222 "dLL3<->4",

223 "dLL4<->7",

224 "dLL5<->6",

225 "dLL6<->7",

226 "dLL7<->9","dLL7<->10","dLL7<->11","dLL7<->12",

227 "dLL8<->7",

228 "dLL9<->12",

229 "dLL10<->12",

230 "dLL11<->12",

231 "dLL16<->12",

232 "dLL18<->12",

233 "dLL20<->12"

234 

235 )

236

237 analysis <- function(df,drop=NULL) {

238 if (!is.null(drop)) {

239 df <- df[-drop,]

240 }

241 out1 <- rep(as.numeric(NA),times=length(NMs))

242 names(out1) <- NMs

243 

244 ## Model 0

245 md0 <- mvtnormfit(df, 0, lP.cov)

246 if (!is(md0, "try-error")) {

247 out1["intercept_0"] <- unname(coef(md0)["intercept"])
248 out1["sigma_0"] <- exp(unname(coef(md0)["lsigma"]))

249 out1["AIC_0"] <- md0$AIC

250 out1["LL_0"] <- md0$LogLik

251 }

252 

253 ## Model 1: 

254 md1 <- mvtnormfit(df,1,lP.cov)

255 if (!is(md1, "try-error")) {

256 out1["intercept_1"] <- unname(coef(md1)["intercept"])

257 out1["Field_1"] <- unname(coef(md1)["NearField"])

258 out1["sigma_1"] <- exp(unname(coef(md1)["lsigma"]))

259 out1["AIC_1"] <- md1$AIC

260 out1["LL_1"] <- md1$LogLik

261 out1["R2_1"] <- md1$R2

262 }

263 

264 ## Model 2

265 md2 <- mvtnormfit(df,2,lP.cov)

266 if (!is(md2, "try-error")) {

267 out1["intercept_2"] <- unname(coef(md2)["intercept"])

268 out1["Field_2"] <- unname(coef(md2)["NearField"])

269 out1["m.lV_2"] <- unname(coef(md2)["m.lV"])

270 out1["sigma_2"] <- exp(unname(coef(md2)["lsigma"]))

271 out1["AIC_2"] <- md2$AIC

272 out1["LL_2"] <- md2$LogLik

273 out1["R2_2"] <- md2$R2

274 }

275 

276 ## Model 3

277 md3 <- mvtnormfit(df, 3, lP.cov)

278 if (!is(md3, "try-error")){

279 out1["intercept_3"] <- unname(coef(md3)["intercept"])

280 out1["Field_3"] <- unname(coef(md3)["NearField"])

281 out1["m.lR_3"] <- unname(coef(md3)["m.lR"])

282 out1["sigma_3"] <- exp(unname(coef(md3)["lsigma"]))

283 out1["AIC_3"] <- md3$AIC

284 out1["LL_3"] <- md3$LogLik

285 out1["R2_3"] <- md3$R2

286 }

287 

288 ## Model 4

289 md4 <- mvtnormfit(df, 4, lP.cov)

290 if (!is(md4, "try-error")){

291 out1["intercept_4"] <- unname(coef(md4)["intercept"])

292 out1["Field_4"] <- unname(coef(md4)["NearField"])

293 out1["m.luR_4"] <- unname(coef(md4)["m.luR"])

294 out1["m.lV_4"] <- unname(coef(md4)["m.lV"])

295 out1["sigma_4"] <- exp(unname(coef(md4)["lsigma"]))

296 out1["AIC_4"] <- md4$AIC

297 out1["LL_4"] <- md4$LogLik

298 out1["R2_4"] <- md4$R2

299 }

300 ## out2 <- influence.peddling(1, df, md4, lP.cor)

301 

302 ## Model 5

303 md5 <- mvtnormfit(df, 5, lP.cov)

304 if (!is(md5, "try-error")) {

305 out1["intercept_5"] <- unname(coef(md5)["intercept"])

306 out1["Field_5"] <- unname(coef(md5)["NearField"])

307 out1["m.lU_5"] <- unname(coef(md5)["m.lU"])

308 out1["sigma_5"] <- exp(unname(coef(md5)["lsigma"]))

309 out1["AIC_5"] <- md5$AIC

310 out1["LL_5"] <- md5$LogLik
311 out1["R2_5"] <- md5$R2

312 }

313 

314 ## Model 6

315 md6 <- mvtnormfit(df, 6, lP.cov)

316 if (!is(md6, "try-error")){

317 out1["intercept_6"] <- unname(coef(md6)["intercept"])

318 out1["Field_6"] <- unname(coef(md6)["NearField"])

319 out1["m.luU_6"] <- unname(coef(md6)["m.luU"])

320 out1["m.lV_6"] <- unname(coef(md6)["m.lV"])

321 out1["sigma_6"] <- exp(unname(coef(md6)["lsigma"]))

322 out1["AIC_6"] <- md6$AIC

323 out1["LL_6"] <- md6$LogLik

324 out1["R2_6"] <- md6$R2

325 }

326 ## out3 <- influence.peddling(2, df, md6, lP.cor)

327 

328 ## Model 7

329 md7 <- mvtnormfit(df, 7, lP.cov)

330 if (!is(md7, "try-error")) {

331 out1["intercept_7"] <- unname(coef(md7)["intercept"])

332 out1["Field_7"] <- unname(coef(md7)["NearField"])

333 out1["m.luR_7"] <- unname(coef(md7)["m.luR"])

334 out1["m.luU_7"] <- unname(coef(md7)["m.luU"])

335 out1["m.lV_7"] <- unname(coef(md7)["m.lV"])

336 out1["sigma_7"] <- exp(unname(coef(md7)["lsigma"]))

337 out1["AIC_7"] <- md7$AIC

338 out1["LL_7"] <- md7$LogLik

339 out1["R2_7"] <- md7$R2

340 }

341

342 ## Model 8

343 md8 <- mvtnormfit(df, 8, lP.cov)

344 if (!is(md8, "try-error")) {

345 out1["intercept_8"] <- unname(coef(md8)["intercept"])

346 out1["Field_8"] <- unname(coef(md8)["NearField"])

347 out1["m.lR_8"] <- unname(coef(md8)["m.lR"])

348 out1["m.lU_8"] <- unname(coef(md8)["m.lU"])

349 out1["sigma_8"] <- exp(unname(coef(md8)["lsigma"]))

350 out1["AIC_8"] <- md8$AIC

351 out1["LL_8"] <- md8$LogLik

352 out1["R2_8"] <- md8$R2

353 }

354 

355 ## Model 9

356 md9 <- mvtnormfit(df, 9, lP.cov)

357 if (!is(md9, "try-error")){

358 out1["intercept_9"] <- unname(coef(md9)["intercept"])

359 out1["Field_9"] <- unname(coef(md9)["NearField"])

360 out1["m.luR_9"] <- unname(coef(md9)["m.luR"])

361 out1["m.luU_9"] <- unname(coef(md9)["m.luU"])

362 out1["m.lV_9"] <- unname(coef(md9)["m.lV"])

363 out1["m.luR.N_9"] <- unname(coef(md9)["m.luR.N"])

364 out1["sigma_9"] <- exp(unname(coef(md9)["lsigma"]))

365 out1["AIC_9"] <- md9$AIC

366 out1["LL_9"] <- md9$LogLik

367 out1["R2_9"] <- md9$R2

368 }

369

370 ## Model 10

371 md10 <- mvtnormfit(df, 10, lP.cov)

372 if (!is(md10, "try-error")){

373 out1["intercept_10"] <- unname(coef(md10)["intercept"])
374 out1["Field_10"] <- unname(coef(md10)["NearField"])

375 out1["m.luR_10"] <- unname(coef(md10)["m.luR"])

376 out1["m.luU_10"] <- unname(coef(md10)["m.luU"])

377 out1["m.lV_10"] <- unname(coef(md10)["m.lV"])

378 out1["m.luU.N_10"] <- unname(coef(md10)["m.luU.N"])

379 out1["sigma_10"] <- exp(unname(coef(md10)["lsigma"]))

380 out1["AIC_10"] <- md10$AIC

381 out1["LL_10"] <- md10$LogLik

382 out1["R2_10"] <- md10$R2

383 }

384

385 ## Model 11

386 md11 <- mvtnormfit(df, 11, lP.cov)

387 if (!is(md11, "try-error")){

388 out1["intercept_11"] <- unname(coef(md11)["intercept"])

389 out1["Field_11"] <- unname(coef(md11)["NearField"])

390 out1["m.luR_11"] <- unname(coef(md11)["m.luR"])

391 out1["m.luU_11"] <- unname(coef(md11)["m.luU"])

392 out1["m.lV_11"] <- unname(coef(md11)["m.lV"])

393 out1["m.lV.N_11"] <- unname(coef(md11)["m.lV.N"])

394 out1["sigma_11"] <- exp(unname(coef(md11)["lsigma"]))

395 out1["AIC_11"] <- md11$AIC

396 out1["LL_11"] <- md11$LogLik

397 out1["R2_11"] <- md11$R2

398 }

399

400 ## Model 12

401 md12 <- mvtnormfit(df, 12, lP.cov)

402 if (!is(md12, "try-error")){

403 out1["intercept_12"] <- unname(coef(md12)["intercept"])

404 out1["Field_12"] <- unname(coef(md12)["NearField"])

405 out1["m.luR_12"] <- unname(coef(md12)["m.luR"])

406 out1["m.luU_12"] <- unname(coef(md12)["m.luU"])

407 out1["m.lV_12"] <- unname(coef(md12)["m.lV"])

408 out1["m.luR.N_12"] <- unname(coef(md12)["m.luR.N"])

409 out1["m.luU.N_12"] <- unname(coef(md12)["m.luU.N"])

410 out1["m.lV.N_12"] <- unname(coef(md12)["m.lV.N"])

411 out1["sigma_12"] <- exp(unname(coef(md12)["lsigma"]))

412 out1["AIC_12"] <- md12$AIC

413 out1["LL_12"] <- md12$LogLik

414 out1["R2_12"] <- md12$R2

415 }

416 ## Model 13

417 md13 <- mvtnormfit(df, 13, lP.cov)

418 if (!is(md13, "try-error")){

419 out1["intercept_13"] <- unname(coef(md13)["intercept"])

420 out1["Field_13"] <- unname(coef(md13)["NearField"])

421 out1["m.luR_13"] <- unname(coef(md13)["m.luR"])

422 out1["m.luU_13"] <- unname(coef(md13)["m.luU"])

423 out1["m.lV_13"] <- unname(coef(md13)["m.lV"])

424 out1["m.luR.N_13"] <- unname(coef(md13)["m.luR.N"])

425 out1["m.luU.N_13"] <- unname(coef(md13)["m.luU.N"])

426 out1["sigma_13"] <- exp(unname(coef(md13)["lsigma"]))

427 out1["AIC_13"] <- md13$AIC

428 out1["LL_13"] <- md13$LogLik

429 out1["R2_13"] <- md13$R2

430 }

431 ## Model 14

432 md14 <- mvtnormfit(df, 14, lP.cov)

433 if (!is(md14, "try-error")){

434 out1["intercept_14"] <- unname(coef(md14)["intercept"])

435 out1["Field_14"] <- unname(coef(md14)["NearField"])

436 out1["m.luR_14"] <- unname(coef(md14)["m.luR"])
437 out1["m.luU_14"] <- unname(coef(md14)["m.luU"])

438 out1["m.lV_14"] <- unname(coef(md14)["m.lV"])

439 out1["m.luU.N_14"] <- unname(coef(md14)["m.luU.N"])

440 out1["m.lV.N_14"] <- unname(coef(md14)["m.lV.N"])

441 out1["sigma_14"] <- exp(unname(coef(md14)["lsigma"]))

442 out1["AIC_14"] <- md14$AIC

443 out1["LL_14"] <- md14$LogLik

444 out1["R2_14"] <- md14$R2

445 }

446 

447 ## Model 15

448 md15 <- mvtnormfit(df, 15, lP.cov)

449 if (!is(md15, "try-error")){

450 out1["intercept_15"] <- unname(coef(md15)["intercept"])

451 out1["Field_15"] <- unname(coef(md15)["NearField"])

452 out1["m.luR_15"] <- unname(coef(md15)["m.luR"])

453 out1["m.luU_15"] <- unname(coef(md15)["m.luU"])

454 out1["m.lV_15"] <- unname(coef(md15)["m.lV"])

455 out1["m.luR.N_15"] <- unname(coef(md15)["m.luR.N"])

456 out1["m.lV.N_15"] <- unname(coef(md15)["m.lV.N"])

457 out1["sigma_15"] <- exp(unname(coef(md15)["lsigma"]))

458 out1["AIC_15"] <- md15$AIC

459 out1["LL_15"] <- md15$LogLik

460 out1["R2_15"] <- md15$R2

461 }

462 

463 ## Model 16

464 md16 <- mvtnormfit(df, 16, lP.cov)

465 if (!is(md16, "try-error")){

466 out1["intercept_16"] <- unname(coef(md16)["intercept"])

467 out1["Field_16"] <- unname(coef(md16)["NearField"])

468 out1["m.luR.N_16"] <- unname(coef(md16)["m.luR.N"])

469 out1["m.luU.N_16"] <- unname(coef(md16)["m.luU.N"])

470 out1["m.lV.N_16"] <- unname(coef(md16)["m.lV.N"])

471 out1["sigma_16"] <- exp(unname(coef(md16)["lsigma"]))

472 out1["AIC_16"] <- md16$AIC

473 out1["LL_16"] <- md16$LogLik

474 out1["R2_16"] <- md16$R2

475 }

476 

477 ## Model 17

478 md17 <- mvtnormfit(df, 17, lP.cov)

479 if (!is(md17, "try-error")){

480 out1["intercept_17"] <- unname(coef(md17)["intercept"])

481 out1["Field_17"] <- unname(coef(md17)["NearField"])

482 out1["m.luR.N_17"] <- unname(coef(md17)["m.luR.N"])

483 out1["m.luU.N_17"] <- unname(coef(md17)["m.luU.N"])

484 out1["sigma_17"] <- exp(unname(coef(md17)["lsigma"]))

485 out1["AIC_17"] <- md17$AIC

486 out1["LL_17"] <- md17$LogLik

487 out1["R2_17"] <- md17$R2

488 }

489 ## Model 18

490 md18 <- mvtnormfit(df, 18, lP.cov)

491 if (!is(md18, "try-error")){

492 out1["intercept_18"] <- unname(coef(md18)["intercept"])

493 out1["Field_18"] <- unname(coef(md18)["NearField"])

494 out1["m.luR.N_18"] <- unname(coef(md18)["m.luR.N"])

495 out1["m.lV.N_18"] <- unname(coef(md18)["m.lV.N"])

496 out1["sigma_18"] <- exp(unname(coef(md18)["lsigma"]))

497 out1["AIC_18"] <- md18$AIC

498 out1["LL_18"] <- md18$LogLik

499 out1["R2_18"] <- md18$R2
500 }

501

502 ## Model 19

503 md19 <- mvtnormfit(df, 19, lP.cov)

504 if (!is(md19, "try-error")){

505 out1["intercept_19"] <- unname(coef(md19)["intercept"])

506 out1["Field_19"] <- unname(coef(md19)["NearField"])

507 out1["m.luU.N_19"] <- unname(coef(md19)["m.luU.N"])

508 out1["m.lV.N_19"] <- unname(coef(md19)["m.lV.N"])

509 out1["sigma_19"] <- exp(unname(coef(md19)["lsigma"]))

510 out1["AIC_19"] <- md19$AIC

511 out1["LL_19"] <- md19$LogLik

512 out1["R2_19"] <- md19$R2

513 }

514 

515 ## Model 20

516 md20 <- mvtnormfit(df, 20, lP.cov)

517 if (!is(md20, "try-error")){

518 out1["intercept_20"] <- unname(coef(md20)["intercept"])

519 out1["Field_20"] <- unname(coef(md20)["NearField"])

520 out1["m.luR.N_20"] <- unname(coef(md20)["m.luR.N"])

521 out1["sigma_20"] <- exp(unname(coef(md20)["lsigma"]))

522 out1["AIC_20"] <- md20$AIC

523 out1["LL_20"] <- md20$LogLik

524 out1["R2_20"] <- md20$R2

525 }

526

527 ## Model 21

528 md21 <- mvtnormfit(df, 21, lP.cov)

529 if (!is(md21, "try-error")){

530 out1["intercept_21"] <- unname(coef(md21)["intercept"])

531 out1["m.lR_21"] <- unname(coef(md21)["m.lR"])

532 out1["sigma_21"] <- exp(unname(coef(md21)["lsigma"]))

533 out1["AIC_21"] <- md21$AIC

534 out1["LL_21"] <- md21$LogLik

535 out1["R2_21"] <- md21$R2

536 }

537

538 ## Model 22

539 md22 <- mvtnormfit(df, 22, lP.cov)

540 if (!is(md22, "try-error")) {

541 out1["intercept_22"] <- unname(coef(md22)["intercept"])

542 out1["m.lU_22"] <- unname(coef(md22)["m.lU"])

543 out1["sigma_22"] <- exp(unname(coef(md22)["lsigma"]))

544 out1["AIC_22"] <- md22$AIC

545 out1["LL_22"] <- md22$LogLik

546 out1["R2_22"] <- md22$R2

547 }

548 

549 ## Model Comparisons: record the null-hypothetical delta-R2 and

550 ## delta log-likelihood for the corresponding two models. Construct

551 ## a new null-hypothetical dataset for each test.

552

553 ## Try each set of fits 3 times before giving up

554 df0 <- df

555 ## newDataset from md0

556 if (!is(md0, "try-error")) {

557 muhat <- rep(coef(md0)["intercept"], nrow(df))

558 Sigma <- lP.cov +

559 diag(exp(2*coef(md0)["lsigma"]),

560 nrow(lP.cov), ncol(lP.cov))

561 Nrep <- 0

562 dorep <- TRUE
563 while(dorep && Nrep < 3){

564 Nrep <- Nrep + 1

565 lP <- rmvnorm(1,muhat,sigma=Sigma)

566 df0$lP <- lP[1,]

567 md0a <- mvtnormfit(df0, 0, lP.cov)

568 md1a <- mvtnormfit(df0, 1, lP.cov)

569 if (!is(md0a, "try-error") && !is(md1a,"try-error")) dorep <- FALSE

570 }

571 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

572 out1["dR0<->1"] <- md1a$R2 ## md0a$R2 should be 0

573 out1["dLL0<->1"] <- md1a$LogLik - md0a$LogLik

574 }

575

576 Nrep <- 0

577 dorep <- TRUE

578 while(dorep && Nrep < 3){

579 Nrep <- Nrep + 1

580 lP <- rmvnorm(1,muhat,sigma=Sigma)

581 df0$lP <- lP[1,]

582 md0a <- mvtnormfit(df0, 0, lP.cov)

583 md1a <- mvtnormfit(df0, 21, lP.cov)

584 if (!is(md0a, "try-error") && !is(md1a,"try-error")) dorep <- FALSE

585 }

586 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

587 out1["dR0<->21"] <- md1a$R2

588 out1["dLL0<->21"] <- md1a$LogLik - md0a$LogLik

589 }

590

591 Nrep <- 0

592 dorep <- TRUE

593 while(dorep && Nrep < 3){

594 Nrep <- Nrep + 1

595 lP <- rmvnorm(1,muhat,sigma=Sigma)

596 df0$lP <- lP[1,]

597 md0a <- mvtnormfit(df0, 0, lP.cov)

598 md1a <- mvtnormfit(df0, 22, lP.cov)

599 if (!is(md0a, "try-error") && !is(md1a,"try-error")) dorep <- FALSE

600 }

601 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

602 out1["dR0<->22"] <- md1a$R2

603 out1["dLL0<->22"] <- md1a$LogLik - md0a$LogLik

604 }

605 }

606 ## New datasets from md1

607 if (!is(md1, "try-error")) {

608 muhat <- coef(md1)["intercept"] + coef(md1)["NearField"]*df$NearField

609 Sigma <- lP.cov + diag(exp(2*coef(md1)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

610 Nrep <- 0

611 dorep <- TRUE

612 while(dorep && Nrep < 3){

613 Nrep <- Nrep + 1

614 lP <- rmvnorm(1, muhat, sigma=Sigma)

615 df0$lP <- lP[1,]

616 md0a <- mvtnormfit(df0, 1, lP.cov)

617 md1a <- mvtnormfit(df0, 2, lP.cov)

618 }

619 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

620 out1["dR1<->2"] <- md1a$R2 - md0a$R2

621 out1["dLL1<->2"] <- md1a$LogLik - md0a$LogLik

622 }

623 

624 Nrep <- 0

625 dorep <- TRUE
626 while(dorep && Nrep < 3){

627 Nrep <- Nrep + 1

628 lP <- rmvnorm(1, muhat, sigma=Sigma)

629 df0$lP <- lP[1,]

630 md0a <- mvtnormfit(df0, 1, lP.cov)

631 md1a <- mvtnormfit(df0, 3, lP.cov)

632 }

633 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

634 out1["dR1<->3"] <- md1a$R2 - md0a$R2

635 out1["dLL1<->3"] <- md1a$LogLik - md0a$LogLik

636 }

637 

638 Nrep <- 0

639 dorep <- TRUE

640 while(dorep && Nrep < 3){

641 Nrep <- Nrep + 1

642 lP <- rmvnorm(1, muhat, sigma=Sigma)

643 df0$lP <- lP[1,]

644 md0a <- mvtnormfit(df0, 1, lP.cov)

645 md1a <- mvtnormfit(df0, 4, lP.cov)

646 }

647 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

648 out1["dR1<->4"] <- md1a$R2 - md0a$R2

649 out1["dLL1<->4"] <- md1a$LogLik - md0a$LogLik

650 }

651 

652 Nrep <- 0

653 dorep <- TRUE

654 while(dorep && Nrep < 3){

655 Nrep <- Nrep + 1

656 lP <- rmvnorm(1, muhat, sigma=Sigma)

657 df0$lP <- lP[1,]

658 md0a <- mvtnormfit(df0, 1, lP.cov)

659 md1a <- mvtnormfit(df0, 5, lP.cov)

660 }

661 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

662 out1["dR1<->5"] <- md1a$R2 - md0a$R2

663 out1["dLL1<->5"] <- md1a$LogLik - md0a$LogLik

664 }

665 

666 Nrep <- 0

667 dorep <- TRUE

668 while(dorep && Nrep < 3){

669 Nrep <- Nrep + 1

670 lP <- rmvnorm(1, muhat, sigma=Sigma)

671 df0$lP <- lP[1,]

672 md0a <- mvtnormfit(df0, 1, lP.cov)

673 md1a <- mvtnormfit(df0, 6, lP.cov)

674 }

675 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

676 out1["dR1<->6"] <- md1a$R2 - md0a$R2

677 out1["dLL1<->6"] <- md1a$LogLik - md0a$LogLik

678 }

679 

680 Nrep <- 0

681 dorep <- TRUE

682 while(dorep && Nrep < 3){

683 Nrep <- Nrep + 1

684 lP <- rmvnorm(1, muhat, sigma=Sigma)

685 df0$lP <- lP[1,]

686 md0a <- mvtnormfit(df0, 1, lP.cov)

687 md1a <- mvtnormfit(df0, 12, lP.cov)

688 }
689 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

690 out1["dR1<->12"] <- md1a$R2 - md0a$R2

691 out1["dLL1<->12"] <- md1a$LogLik - md0a$LogLik

692 }

693 }

694 

695 ## New datasets from M2

696 if (!is(md2, "try-error")){

697 muhat <- coef(md2)["intercept"] + coef(md2)["NearField"]*df$NearField + 

698 coef(md2)["m.lV"] * df$lV

699 Sigma <- lP.cov + diag(exp(2*coef(md2)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

700

701 Nrep <- 0

702 dorep <- TRUE

703 while(dorep && Nrep < 3){

704 Nrep <- Nrep + 1

705 lP <- rmvnorm(1, muhat, sigma=Sigma)

706 df0$lP <- lP[1,]

707 md0a <- mvtnormfit(df0, 2, lP.cov)

708 md1a <- mvtnormfit(df0, 4, lP.cov)

709 }

710 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

711 out1["dR2<->4"] <- md1a$R2 - md0a$R2

712 out1["dLL2<->4"] <- md1a$LogLik - md0a$LogLik

713 }

714 

715 Nrep <- 0

716 dorep <- TRUE

717 while(dorep && Nrep < 3){

718 Nrep <- Nrep + 1

719 lP <- rmvnorm(1, muhat, sigma=Sigma)

720 df0$lP <- lP[1,]

721 md0a <- mvtnormfit(df0, 2, lP.cov)

722 md1a <- mvtnormfit(df0, 6, lP.cov)

723 }

724 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

725 out1["dR2<->6"] <- md1a$R2 - md0a$R2

726 out1["dLL2<->6"] <- md1a$LogLik - md0a$LogLik

727 }

728 }

729 

730 ## One new dataset from M3

731 if (!is(md3, "try-error")) {

732 muhat <- coef(md3)["intercept"] + coef(md3)["NearField"]*df$NearField + 

733 coef(md3)["m.lR"] * df$lR

734 Sigma <- lP.cov + diag(exp(2*coef(md3)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

735 Nrep <- 0

736 dorep <- TRUE

737 while(dorep && Nrep < 3){

738 Nrep <- Nrep + 1

739 lP <- rmvnorm(1, muhat, sigma= Sigma)

740 df0$lP <- lP[1,]

741 md0a <- mvtnormfit(df0, 3, lP.cov)

742 md1a <- mvtnormfit(df0, 4, lP.cov)

743 }

744 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

745 out1["dR3<->4"] <- md1a$R2 - md0a$R2

746 out1["dLL3<->4"] <- md1a$LogLik - md0a$LogLik

747 }

748 }

749 

750 ## One new dataset from M4

751 if (!is(md4, "try-error")){
752 muhat <- coef(md4)["intercept"] + coef(md4)["NearField"]*df$NearField +

753 coef(md4)["m.luR"] * df$ls.R + coef(md4)["m.lV"] * df$lV

754 Sigma <- lP.cov + diag(exp(2*coef(md4)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

755 

756 Nrep <- 0

757 dorep <- TRUE

758 while(dorep && Nrep < 3){

759 Nrep <- Nrep + 1

760 lP <- rmvnorm(1, muhat, sigma=Sigma)

761 df0$lP <- lP[1,]

762 md0a <- mvtnormfit(df0, 4, lP.cov)

763 md1a <- mvtnormfit(df0, 7, lP.cov)

764 }

765 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

766 out1["dR4<->7"] <- md1a$R2 - md0a$R2

767 out1["dLL4<->7"] <- md1a$LogLik - md0a$LogLik

768 }

769 }

770 

771 ## One new dataset from M5

772 if (!is(md5, "try-error")) {

773 muhat <- coef(md5)["intercept"] + coef(md5)["NearField"]*df$NearField + 

774 coef(md5)["m.lU"] * df$lU

775 Sigma <- lP.cov + diag(exp(2*coef(md5)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

776 

777 Nrep <- 0

778 dorep <- TRUE

779 while(dorep && Nrep < 3){

780 Nrep <- Nrep + 1

781 lP <- rmvnorm(1, muhat, sigma = Sigma)

782 

783 df0$lP <- lP[1,]

784 md0a <- mvtnormfit(df0, 5, lP.cov)

785 md1a <- mvtnormfit(df0, 6, lP.cov)

786 }

787 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

788 out1["dR5<->6"] <- md1a$R2 - md0a$R2

789 out1["dLL5<->6"] <- md1a$LogLik - md0a$LogLik

790 }

791 }

792 

793 ## One new dataset from M6

794 if (!is(md6, "try-error")) {

795 muhat <- coef(md6)["intercept"] + coef(md6)["NearField"]*df$NearField +

796 coef(md6)["m.luU"] * df$ls.U + coef(md6)["m.lV"] * df$lV

797 Sigma <- lP.cov + diag(exp(2*coef(md6)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

798 

799 Nrep <- 0

800 dorep <- TRUE

801 while(dorep && Nrep < 3){

802 Nrep <- Nrep + 1

803 lP <- rmvnorm(1, muhat, sigma=Sigma)

804

805 df0$lP <- lP[1,]

806 md0a <- mvtnormfit(df0, 6, lP.cov)

807 md1a <- mvtnormfit(df0, 7, lP.cov)

808 }

809 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

810 out1["dR6<->7"] <- md1a$R2 - md0a$R2

811 out1["dLL6<->7"] <- md1a$LogLik - md0a$LogLik

812 }

813 }

814 
815 ## New dataset from M7

816 if (!is(md7, "try-error")) {

817 muhat <- coef(md7)["intercept"] + coef(md7)["NearField"]*df$NearField +

818 coef(md7)["m.luR"] * df$ls.R +

819 coef(md7)["m.luU"] * df$ls.U + coef(md7)["m.lV"] * df$lV

820 Sigma <- lP.cov + diag(exp(2*coef(md7)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

821 

822 Nrep <- 0

823 dorep <- TRUE

824 while(dorep && Nrep < 3){

825 Nrep <- Nrep + 1

826 lP <- rmvnorm(1, muhat, sigma=Sigma)

827 

828 df0$lP <- lP[1,]

829 md0a <- mvtnormfit(df0, 7, lP.cov)

830 md1a <- mvtnormfit(df0, 9, lP.cov)

831 }

832 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

833 out1["dR7<->9"] <- md1a$R2 - md0a$R2

834 out1["dLL7<->9"] <- md1a$LogLik - md0a$LogLik

835 }

836

837 Nrep <- 0

838 dorep <- TRUE

839 while(dorep && Nrep < 3){

840 Nrep <- Nrep + 1

841 lP <- rmvnorm(1, muhat, sigma=Sigma)

842 

843 df0$lP <- lP[1,]

844 md0a <- mvtnormfit(df0, 7, lP.cov)

845 md1a <- mvtnormfit(df0, 10, lP.cov)

846 }

847 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

848 out1["dR7<->10"] <- md1a$R2 - md0a$R2

849 out1["dLL7<->10"] <- md1a$LogLik - md0a$LogLik

850 }

851

852 Nrep <- 0

853 dorep <- TRUE

854 while(dorep && Nrep < 3){

855 Nrep <- Nrep + 1

856 lP <- rmvnorm(1, muhat, sigma=Sigma)

857 

858 df0$lP <- lP[1,]

859 md0a <- mvtnormfit(df0, 7, lP.cov)

860 md1a <- mvtnormfit(df0, 11, lP.cov)

861 }

862 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

863 out1["dR7<->11"] <- md1a$R2 - md0a$R2

864 out1["dLL7<->11"] <- md1a$LogLik - md0a$LogLik

865 }

866 

867 Nrep <- 0

868 dorep <- TRUE

869 while(dorep && Nrep < 3){

870 Nrep <- Nrep + 1

871 lP <- rmvnorm(1, muhat, sigma=Sigma)

872 

873 df0$lP <- lP[1,]

874 md0a <- mvtnormfit(df0, 7, lP.cov)

875 md1a <- mvtnormfit(df0, 12, lP.cov)

876 }

877 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {
878 out1["dR7<->12"] <- md1a$R2 - md0a$R2

879 out1["dLL7<->12"] <- md1a$LogLik - md0a$LogLik

880 }

881 }

882 

883 ## One new dataset from M8

884 if (!is(md8, "try-error")) {

885 muhat <- coef(md8)["intercept"] + coef(md8)["NearField"]*df$NearField +

886 coef(md8)["m.lR"] * df$lR + coef(md8)["m.lU"] * df$lU

887 Sigma <- lP.cov + diag(exp(2*coef(md8)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

888 Nrep <- 0

889 dorep <- TRUE

890 while(dorep && Nrep < 3){

891 Nrep <- Nrep + 1

892 lP <- rmvnorm(1, muhat, sigma=Sigma)

893 

894 df0$lP <- lP[1,]

895 md0a <- mvtnormfit(df0, 8, lP.cov)

896 md1a <- mvtnormfit(df0, 7, lP.cov)

897 }

898 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

899 out1["dR8<->7"] <- md1a$R2 - md0a$R2

900 out1["dLL8<->7"] <- md1a$LogLik - md0a$LogLik

901 }

902 }

903 ## One new dataset from M9

904 if (!is(md9, "try-error")) {

905 muhat <- coef(md9)["intercept"] + coef(md9)["NearField"]*df$NearField +

906 coef(md9)["m.luR"] * df$ls.R + coef(md9)["m.luU"] * df$ls.U +

907 coef(md9)["m.lV"] * df$lV + coef(md9)["m.luR.N"] * df$NearField * df$ls.R

908 Sigma <- lP.cov + diag(exp(2*coef(md9)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

909 

910 Nrep <- 0

911 dorep <- TRUE

912 while(dorep && Nrep < 3){

913 Nrep <- Nrep + 1

914 lP <- rmvnorm(1, muhat, sigma=Sigma)

915 

916 df0$lP <- lP[1,]

917 md0a <- mvtnormfit(df0, 9, lP.cov)

918 md1a <- mvtnormfit(df0, 12, lP.cov)

919 }

920 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

921 out1["dR9<->12"] <- md1a$R2 - md0a$R2

922 out1["dLL9<->12"] <- md1a$LogLik - md0a$LogLik

923 }

924 }

925 

926 ## One new dataset from M10

927 if (!is(md10, "try-error")) {

928 muhat <- coef(md10)["intercept"] + coef(md10)["NearField"]*df$NearField +

929 coef(md10)["m.luR"] * df$ls.R + coef(md10)["m.luU"] * df$ls.U +

930 coef(md10)["m.lV"] * df$lV + coef(md10)["m.luU.N"] * df$NearField * df$ls.U

931 Sigma <- lP.cov + diag(exp(2*coef(md10)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

932 

933 Nrep <- 0

934 dorep <- TRUE

935 while(dorep && Nrep < 3){

936 Nrep <- Nrep + 1

937 lP <- rmvnorm(1, muhat, sigma=lP.cov + diag(exp(2*coef(md10)["lsigma"]), 

938 nrow(lP.cov), ncol(lP.cov)))

939 

940 df0$lP <- lP[1,]
941 md0a <- mvtnormfit(df0, 10, lP.cov)

942 md1a <- mvtnormfit(df0, 12, lP.cov)

943 }

944 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

945 out1["dR10<->12"] <- md1a$R2 - md0a$R2

946 out1["dLL10<->12"] <- md1a$LogLik - md0a$LogLik

947 }

948 }

949 

950 ## One new dataset from M11

951 if (!is(md11, "try-error")) {

952 muhat <- coef(md11)["intercept"] + coef(md11)["NearField"]*df$NearField +

953 coef(md11)["m.luR"] * df$ls.R + coef(md11)["m.luU"] * df$ls.U +

954 coef(md11)["m.lV"] * df$lV + coef(md11)["m.lV.N"] * df$NearField * df$lV

955 Sigma <- lP.cov + diag(exp(2*coef(md11)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

956 Nrep <- 0

957 dorep <- TRUE

958 while(dorep && Nrep < 3){

959 Nrep <- Nrep + 1

960 lP <- rmvnorm(1, muhat, sigma=Sigma)

961 

962 df0$lP <- lP[1,]

963 md0a <- mvtnormfit(df0, 11, lP.cov)

964 md1a <- mvtnormfit(df0, 12, lP.cov)

965 }

966 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

967 out1["dR11<->12"] <- md1a$R2 - md0a$R2

968 out1["dLL11<->12"] <- md1a$LogLik - md0a$LogLik

969 }

970 }

971 

972 ## dataset from M16

973 if (!is(md16, "try-error")){

974 muhat <- coef(md16)["intercept"] +

975 df$NearField * (coef(md16)["NearField"] +

976 coef(md16)["m.luR.N"] * df$ls.R +

977 coef(md16)["m.luU.N"] * df$ls.U +

978 coef(md16)["m.lV.N"] * df$lV)

979 Sigma <- lP.cov +

980 diag(exp(2*coef(md16)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

981 

982 Nrep <- 0

983 dorep <- TRUE

984 while(dorep && Nrep < 3){

985 Nrep <- Nrep + 1

986 lP <- rmvnorm(1, muhat,

987 sigma=Sigma)

988 df0$lP <- lP[1,]

989 md0a <- mvtnormfit(df0, 16, lP.cov)

990 md1a <- mvtnormfit(df0, 12, lP.cov)

991 }

992 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

993 out1["dR16<->12"] <- md1a$R2 - md0a$R2

994 out1["dLL16<->12"] <- md1a$LogLik - md0a$LogLik

995 }

996 }

997 

998 ## dataset from M18

999 if (!is(md18, "try-error")){

1000 muhat <- coef(md18)["intercept"] +

1001 df$NearField * (coef(md18)["NearField"] +

1002 coef(md18)["m.luR.N"] * df$ls.R +

1003 coef(md18)["m.lV.N"] * df$lV)
1004 Sigma <- lP.cov +

1005 diag(exp(2*coef(md18)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

1006 

1007 Nrep <- 0

1008 dorep <- TRUE

1009 while(dorep && Nrep < 3){

1010 Nrep <- Nrep + 1

1011 lP <- rmvnorm(1, muhat,

1012 sigma=Sigma)

1013 df0$lP <- lP[1,]

1014 md0a <- mvtnormfit(df0, 18, lP.cov)

1015 md1a <- mvtnormfit(df0, 12, lP.cov)

1016 }

1017 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

1018 out1["dR18<->12"] <- md1a$R2 - md0a$R2

1019 out1["dLL18<->12"] <- md1a$LogLik - md0a$LogLik

1020 }

1021 }

1022 

1023 ## dataset from M20

1024 if (!is(md20, "try-error")){

1025 muhat <- coef(md20)["intercept"] +

1026 df$NearField * (coef(md20)["NearField"] +

1027 coef(md20)["m.luR.N"] * df$ls.R)

1028 Sigma <- lP.cov +

1029 diag(exp(2*coef(md20)["lsigma"]), nrow(lP.cov), ncol(lP.cov))

1030 

1031 Nrep <- 0

1032 dorep <- TRUE

1033 while(dorep && Nrep < 3){

1034 Nrep <- Nrep + 1

1035 lP <- rmvnorm(1, muhat,

1036 sigma=)

1037 df0$lP <- lP[1,]

1038 md0a <- mvtnormfit(df0, 20, lP.cov)

1039 md1a <- mvtnormfit(df0, 12, lP.cov)

1040 }

1041 if (!is(md0a, "try-error") && !is(md1a,"try-error")) {

1042 out1["dR20<->12"] <- md1a$R2 - md0a$R2

1043 out1["dLL20<->12"] <- md1a$LogLik - md0a$LogLik

1044 }

1045 }

1046

1047 out1

1048 }

1049

1050 analysisC <- cmpfun(analysis)

1051

1052 df <- data.frame(NearField = pred.data$NearField,

1053 lR = log(pred.data$s.R) + lV.mx[1,],

1054 lU = log(pred.data$s.U) + lV.mx[1,],

1055 ls.R = log(pred.data$s.R),

1056 ls.U = log(pred.data$s.U),

1057 CAS = pred.data$CAS,

1058 Name = pred.data$Name,

1059 MW = pred.data$MW,

1060 lP = lP.mx[1,],

1061 lV = lV.mx[1,],

1062 vars = lP.vars)

1063

1064 out <- mclapply(RegSampIndx, function(j) {

1065 df$lV <- lV.mx[j,]

1066 df$lP <- lP.mx[j,]
1067 analysisC(df)

1068 },mc.preschedule=FALSE, mc.cores=6)

1069 Betas <- do.call("rbind", out)

1070 save(Betas, file="AnalysisOmnibusFitssave.RData")

1071 }

1072

1073 if (!FittingModels) load("AnalysisOmnibusFitssave.RData")

1074

1075 doSummarizeAnalysis <- function(Betas) {

1076

1077 ## 

1078 ## Which model works best overall?

1079 ## Truncate Betas to the first 1500 values. The rest are all NA

1080 AICs <- Betas[,grep("AIC",colnames(Betas))]

1081

1082

1083 ## Summary of the Betas. Primarily of interest to track NA's.

1084 writeLines(c("",

1085 

1086 "========================================================================",

1087 "Summary of All Model Fits"))

1088 print(summary(Betas))

1089 

1090 ## Mean AIC

1091 writeLines(c("",

1092 

1093 "========================================================================",

1094 "Mean AIC"))

1095 z <- apply(AICs,2,mean, na.rm=TRUE)

1096 print(z - z[which.min(z)], digits=3)

1097

1098 writeLines(c("",

1099 

1100 "==============================================================================",

1101 "Distribution of which models have best AIC"))

1102

1103 print(table(apply(AICs,1,which.min)-1))

1104

1105

1106 ##

1107 ## Another way to look at this is with mean ranks: Rank

1108 ## the AICs across models within replicates, then

1109 ## average the ranks across replicates. This will let

1110 ## us see, for instance, if a given model is often

1111 ## second-best

1112

1113 writeLines(c("",

1114 

1115 "=============================================================================",

1116 "Mean Rank of each model's AIC"))

1117

1118 tmp <- apply(AICs,1,rank, na.last="keep")

1119 mnrank <- apply(tmp,1,mean, na.rm=TRUE)

1120 print(mnrank)

1121

1122

1123

1124 AICsummaries <-

1125 rbind(

1126 "0 <-> 1" = mysummary(AICs[,"AIC_0"] - AICs[,"AIC_1"]),

1127 "0 <-> 21" = mysummary(AICs[,"AIC_0"] - AICs[,"AIC_21"]),

1128 "0 <-> 22" = mysummary(AICs[,"AIC_0"] - AICs[,"AIC_22"]),

1129 "1 <-> 2" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_2"]),
1130 "1 <-> 3" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_3"]),

1131 "1 <-> 4" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_4"]),

1132 "1 <-> 5" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_5"]),

1133 "1 <-> 6" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_6"]),

1134 "1 <-> 12" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_12"]),

1135 "2 <-> 4" = mysummary(AICs[,"AIC_2"] - AICs[,"AIC_4"]),

1136 "2 <-> 6" = mysummary(AICs[,"AIC_2"] - AICs[,"AIC_6"]),

1137 "3 <-> 4" = mysummary(AICs[,"AIC_3"] - AICs[,"AIC_4"]),

1138 "4 <-> 7" = mysummary(AICs[,"AIC_4"] - AICs[,"AIC_7"]),

1139 "5 <-> 6" = mysummary(AICs[,"AIC_5"] - AICs[,"AIC_6"]),

1140 "6 <-> 7" = mysummary(AICs[,"AIC_6"] - AICs[,"AIC_7"]),

1141 "7 <-> 12" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_12"]),

1142 "7 <-> 9" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_9"]),

1143 "7 <-> 10" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_10"]),

1144 "7 <-> 11" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_11"]),

1145 "8 <-> 7" = mysummary(AICs[,"AIC_8"] - AICs[,"AIC_7"]),

1146 "9 <-> 12" = mysummary(AICs[,"AIC_9"] - AICs[,"AIC_12"]),

1147 "10 <-> 12" = mysummary(AICs[,"AIC_10"] - AICs[,"AIC_12"]),

1148 "11 <-> 12" = mysummary(AICs[,"AIC_11"] - AICs[,"AIC_12"]),

1149 "16 <-> 12" = mysummary(AICs[,"AIC_16"] - AICs[,"AIC_12"]),

1150 "18 <-> 12" = mysummary(AICs[,"AIC_18"] - AICs[,"AIC_12"]),

1151 "20 <-> 12" = mysummary(AICs[,"AIC_20"] - AICs[,"AIC_12"])

1152 )

1153

1154

1155 writeLines(c("",

1156 

1157 "=============================================================================",

1158 "Distributions of AIC differences for Model Comparisons"))

1159 print(AICsummaries, digits=3) 

1160

1161 ## Test Results. Compare log-likelihood differences with

1162 ## differences generated from the simpler model for all these

1163 ## comparisons. Tabulate the fraction of replicates for which

1164 ## LL_j - LL_k > dLLjk

1165

1166 ## Also, What is the change in R2 for each treatment? Put them both in

1167 ## the same table

1168

1169

1170

1171 Tests <-

1172 matrix(c(

1173 mycomps(Betas,c(0,1)),

1174 mycomps(Betas,c(0,21)),

1175 mycomps(Betas,c(0,22)),

1176 mycomps(Betas,c(1,2)),

1177 mycomps(Betas,c(1,3)),

1178 mycomps(Betas,c(1,4)),

1179 mycomps(Betas,c(1,5)),

1180 mycomps(Betas,c(1,6)),

1181 mycomps(Betas,c(1,12)),

1182 mycomps(Betas,c(2,4)),

1183 mycomps(Betas,c(2,6)),

1184 mycomps(Betas,c(3,4)),

1185 mycomps(Betas,c(4,7)),

1186 mycomps(Betas,c(5,6)),

1187 mycomps(Betas,c(6,7)),

1188 mycomps(Betas,c(7,9)),

1189 mycomps(Betas,c(7,10)),

1190 mycomps(Betas,c(7,11)),

1191 mycomps(Betas,c(7,12)),

1192 mycomps(Betas,c(8,7)),
1193 mycomps(Betas,c(9,12)),

1194 mycomps(Betas,c(10,12)),

1195 mycomps(Betas,c(11,12)),

1196 mycomps(Betas,c(16,12)),

1197 mycomps(Betas,c(18,12)),

1198 mycomps(Betas,c(20,12))

1199 ),

1200 ncol=2,

1201 dimnames=list(c(

1202 "0 <-> 1","0 <-> 21","0 <-> 22",

1203 "1 <-> 2","1 <-> 3","1 <-> 4","1 <-> 5","1 <-> 6","1 <-> 12",

1204 "2 <-> 4","2 <-> 6",

1205 "3 <-> 4",

1206 "4 <-> 7",

1207 "5 <-> 6",

1208 "6 <-> 7",

1209 "7 <-> 9","7 <-> 10","7 <-> 11","7 <-> 12",

1210 "8 <-> 7",

1211 "9 <-> 12","10 <-> 12","11 <-> 12",

1212 "16 <-> 12","18 <-> 12","20 <-> 12"

1213 ),

1214 c("P-value","delta R2")),

1215 byrow=TRUE)

1216

1217

1218 writeLines(c("",

1219 

1220 "=============================================================================",

1221 "Significance of Model Comparisons"))

1222 print(Tests)

1223

1224

1225 ## -------------------------------------------------------------------------------

1226 ## Estimates

1227

1228 writeLines(c("",

1229 

1230 "=============================================================================",

1231 "Estimates of Model Parameters, with 95% CI",""))

1232 for (i in 1:22)

1233 printests(i, Betas)

1234

1235 }

1236

1237 try(doSummarizeAnalysis(Betas))
QFunctions.R: Supplemental Functions for Previous Section

1 ## It seems natural to use gls() to do the fitting, but I cannot

2 ## figure out how to get it to handle the covariance structure I need

3 ## for this problem. The model is lP = f(predictions, parameters) +

4 ## epsilon + zeta epsilon is normal(0, sigma^2), and zeta is MVnorm(0,

5 ## defined covariance), where the defined covariance is estimated from

6 ## the MCMC samples. To carry this out, I make a new fitting function

7 ## mvtnormfit(), which takes as arguments:

8 ## y

9 ## x1,...,x3 (not all required)

10 ## M (model number, according to the above scheme)

11 ## lPcov: the defined covariance matrix

12 ## Return value is:

13 ## loglik: log likelihood

14 ## AIC: AIC

15 ## parameters: named appropriately

16 ## R2: defined as 1 - var(y - pred) / var(y)

17

18 ## getparms just returns a named parameter vector for the specified model

19

20 getparms <- function(M) {

21 switch(M+1,

22 ## ------- 0

23 structure(rep(as.numeric(NA),2), names=c("intercept","lsigma")),

24 ## ------- 1

25 structure(rep(as.numeric(NA),3),

26 names=c("intercept","NearField","lsigma")),

27 ## ------- 2

28 structure(rep(as.numeric(NA),4),

29 names=c("intercept","NearField","m.lV","lsigma")),

30 ## ------- 3

31 structure(rep(as.numeric(NA),4),

32 names=c("intercept","NearField","m.lR","lsigma")),

33 ## ------- 4

34 structure(rep(as.numeric(NA),5),

35 names=c("intercept","NearField","m.luR","m.lV",

36 "lsigma")),

37 ## ------- 5

38 structure(rep(as.numeric(NA),4),

39 names=c("intercept","NearField","m.lU","lsigma")),

40 ## ------- 6

41 structure(rep(as.numeric(NA),5),

42 names=c("intercept","NearField","m.luU","m.lV","lsigma")),

43 ## ------- 7

44 structure(rep(as.numeric(NA),6),

45 names=c("intercept","NearField","m.luR","m.luU","m.lV",

46 "lsigma")),

47 ## ------- 8

48 structure(rep(as.numeric(NA),5),

49 names=c("intercept","NearField","m.lR","m.lU","lsigma")),

50 ## ------- 9

51 structure(rep(as.numeric(NA),7),

52 names=c("intercept","NearField",

53 "m.luR","m.luU","m.lV",

54 "m.luR.N",

55 "lsigma")),

56 ## ------- 10

57 structure(rep(as.numeric(NA),7),

58 names=c("intercept","NearField",

59 "m.luR","m.luU","m.lV",

60 "m.luU.N",
61 "lsigma")),

62 ## ------- 11

63 structure(rep(as.numeric(NA),7),

64 names=c("intercept","NearField",

65 "m.luR","m.luU","m.lV",

66 "m.lV.N",

67 "lsigma")),

68 ## ------- 12

69 structure(rep(as.numeric(NA),9),

70 names=c("intercept","NearField",

71 "m.luR","m.luU","m.lV",

72 "m.luR.N","m.luU.N","m.lV.N",

73 "lsigma")),

74 ## ------- 13

75 structure(rep(as.numeric(NA),8),

76 names=c("intercept","NearField",

77 "m.luR","m.luU","m.lV",

78 "m.luR.N","m.luU.N",

79 "lsigma")),

80 ## ------- 14

81 structure(rep(as.numeric(NA),8),

82 names=c("intercept","NearField",

83 "m.luR","m.luU","m.lV",

84 "m.luU.N","m.lV.N",

85 "lsigma")),

86 ## ------- 15

87 structure(rep(as.numeric(NA),8),

88 names=c("intercept","NearField",

89 "m.luR","m.luU","m.lV",

90 "m.luR.N","m.lV.N",

91 "lsigma")),

92 ## ------- 16

93 structure(rep(as.numeric(NA),6),

94 names=c("intercept","NearField",

95 "m.luR.N","m.luU.N","m.lV.N",

96 "lsigma")),

97 ## ------- 17

98 structure(rep(as.numeric(NA),5),

99 names=c("intercept","NearField",

100 "m.luR.N","m.luU.N",

101 "lsigma")),

102 ## ------- 18

103 structure(rep(as.numeric(NA),5),

104 names=c("intercept","NearField",

105 "m.luR.N","m.lV.N",

106 "lsigma")),

107 ## ------- 19

108 structure(rep(as.numeric(NA),5),

109 names=c("intercept","NearField",

110 "m.luU.N","m.lV.N",

111 "lsigma")),

112 ## ------- 20

113 structure(rep(as.numeric(NA),4),

114 names=c("intercept","NearField",

115 "m.luR.N",

116 "lsigma")),

117 ## ------- 21

118 structure(rep(as.numeric(NA),3),

119 names=c("intercept","m.lR",

120 "lsigma")),

121 ## ------- 22

122 structure(rep(as.numeric(NA),3),

123 names=c("intercept","m.lU",
124 "lsigma"))

125 )

126 }

127 mvtnormfit <- function(df, M, lPcov) {

128 parms <- getparms(M)

129 n <- nrow(df)

130

131 mufun <- with(c(as.list(df),list(n=n)),

132 switch(M+1,

133 ## ------ M0

134 function(p){

135 rep(p[1],n)

136 },

137 

138 ## ------ M1

139 function(p){

140 p[1] + p[2] * NearField

141 },

142 ## ------ M2

143 function(p){

144 p[1] + p[2] * NearField + p[3] * lV

145 },

146 ## ------ M3

147 function(p){

148 p[1] + p[2] * NearField + p[3] * lR

149 },

150 ## ------ M4

151 function(p){

152 p[1] + p[2] * NearField + p[3] * ls.R + p[4] * lV

153 },

154 ## ------ M5

155 function(p){

156 p[1] + p[2] * NearField + p[3] * lU

157 },

158 ## ------ M6

159 function(p){

160 p[1] + p[2] * NearField + p[3] * ls.U + p[4] * lV

161 },

162 ## ------ M7

163 function(p){

164 p[1] + p[2] * NearField + p[3] * ls.R + p[4] * ls.U +

165 p[5] * lV

166 },

167 ## ------ M8

168 function(p){

169 p[1] + p[2] * NearField + p[3] * lR + p[4] * lU

170 },

171 ## ------ M9

172 function(p){

173 p[1] + p[2] * NearField +

174 p[3] * ls.R +

175 p[4] * ls.U +

176 p[5] * lV +

177 p[6] * NearField * ls.R

178 },

179 ## ------ M10

180 function(p){

181 p[1] + p[2] * NearField +

182 p[3] * ls.R +

183 p[4] * ls.U +

184 p[5] * lV +

185 p[6] * NearField * ls.U

186 },
187 ## ------ M11

188 function(p){

189 p[1] + p[2] * NearField +

190 p[3] * ls.R +

191 p[4] * ls.U +

192 p[5] * lV +

193 p[6] * NearField * lV

194 },

195 ## ------ M12

196 function(p){

197 p[1] + p[2] * NearField +

198 p[3] * ls.R +

199 p[4] * ls.U +

200 p[5] * lV +

201 p[6] * NearField * ls.R +

202 p[7] * NearField * ls.U +

203 p[8] * NearField * lV

204 },

205 ## ------ M13

206 function(p){

207 p[1] + p[2] * NearField +

208 p[3] * ls.R +

209 p[4] * ls.U +

210 p[5] * lV +

211 p[6] * NearField * ls.R +

212 p[7] * NearField * ls.U

213 },

214 ## ------ M14

215 function(p){

216 p[1] + p[2] * NearField +

217 p[3] * ls.R +

218 p[4] * ls.U +

219 p[5] * lV +

220 p[6] * NearField * ls.U +

221 p[7] * NearField * lV

222 },

223 ## ------ M15

224 function(p){

225 p[1] + p[2] * NearField +

226 p[3] * ls.R +

227 p[4] * ls.U +

228 p[5] * lV +

229 p[6] * NearField * ls.R +

230 p[7] * NearField * lV

231 },

232 ## ------ M16

233 function(p){

234 p[1] + p[2] * NearField +

235 p[3] * NearField * ls.R +

236 p[4] * NearField * ls.U +

237 p[5] * NearField * lV

238 },

239 ## ------ M17

240 function(p){

241 p[1] + p[2] * NearField +

242 p[3] * NearField * ls.R +

243 p[4] * NearField * ls.U

244 },

245 ## ------ M18

246 function(p){

247 p[1] + p[2] * NearField +

248 p[3] * NearField * ls.R +

249 p[4] * NearField * lV
250 },

251 ## ------ M19

252 function(p){

253 p[1] + p[2] * NearField +

254 p[3] * NearField * ls.U +

255 p[4] * NearField * lV

256 },

257 ## ------ M20

258 function(p){

259 p[1] + p[2] * NearField +

260 p[3] * NearField * ls.R

261 },

262 ## ------ M21

263 function(p){

264 p[1] + p[2] * lR

265 },

266 ## ------ M22

267 function(p){

268 p[1] + p[2] * lU

269 },

270

271 ))

272 of <- function(p) {

273 mu <- mufun(p)

274 dmvnorm(df$lP, mu, lPcov + diag(exp(p[length(p)]*2), nrow(df), nrow(df)),

275 log=TRUE)

276 }

277

278 ## For these models, it works better if NearField is a factor

279 df0 <- df

280 df0$NearFieldn <- df0$NearField

281 df0$NearField <- factor(df0$NearField)

282 ## Get initial value of p. Use gls()

283 out <- switch(M+1,

284 ## ------ M0

285 gls(lP ~ 1,

286 data=df0, 

287 correlation=corSymm(value=corparms, fixed=TRUE),

288 weights=varFixed(value = ~vars), method="ML"),

289 ## ------ M1

290 gls(lP ~ NearField,

291 data=df0, 

292 correlation=corSymm(value=corparms, fixed=TRUE),

293 weights=varFixed(value = ~vars), method="ML"),

294 ## ------ M2

295 gls(lP ~ NearField + lV,

296 data=df0, 

297 correlation=corSymm(value=corparms, fixed=TRUE),

298 weights=varFixed(value = ~vars), method="ML"),

299 ## ------ M3

300 gls(lP ~ NearField + lR,

301 data=df0, 

302 correlation=corSymm(value=corparms, fixed=TRUE),

303 weights=varFixed(value = ~vars), method="ML"),

304 ## ------ M4

305 gls(lP ~ NearField + ls.R + lV,

306 data=df0, 

307 correlation=corSymm(value=corparms, fixed=TRUE),

308 weights=varFixed(value = ~vars), method="ML"),

309 ## ------ M5

310 gls(lP ~ NearField + lU,

311 data=df0, 

312 correlation=corSymm(value=corparms, fixed=TRUE),
313 weights=varFixed(value = ~vars), method="ML"),

314 ## ------ M6

315 gls(lP ~ NearField + ls.U + lV,

316 data=df, 

317 correlation=corSymm(value=corparms, fixed=TRUE),

318 weights=varFixed(value = ~vars), method="ML"),

319 ## ------ M7

320 gls(lP ~ NearField + ls.R + ls.U + lV,

321 data=df0, 

322 correlation=corSymm(value=corparms, fixed=TRUE),

323 weights=varFixed(value = ~vars), method="ML"),

324 ## ------ M8

325 gls(lP ~ NearField + lR + lU,

326 data=df0, 

327 correlation=corSymm(value=corparms, fixed=TRUE),

328 weights=varFixed(value = ~vars), method="ML"),

329 ## ------ M9

330 gls(lP ~ NearField * ls.R + ls.U + lV,

331 data=df0, 

332 correlation=corSymm(value=corparms, fixed=TRUE),

333 weights=varFixed(value = ~vars), method="ML"),

334 ## ------ M10

335 gls(lP ~ NearField + ls.R + NearField * ls.U + lV,

336 data=df0, 

337 correlation=corSymm(value=corparms, fixed=TRUE),

338 weights=varFixed(value = ~vars), method="ML"),

339 ## ------ M11

340 gls(lP ~ NearField + ls.R + ls.U + NearField * lV,

341 data=df0, 

342 correlation=corSymm(value=corparms, fixed=TRUE),

343 weights=varFixed(value = ~vars), method="ML"),

344 ## ------ M12

345 gls(lP ~ NearField * (ls.R + ls.U + lV),

346 data=df0, 

347 correlation=corSymm(value=corparms, fixed=TRUE),

348 weights=varFixed(value = ~vars), method="ML"),

349 ## ------ M13

350 gls(lP ~ NearField + NearField * (ls.R + ls.U) + lV,

351 data=df0, 

352 correlation=corSymm(value=corparms, fixed=TRUE),

353 weights=varFixed(value = ~vars), method="ML"),

354 ## ------ M14

355 gls(lP ~ NearField + ls.R + NearField * (ls.U + lV),

356 data=df0, 

357 correlation=corSymm(value=corparms, fixed=TRUE),

358 weights=varFixed(value = ~vars), method="ML"),

359 ## ------ M15

360 gls(lP ~ NearField + NearField * ls.R + ls.U + NearField * lV,

361 data=df0, 

362 correlation=corSymm(value=corparms, fixed=TRUE),

363 weights=varFixed(value = ~vars), method="ML"),

364 ## ------ M16

365 gls(lP ~ NearField + I(NearFieldn*ls.R) +

366 I(NearFieldn*ls.U) + I(NearFieldn*lV),

367 data=df0, 

368 correlation=corSymm(value=corparms, fixed=TRUE),

369 weights=varFixed(value = ~vars), method="ML"),

370 ## ------ M17

371 gls(lP ~ NearField + I(NearFieldn*ls.R) + I(NearFieldn*ls.U),

372 data=df0, 

373 correlation=corSymm(value=corparms, fixed=TRUE),

374 weights=varFixed(value = ~vars), method="ML"),

375 ## ------ M18
376 gls(lP ~ NearField + I(NearFieldn*ls.R) + I(NearFieldn*lV),

377 data=df0, 

378 correlation=corSymm(value=corparms, fixed=TRUE),

379 weights=varFixed(value = ~vars), method="ML"),

380 ## ------ M19

381 gls(lP ~ NearField + I(NearFieldn*ls.U) + I(NearFieldn*lV),

382 data=df0, 

383 correlation=corSymm(value=corparms, fixed=TRUE),

384 weights=varFixed(value = ~vars), method="ML"),

385 ## ------ M20

386 gls(lP ~ NearField + I(NearFieldn*ls.R),

387 data=df0, 

388 correlation=corSymm(value=corparms, fixed=TRUE),

389 weights=varFixed(value = ~vars), method="ML"),

390 ## ------ M21

391 gls(lP ~ lR,

392 data=df0, 

393 correlation=corSymm(value=corparms, fixed=TRUE),

394 weights=varFixed(value = ~vars), method="ML"),

395 ## ------ M22

396 gls(lP ~ lU,

397 data=df0, 

398 correlation=corSymm(value=corparms, fixed=TRUE),

399 weights=varFixed(value = ~vars), method="ML"),

400 )

401

402 SS <- median((out$sigma^2 - 1)*diag(lPcov))

403 pstart <- c(coef(out),log(SS)/2)

404 out2 <- try(optim(pstart, of, method="BFGS",control=list(fnscale=-1)))

405 ## If try-error, first try using simplex method

406 if (is(out2, "try-error")) {

407 out2 <- try(optim(pstart, of, control=list(fnscale=-1, maxit=2000)))

408 }

409 if (!is(out2, "try-error")) {

410 muhat <- with(df, mufun(out2$par))

411 par <- structure(out2$par, names=names(parms))

412 list(LogLik = out2$value,

413 ## AIC = -2*(out2$value) -2*length(pstart), WRONG! Corrected 7/5/2012

414 ## RWS

415 AIC = -2*(out2$value) + 2*length(pstart),

416 R2 = 1 - var(df$lP - muhat)/var(df$lP),

417 coefficients = par)

418 } else {

419 out2

420 }

421 }

422

423

424 ## Distributions of AIC values across replicate samples

425 mysummary <- function(x){

426 qq <- stats::quantile(x, pr=c(0,0.025,0.25,0.5,0.75,0.975,1), na.rm=TRUE)

427 c(qq[1:4],mean=mean(x, na.rm=TRUE),qq[5:7])

428 }

429

430 ## Changes 7/5/2012: separator for dll is now "<->" RWS

431 ## Also, comparison between delta LL and null delta LL was wrong. RWS

432 mycomps <- function(z, comp=c(0,1),NEW = TRUE) {

433 dll <- if (NEW) {

434 paste("dLL",comp[1],"<->",comp[2],sep="")

435 } else {

436 paste("dLL",comp[1],comp[2],sep="")

437 }

438 LLsimp <- paste("LL",comp[1],sep="_")
439 LLcompl <- paste("LL",comp[2],sep="_")

440 R2simp <- paste("R2",comp[1],sep="_")

441 R2compl <- paste("R2",comp[2],sep="_")

442 c(

443 sum(as.numeric(abs(z[,LLcompl] - z[,LLsimp]) < abs(z[,dll])))/nrow(z),

444 if(comp[1] > 0) {

445 mean(z[,R2compl] - z[,R2simp], na.rm=TRUE)

446 } else {

447 mean(z[,R2compl])

448 }

449 )

450 }

451

452 printests <- function(n, Betas) {

453 pat <- paste("_",n,"$",sep="")

454 indx <- grep(pat, colnames(Betas))

455 mns <- apply(Betas[,indx],2,mean)

456 lcl <- apply(Betas[,indx],2,quantile,pr=0.025)

457 ucl <- apply(Betas[,indx],2,quantile,pr=0.975)

458

459 out <- cbind(mns, lcl, ucl)

460 rnames <- names(mns)

461 rnames[1] <- "Intercept"

462 rownames(out) <- rnames

463 colnames(out) <- c("Mean","LCL","UCL")

464 writeLines(c("",paste("Bootstrap Distributions: Model",n)))

465 print(out)

466 }
Fitmodel16.R: R code to fit Model 16 to NHANES Parent Chemical Data

1 library(doMC)

2 registerDoMC(cores=3)

3 library(rjags)

4 load.module("lecuyer")

5 library(random)

6 library(R2jags)

7

8 ## JAGS model is in model16.jag

9 ## It needs the following data:

10 ## R.R : RAIDAR predicted exposure (1:N),

11 ## R.U : USEtox predicted exposure (1:N),

12 ## y : observed metabolite rates (NA for below LOD) (1:M),

13 ## N = number of parents,

14 ## N.nontriv = number of parents with > 1 metabolite (sorted so they come first

15 ## in R.x)

16 ## N.mets[1:N.nontriv] = number of metabolites for parent i

17 ## indx[N,M]: indx[i,1:Nmets[i]] gives the indices in U for the metabolites

18 ## i contributes. Corresponding values of Phi[i,] are non-zero.

19 ## If i > N.nontriv, indx[i,1] gives the entry in Phi[i,] that is 1

20 ## For j > N.nontriv, indx[i,j] gives column indices in Phi[i,] for 0's.

21 ## M = number of metabolites,

22 ## tau[1:M] : precision for observed metabolite rates.

23 ## lod = vector of limits of detection, tau = vector of precisions for observed

24 ## metabolite rates.

25 ## dstart[1:N.nontriv], dstop[1:N.nontriv]: indices into delta and phi

26 ## such that, e.g., delta[dstart[i]:dstop[i]] are the deltas for row i of Phi.

27 ## dstart[1] == 1, dstop[N.nontriv] == Ndelta.

28 ##

29 ## The following variables need initializing

30 ## P[i], i in 1:N

31 ## a

32 ## m.R

33 ## m.U

34 ## X1

35 ## X2

36 ## delta[1:Ndelta]: the non-trivial deltas corresponding to the non-trivial phis

37

38 ## Data are based on John W.'s RData file DataTables-040212.RData, and successors.

39

40 attach(file.path("..","MCMCdata.RData"))

41 z <- runif(1)

42 doONEchain <- TRUE

43 if (doONEchain) { ## this is a restart; use the chains we already have

44 ## to initialize stuff.

45

46 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

47 ## JAGS inputs: nhanesdata and inits

48 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

49

50 nhanesdata <- list()

51 ## The Jags models need these data:

52 ## IURcuts - cutpoints for the IUR production volumes (cutpoints above)

53 nhanesdata$IURcuts <- log(cutpoints)

54 ## IURindx - vector indexing the rows of pred.data with IUR production volumes

55 nhanesdata$IURindx <- grep("IUR HPV",pred.data$PUsource)

56 ## IURvols - indicator for which interval an IUR production volume falls

57 nhanesdata$IURvols <-

58 as.numeric(pred.data$V.cat[nhanesdata$IURindx[1:length(nhanesdata$IURindx)]])

59 ## CPRIindx indexes the rows of pred.data with CPRI use volumes

60 nhanesdata$CPRIindx <- which(pred.data$PUsource == "CPRI 2002")
61 nhanesdata$lV.CPRI <- log(pred.data[nhanesdata$CPRIindx,"V"])

62 ## ls.R, ls.U are the log unit predictions from RADAR and USEtox, respectively

63 nhanesdata$ls.R <- log(pred.data$s.R)

64 nhanesdata$ls.U <- log(pred.data$s.U)

65 nhanesdata$NearField <- pred.data$NearField

66

67 ## ly - the log of the observed metabolite quantiles

68 nhanesdata$ly <-

69 ly <-

70 t(log(data.matrix(Measured[,c("Value50","Value75","Value90","Value95")])))

71 ## Q, the proportions associated with the quantiles

72 nhanesdata$Q <- c(0.5,0.75,0.9,0.95)

73 ## se - precision of measured metabolite values, from reported confidence limits on 

74 quantiles

75 nhanesdata$se <- Measured$se

76 ## llod - log(Limit of Detection)

77 nhanesdata$llod <- log(Measured$LOQ)

78 ## is.above.LOD - vector with 0 for below lod, 1 for above it

79 nhanesdata$is.above.LOD <- 0 + !is.na(ly)

80 ## Some counts:

81 ## N is the number of parents

82 nhanesdata$N <- nrow(pred.data)

83 ## M is the number of metabolites

84 nhanesdata$M <- nrow(Measured)

85 ## N.IUR is the number of IUR reported production volumes (length of IURindx)

86 nhanesdata$N.IUR <- length(nhanesdata$IURindx)

87 ## N.CPRI is the number of CPRI reported use volumes (length of CPRIindx)

88 nhanesdata$N.CPRI <- length(nhanesdata$CPRIindx)

89 ## N.nontriv is the number of parents that result in >1 metabolite

90 nhanesdata$N.nontriv <- sum(N.metabs > 1)

91 ## Nmets - number of metabolites for parent[i], when > 1

92 nhanesdata$Nmets <- N.metabs[N.metabs > 1]

93 ## index - index[i,j] is the index into the metabolite list for the jth metabolite

94 ## of the ith parent

95 nhanesdata$indx <-

96 t(apply(Xref,1,function(x)c((1:length(x))[is.na(x)],

97 (1:length(x))[!is.na(x)])))

98 ## begin and end index into phi[] for each row of Phi

99 tmp <- cumsum(nhanesdata$Nmets)

100 nhanesdata$dstart <- 1+c(0,tmp[-nhanesdata$N.nontriv])

101 nhanesdata$dstop <- tmp

102

103 ## Phi

104 nhanesdata$Phi <- Xref

105 nhanesdata$Ndelta <- sum(is.na(Xref))

106

107 ## ---------------------------------------------------------------------

108 ## Initialization

109

110 inits <- list(list())

111 inits[[1]]$lV.IUR <- log(pred.data$V[nhanesdata$IURindx])

112 inits[[1]]$sdlVmu.dum <- sd(log(pred.data$V))

113 inits[[1]]$lVmu <- mean(log(pred.data$V))

114 inits[[1]]$lV.CPRI.mu <- log(pred.data$V[nhanesdata$CPRIindx])

115 inits[[1]]$lP <- rep(mean(nhanesdata$ly,na.rm=TRUE), nrow(pred.data))

116

117 inits[[1]]$a.F <- -7 ## mean(nhanesdata$ly, na.rm=TRUE)

118 inits[[1]]$a.N <- inits[[1]]$a.F

119 inits[[1]]$m.R.N <- -1.2

120 inits[[1]]$m.U.N <- 0.60

121 inits[[1]]$m.V.N <- 0.20

122 inits[[1]]$sd.dum <- 3

123 inits[[1]]$sd <- rep(1.0, ncol(ly))
124 ## for the censored values in ly

125 ily <- is.na(ly)

126 loqmx <- outer(rep(1,4),log(Measured$LOQ),"*")

127 lyinit <- ifelse(ily,loqmx-2,NA)

128 inits[[1]]$ly <- lyinit

129 inits[[1]]$lsd <- 1

130 inits[[1]]$lsd.dum <- 2

131

132 inits[[1]]$delta <- rep(1, nhanesdata$Ndelta)

133

134

135 doinits <- function(chain) {inits[[chain]]}

136

137 assign("inits", inits, envir=environment(doinits))

138

139 ## JAGS code

140

141 ## model <- jags.model("model3.jag",data=nhanesdata,inits=inits[[1]],

142 ## n.chains=1,n.adapt=1000)

143

144

145 ## Now, run for 300000 samples, thinned every 100, and dump lP, lU, and delta

146 ## as well. We will use these to get starting values for three new chains.

147

148 out.samples <- jags(data=nhanesdata, inits=doinits,

149 parameters.to.save=c("a.F","a.N",

150 "m.R.N","m.U.N","m.V.N",

151 "lV","sd.V","lP","phi","lU","sd.lVmu","lVmu","Z",

152 "sd","lsd","lsd.sd"),

153 model.file="model16.jag",

154 n.chains=1, n.iter=3000, n.burnin=1000, n.thin=2,

155 DIC=FALSE,progress.bar="text")

156

157

158 out.coda <- as.mcmc(out.samples)

159 nms <- varnames(out.coda)

160 for (nm in nms) {

161 writeLines(c("","-----------------",nm,""))

162 print(try(heidel.diag(out.coda[,nm])))

163 print(try(summary(out.coda[,nm])))

164 }

165

166 save(out.samples, out.coda, nhanesdata, inits, file="m3parms1.RData")

167 }

168

169 load("m3parms1.RData")

170

171 ## --------------------------------------------------------------------------------

172 ## This is the initial, hopefully converged, chain. Now, use these to start with

173 ## and run three new chains.

174

175

176 ## getInits finds extreme values 

177 getInits <- function(X, ninits=3, P=0.05) {

178 ## X is an mcmc object, essentially a matrix.

179 dta <- scale(X)

180 ## Ultimately, we are looking for ninits indices into the elements of X

181 outinits <- numeric(ninits)

182 ## Pick ninits random directions. Project the data in dta on each direction

183 ## and pick the Pth quantile (since the original direction u was random

184 ## u and -u are equally likely).

185 for (i in 1:ninits) {

186 u <- rnorm(ncol(dta))
187 u <- u/sqrt(sum(u*u))

188 lngt <- dta %*% u

189 indx <- seq(along=lngt)

190 keep <- round(P*length(indx)) ## this works fine as long as length(indx) > 20

191 oind <- order(lngt)

192 outinits[i] <- oind[keep]

193 }

194 outinits

195 }

196 ## Get 3 hyperdispersed points.

197 Indx <- getInits(out.coda,3)

198

199

200

201 getsamps <- function(nm) {

202 nm <- gsub("\\.","\\\\.",nm)

203 cindx <- grep(nm, colnames(out.coda))

204 out.coda[,cindx,drop=FALSE]

205 }

206

207 inits2 <- list(list(),list(),list())

208 for (i in 1:3) {

209 inits2[[i]]$lV.IUR <- getsamps("^lV\\[")[Indx[i],nhanesdata$IURindx]

210 inits2[[i]]$sdlVmu.dum <- getsamps("^sd.lVmu")[Indx[i],]

211 inits2[[i]]$lVmu <- getsamps("^lVmu")[Indx[i],]

212 inits2[[i]]$lV.CPRI.mu <- getsamps("^lV\\[")[Indx[i],nhanesdata$CPRIindx]

213 inits2[[i]]$lP <- getsamps("^lP\\[")[Indx[i],]

214 ##inits2[[i]]$lP <- inits[[1]]$lP

215 inits2[[i]]$a.F <- getsamps("^a.F")[Indx[i],]

216 ##inits2[[i]]$a.F <- inits[[1]]$a.F

217 inits2[[i]]$a.N <- getsamps("^a.N")[Indx[i],]

218 inits2[[i]]$m.R.N <- getsamps("^m.R.N")[Indx[i],]

219 inits2[[i]]$m.U.N <- getsamps("^m.U.N")[Indx[i],]

220 inits2[[i]]$m.V.N <- getsamps("^m.V.N")[Indx[i],]

221 inits2[[i]]$sd.dum <- getsamps("^sd.V")[Indx[i],]

222 inits2[[i]]$sd <- getsamps("^sd\\[")[Indx[i],]

223 inits2[[i]]$lsd <- getsamps("^lsd$")[Indx[i],]

224 inits2[[i]]$lsd.dum <- getsamps("lsd.sd")[Indx[i],]

225 inits2[[i]]$delta <- getsamps("phi")[Indx[i],]

226 inits2[[i]]$ly<-inits[[1]]$ly

227 inits2[[i]]$.RNG.name="lecuyer::RngStream"

228 inits2[[i]]$.RNG.seed=randomNumbers(n=1, min=1, max=1e+06, col=1)

229 }

230

231 ## quickie tool to compare inits and inits2 visually

232 ## cmpin <- function(nm) {

233 ## cbind(inits[[1]][[nm]], inits2[[1]][[nm]],

234 ## inits2[[2]][[nm]], inits2[[3]][[nm]])

235 ## }

236

237 ## jags.parallel doesn't export all the global environment, so inits2

238 ## gets lost, and doinits cannot find it in their subprocesses.

239 ## The workaround used here is to create an environment for the

240 ## function doinits and explicitly assign inits2 to it. Then, when

241 ## doinits gets exported, inits2 goes with it.

242 doinits <- function(chain){inits2[[chain]]}

243 env1 <- new.env()

244 assign("inits2", inits2, pos=env1)

245 environment(doinits) <- env1

246

247 NIters <- 1000

248 Thin <- 500

249
250 out.samps3R <- foreach(i=1:3) %dopar% {

251 model <- jags.model("model16.jag", data=nhanesdata, inits=inits2[i], n.chains=1,

252 n.adapt=100000)

253 result <- coda.samples(model, variable.names=c("a.F","a.N",

254 "m.R.N","m.U.N","m.V.N",

255 "lV","sd.V","lP","lU","sd.lVmu","lVmu","Z",

256 "sd","lsd","lsd.sd"),

257 n.iter=NIters*Thin, thin=Thin)

258 return(result)

259 }

260 

261 

262

263 out.coda3R <- do.call("mcmc.list", lapply(out.samps3R, function(z) z[[1]]))

264 save(out.samps3R, out.coda3R, file="m3parms3.RData")

265 nms <- varnames(out.coda3R)

266 for (nm in nms) {

267 writeLines(c("","-----------------",nm,""))

268 print(try(heidel.diag(out.coda3R[,nm, drop=FALSE])))

269 print(try(gelman.diag(out.coda3R[,nm,drop=FALSE])))

270 print(summary(out.coda3R[,nm,drop=FALSE]))

271 }

272

273

274

275

276 
model16.jag: JAGS Code for Fitting Model 16 to NHANES Parent 

Chemical Data

1 model

2 {

3 ## Production volumes are crudely estimated or reported. We have two

4 ## sources, with their own error characteristics. The CPRI data are

5 ## estimated from reported use data, as point estimates. The EPA IUR data

6 ## are reported in predetermined bins. We treat the overall distribution

7 ## of use/production as lognormal (since everything here is on a log

8 ## scale, that means lots of normal distributions). IUR data are then

9 ## interval-censored, with predefined bins reported in IURcuts,

10 ## while each CPRI estimate is sampled from a lognormal distribution

11 ## with unknown median and coefficient of variation depending on the

12 ## estimated error of the estimate. The means are distributed according

13 ## to the same distribution as the IUR data.

14 ## Model parental dose rates as a function of RAIDAR and USEtox

15 ## predictions

16 ## ++++++++++++++++++++++++++++++++++++++++++++

17 ## Model the volumes

18 ## ++++++++++++++++++++++++++++++++++++++++++++

19 ## -------------------------------------------

20 ## IUR volumes

21 for (i in 1:N.IUR) {

22 IURvols[i] ~ dinterval(lV.IUR[i], IURcuts)

23 lV.IUR[i] ~ dnorm(lVmu, tau.lVmu)

24 lV[IURindx[i]] <- lV.IUR[i]

25 }

26 lVmu ~ dnorm(0, 0.001)

27 sdlVmu.dum ~ dt(0, 1/625, 1)

28 sd.lVmu <- abs(sdlVmu.dum)

29 tau.lVmu <- 1/pow(sd.lVmu,2)

30 ## -------------------------------------------

31 ## CPRI usage

32 for (i in 1:N.CPRI) {

33 lV.CPRI[i] ~ dnorm(lV.CPRI.mu[i], 25) ## CV of 20%: 25 <- 1/(.2^2)

34 lV.CPRI.mu[i] ~ dnorm(lVmu, tau.lVmu) 

35 lV[CPRIindx[i]] <- lV.CPRI.mu[i]

36 }

37 ## +++++++++++++++++++++++++++++++++++++++++++++

38 ## The heart of the model:

39 ## Model the parent exposures. Fit separate

40 ## coefficients for far- and nearfield chemicals

41 ## +++++++++++++++++++++++++++++++++++++++++++++

42 for (i in 1:N) {

43 Z[i] <-

44 (1 - NearField[i]) * a.F +

45 NearField[i]*(a.N + m.R.N * ls.R[i] +

46 m.U.N * ls.U[i] +

47 m.V.N * lV[i])

48 lP[i] ~ dnorm(Z[i], tau.V)

49 P[i]<- exp(lP[i])

50 }

51 ## Priors for a, m.R, m.U, and m.V

52 a.F ~ dnorm(0, 0.001)

53 a.N ~ dnorm(0, 0.001)

54 m.R.N ~ dnorm(0, 0.001)

55 m.U.N ~ dnorm(0, 0.001)

56 m.V.N ~ dnorm(0, 0.001)

57 ## code so that prior for sd(log(P[i])) is half-Cauchy(25)

58 ## Fit a separate standard deviation for far and near
59 sd.dum ~ dt(0,1/625,1)

60 sd.V <- abs(sd.dum)

61 tau.V <- 1/pow(sd.V,2)

62 ## ++++++++++++++++++++++++++++++++++++++++++++++

63 ## Link the unobserved parent exposure to the observed

64 ## (but censored) metabolite (w/some parent) urinary

65 ## values.

66 ##

67 ## Model metabolite concentrations as (usually) unknown fractions

68 ## of the parent.

69 ## 

70 ## Phi[i,j] say what fraction of parent[i] goes to metabolite[j]

71 ## observed values are lognormally distributed around U, with known

72 ## precision, but censored. Some of the elements of Phi are known to be

73 ## 0 or 1 (thanks to Jim Rabinowitz for some of the 1's).

74 for (j in 1:M) {

75 lU[j]<- log(inprod(Phi[,j],P[]))

76 for (i in 1:4) {

77 ly[i,j] ~ dnorm(qnorm(Q[i],lU[j], tau[j]),tau.se[j])

78 is.above.LOD[i,j] ~ dinterval(ly[i,j],llod[j])

79 }

80 tau.se[j] <- 1/(se[j]*se[j])

81 ## lognormal distribution for chemical-specific population sds

82 sd[j] ~ dlnorm(lsd, tau.lsd)

83 tau[j] <- 1/(sd[j]*sd[j])

84 }

85 lsd ~ dnorm(0, 0.001)

86 lsd.dum ~ dt(0,1/625,1)

87 lsd.sd <- abs(lsd.dum)

88 tau.lsd <- 1/(lsd.sd * lsd.sd)

89 ## P is sorted so that only the first N.nontriv parents result in multiple

90 ## metabolites. For i in (N.nontriv+1:N), there is only one metabolite.

91 ## Nmets[i] gives the number of metabolites of parent i

92 ## indx[i,1:Nmets[i]] give the columns of row i that are being estimated,

93 ## for i in 1:N.nontriv. When i > N.nontriv, indx[i,1] gives the 1

94 ## metabolite matched to parent i.

95 ##

96 ## The elements of phi are essentially unknown. 

97 for (k in 1:Ndelta) {

98 delta[k] ~ dunif(0,1)

99 }

100 for (i in 1:N.nontriv) {

101 for (j in 1:Nmets[i]) {

102 phi[dstart[i] + j - 1] <- delta[dstart[i] + j - 1] / 

103 sum(delta[dstart[i]:dstop[i]])

104 Phi[i,indx[i,j]]<- phi[dstart[i] + j - 1]

105 }

106 }

107 }

108