# readNHANES.R: Data Preprocessing

library(XML) 

### 1 attach(file.path("ExpoCast - 040813.RData")) 

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

### 3 NHANES.parent.metabolite.mapping 

### 4 ## In NHANES: 

### 5 ## Misspellings, etc. 

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

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

### 8 ## 

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

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

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

### 12 13 ## In NHANES.parent.metabo lite.mapping$Name.1 

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

### 15 ## 

### 16 options(stringsAsFactors=FALSE) 

### 17 18 NHANES$Name < - sub("napthalene","naphthalene",NHANES$Name) 

### 19 NHANES.parent.metabolite.mapping$Name.1 < - 

### 20 sub("napthalene","naphthalene", NHANES.parent.metab olite.mapping$Name.1) 

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

### 23 NHANES < - NHANES[ - indx,] 

### 24 25 indx < - grep("9 - Hydroxybenz \ \ [a \ \ ]anthracene", 

### 26 NHANES.parent.metabolite.mapping$Name.1) 

### 27 NHANES.parent.metabolite.mapping < - NHANES.parent.metabolit e.mapping[ - indx,] 

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

### 30 Hydroxybenz[a]anthracene",NHANES$Name) 

### 31 NHANES.parent.metabolite.mapping$Name.1 < - 

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

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

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

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

### 37 ExpoCast.NHANES.compare, 

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

### 39 40 indx < - grep("Triclosan",NHANES$Name) 

### 41 NHANES$Name[indx] < - "Triclosan" 

### 42 ## ==================================================================== 

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

### 44 ## mean production volume 

### 45 geometric.mean.to.IUR < - 

### 46 function (x) 

### 47 { 

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

### 49 1e - 14) { 

### 50 string < - "1 to < 10 million lbs" 

### 51 lower < - 1 * 1e+06/2.204/365.25 

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

### 53 } 

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

### 55 1e+06/2.204/365.25))/2))) < 

### 56 1e - 14) { 

### 57 string < - "500,000 to < 1 million lbs" 

### 58 lower < - 0.5 * 1e+06/2. 204/365.25 

### 59 upper < - 1 * 1e+06/2.204/365.25 

### 60} 

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

### 62 1e - 14) { 

### 63 string < - "25,000 to < 500,000 lbs" 

### 64 lower < - 25000/2.204/365.25 

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

### 66 } 

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

### 68 { 

### 69 string < - "< 25,000 lbs" 

### 70 lower < - 6250/2.204/365.25 

### 71 upper < - 25000/2.204/365.25 

### 72 } 

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

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

### 75 string < - "10 to < 50 million lbs" 

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

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

### 78 } 

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

### 80 1e+06/2.2 04/365.25))/2)) < 1e - 14) { 

### 81 string < - "50 to < 100 million lbs" 

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

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

### 84 } 

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

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

### 87 string < - "100 to < 500 million lbs" 

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

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

### 90 } 

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

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

### 93 string < - "500 million to < 1 billion lbs" 

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

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

### 96 } 

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

### 98 string < - "1 billion lbs and greater" 

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

### 100 upper < - Inf 

### 101 } 

### 102 else { 

### 103 string < - x 

### 104 lower < - NA 

### 105 upper < - NA 

### 106 } 

### 107 out < - data.frame(string) 

### 108 out < - cbind(out, data.frame(lower)) 

### 109 out < - cbind(out, data.frame(upper)) 

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

### 111 return(out) 

### 112 } 

### 113 114 ## ===================================================================== 

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

### 117 118 NHANES < - subset(NHANES, !Inferred) 

### 119 120 get2005Table < - function(name){ 

### 121 ## escape special characters in name 

### 122 name < - gsub(" \ \ ["," \ \ \ \ [",gsub("]"," \ \ \ \ ]",name)) 

### 123name < - gsub(" \ \ ("," \ \ \ \ (",gsub(")"," \ \ \ \ )",name)) 

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

### 125 *.* \ \ (creatinine",sep=""),tmp) 

### 126 if (length(indx) > 0) { 

### 127 indx2 < - grep("SurveySampleyearssizeTo tal, age 6 and 

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

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

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

### 131 ## need to pull out a unit indicator 

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

### 133 indx2[1]]) 

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

### 135 dt < - gsub("< LOD","NA(NA - NA)",dt) 

### 136 dt < - gsub(" \ \ *","NA(NA - NA)",dt) 

### 137 ## split the string at ## - ## 

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

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

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

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

### 142 ## Split the rows of dt into fields. 

### 143 year < - character(length(dt)) 

### 144 for (i in seq(along=year)) { 

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

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

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

### 148 } 

### 149 mytb < - do.call("rbind",str split(dt,")| \ \ (| - ")) 

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

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

### 152 data.frame(sample=year, 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

### 168 } 

### 169 170 } 

### 171 172 173 174 ## --------------------------------------- ----------------------------- -------------- 

### 175 ## Data from 2011 report 

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

### 177 doc < - xmlParse(tmp2) 

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

### 179 tables < - getNodeSet(doc,"//Table") 

### 180 ## The first row contains the chemical labels 

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

### 182 creaIndex < - grep(" \ \ (creatinine",labels) 

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

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

### 187 v < - gsub(" ", "",v) 

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

### 189 n < - strsplit(v," - ") 

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

### 191 } 

### 192 193 194 195 get2011Table < - function(name){ 

### 196 name < - gsub(" \ \ ["," \ \ \ \ [",gsub("]"," \ \ \ \ ]",name)) 

### 197 name < - gsub(" \ \ ("," \ \ \ \ (",gsub(")"," \ \ \ \ )",name)) 

### 198 name < - gsub(" "," *",name) 

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

### 200 if (length(indx) > 0) { 

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

### 202 tb < - tables[creaIndex][[indx]] 

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

### 204 mytb < - do.call("rbind", 

### 205 lapply(Rows, 

### 206 function(z) { 

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

### 208 } 

### 209 ) 

### 210 ) 

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

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

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

### 214 ## first field contains "Age" 

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

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

### 217 ## Now, drop the spacer rows 

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

### 219 mytb < - mytb[ - indx,] 

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

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

### 222 mydf < - 

### 223 data.frame( 

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

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

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

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

### 228 Value 50 = as.numeric(mytb[,4]), 

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

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

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

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

### 233 UCL50 = getcl(my tb[,5],lim=2), 

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

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

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

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

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

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

### 240 ) 

### 241 mydf 

### 242 } 

### 243 } 

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

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

### 247 dta < - list() 

### 248 249for (i in 1:nrow(NHANES)) { 

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

### 251 get2005Table(NHANES$Name[i]) 

### 252 } else { 

### 253 get2011Table(NHANES$Name[i]) 

### 254 } 

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

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

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

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

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

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

### 261 ), 

### 262 dti, stringsAsFactors=FALSE) 

### 263 } 

### 264 NHANES.old < - NHANES 

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

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

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

### 268 ## table labels and data. 

### 269 write.csv(rbind(NHANES[NHANES$Year == 

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

### 271 NHANES[NHANES$Year == 

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

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

### 274 275 ## 

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

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

### 278 08"),ordered=TRUE) 

### 279 280 281 parent < - ExpoCast.NHANES.compare 

### 282 ## parent inadvertantly includes arsenic. Drop that 

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

### 284 parent < - parent[ - indx,] 

### 285 286 rownames(parent) < - as.character(parent$CASRN) 

### 287 288 mapping < - NHANES.parent.metabolite.mapping 

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

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

### 291 colnames(mapping)) 

### 292 mapping < - mapping[, - indx] 

### 293 294 ## This gives three data sets: 

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

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

### 297 NHANES 

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

### 299 ## chemicals measured in NHANES. 

### 300 ## We want to eliminate every parent that either does no t have a product in NHANES or 

### 301 is not in 

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

### 303 either directly 

### 304 ## or through mapping. 

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

### 306 select 

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

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

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

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

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

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

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

### 315 316 toadd < - intersect(prodcas, parcas) 

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

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

### 319 df < - data.frame(Name=tmp$Compound, C AS=tmp$CASRN,MW=rep(NA, nrow(tmp)), 

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

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

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

### 323 MW < - tmp$MW 

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

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

### 326 327 mapping < - rbind(mapping, df) 

### 328 329 ## 2) eliminate rows in mapping for products not in prodcas 

### 330 ## recreate mapprodcas and mapparcas 

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

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

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

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

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

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

### 338 ## Drop those rows: 

### 339 mapping < - mapping[ - indx,] 

### 340 ## rebuild mapprodcas and mapparcas 

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

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

### 343 344 ## 2) What are the chemicals in prodc as that are not in mapprodcas? 

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

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

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

### 348 ## Name CAS MW LOD Year Value sample 

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

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

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

### 352 ## LCL90 UCL90 LCL95 UCL95 

### 353 ## 70 0.507 0.697 0.736 1.12 

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

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

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

### 357 ## Now rebuild mapprodcas and mapparcas again 

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

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

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

### 361 ## length(notinmapprod) == 0 

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

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

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

### 366 setdiff(mapparcas, parcas) 

### 367 ## > setdiff(mapparcas, parcas) 

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

### 369 ## These are: 

### 370 mapping[indx < - which(as.character(mapping$CAS) %in% setdiff(mapparc as, parcas)),] 

### 371 ## Name CAS MW Name.1 CAS.1 MW.1 

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

### 373 ## 94 Isazaphos - methyl 42509 - 83 - 1 313.74 Dimethylphosphate 813 - 78 - 5 126 .050 

### 374 ## 109 Isazaphos - methyl 42509 - 83 - 1 313.74 Dimethylthiophosphate 1112 - 38 - 5 142.110 

### 375## Drop them from mapping 

### 376 mapping < - mapping[ - indx,] 

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

### 378 ## the set mapping$CAS should equal the set paren t$CASRN 

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

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

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

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

### 383 ## Compare parents in the two structures 

### 384 setdiff(parcas, mapparcas) 

### 385 ## > setdiff(parcas, mapparcas) 

### 386 ## character(0) 

### 387 setdiff(mapparcas, parcas) 

### 388 ## > setdiff(mapparcas, parcas) 

### 389 ## character(0) 

### 390 setdiff(mapprodcas, prodcas) 

### 391 ## > setdiff(mapprodcas, prodcas) 

### 392 ## character(0) 

### 393 setdiff(prodcas, mapprodcas) 

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

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

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

### 397 NHANES 

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

### 399 NHANES < - NHANES[ - indx,] 

### 400 401 ## Add PUsource to mapping 

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

### 403 ## ============================================================= ===================== 

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

### 405 using 01 - 02, 

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

### 407 NHANES chemicals 

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

### 409 ## For each product CAS number (mapping$CAS.1) write down which year range to r etain. 

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

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

### 412 02")[x])) 

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

### 416 tst < - NHANES 

### 417 tst$keep < - keep[tst$CAS] 

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

### 419 420 ## ------------------------------------------------------- ------------------------- 

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

### 423 ## keep the most recent earlier record. 

### 424 425 tst2 < - do.call("rbind", 

### 426 unclass(by(tst,tst$CAS, 

### 427 function(x) { 

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

### 429 if (length(indx ) > 0) { 

### 430 ## exact match. Easy 

### 431 z < - x[indx,,drop=FALSE] 

### 432 z$match < - "==" 

### 433 z 

### 434 435 } else { 

### 436 if (length(x$sample) == 1) { 

### 437 ## only 1 row, so just use it. 

### 438x$match < - "1" 

### 439 x 

### 440 } else { 

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

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

### 443 levels=levels(tst$sample), ordered=TRUE) 

### 444 kp < - unique(x$keep) 

### 445 if (any(z1 < kp)) { 

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

### 447 z < - x[indx,,drop=FALSE] 

### 448 z$match < - "<" 

### 449 z 

### 450 } else { 

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

### 452 z$match < - ">" 

### 453 z 

### 454 } 

### 455 } 

### 456 } 

### 457 } 

### 458 ))) 

### 459 Measured < - tst2 

### 460 461 ## ----------------------------------------------------- --------------------- 

### 462 ## Create pred.data 

### 463 464 pred.data < - 

### 465 data.frame(CAS=parent$CASRN, 

### 466 Name=parent$Compound, 

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

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

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

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

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

### 472 473 V.cat < - character(nrow(pred.data)) 

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

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

### 476 V.cat[i] < - NA 

### 477 } else { 

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

### 479 } 

### 480 } 

### 481 levels < - c("< 25,000 lbs", 

### 482 "25,000 to < 500,000 lbs", 

### 483 "500,000 to < 1 million lbs", 

### 484 "1 to < 10 million lbs", 

### 485 "10 to < 50 million lbs", 

### 486 "50 to < 100 million lbs", 

### 487 "100 to < 500 million lbs", 

### 488 "500 million to < 1 billion lbs", 

### 489 "1 billion lbs and greater") 

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

### 491 492 ## Create mapping 

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

### 495 ## save it 

### 496 oldmap < - mapping 

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

### 499 mapping < - subset(mapping, 

### 500 select=c("Name","CAS","Name.1","CAS.1")) 

### 501

### 502 ## ----------------------------------------------------------------------- 

### 503 ## Convert measurements 

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

### 506 507 NHANES$gm < - convert.urine. to.exposure(NHANES$gm) 

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

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

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

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

### 512 N HANES$Value90 < - convert.urine.to.exposure(NHANES$Value90) 

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

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

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

### 516 NHANES$LCL90 < - convert.u rine.to.exposure(NHANES$LCL90) 

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

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

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

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

### 521 N HANES$UCL95 < - convert.urine.to.exposure(NHANES$UCL95) 

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

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

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

### 526 ## and match to th e value of 'sample' in Measured. This value 

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

### 528 ## by average g creatinine / L urine 

### 529 530 LOQ2005 < - read.table("LOQ2005.txt",header=TRUE, 

### 531 532 colClasses=c("character","character","character","numeric","numeric"), 

### 533 sep="|") 

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

### 535 colClasses=c("character","character","character", 

### 536 "numeric","numeric","numeric","numeric","numeric"), 

### 537 sep="|") 

### 538 LOQ < - numeric(nrow(Measured)) 

### 539 for (i in 1:nrow(Measured)) { 

### 540 name < - Measured$Name[i] 

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

### 542 name < - gsub(" \ \ ["," \ \ \ \ [",gsub("] "," \ \ \ \ ]",name)) 

### 543 name < - gsub(" \ \ ("," \ \ \ \ (",gsub(")"," \ \ \ \ )",name)) 

### 544 name < - paste("^",name,sep="") 

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

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

### 547 LOQ2005",sep="")) 

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

### 549 0.001) 

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

### 551 ## LOQ in ug / g creatinine 

### 552 ## In 99 - 00, that was 1.364 g/L (creatinine / uri n e) (page 5 in NHANES 2005 

### 553 ## document 

### 554 ## in 01 - 02 that was 1.306 

### 555 LOQ[i] < - 

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

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

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

### 559 } else { 

### 560 name < - gsub(" \ \ ["," \ \ \ \ [",gsub("]"," \ \ \ \ ]",name)) 

### 561 name < - gsub(" \ \ ("," \ \ \ \ (",gsub(")"," \ \ \ \ )",name)) 

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

### 563 name < - paste("^",name,sep="") 

### 564indx < - grep(name, LOQ2011$Name, ignore.case=TRUE) 

### 565 if (len gth(indx) == 0) 

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

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

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

### 569 ## value from NHANES III 

### 570 LOQ[i] < - 

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

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

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

### 574 "05 - 06"=LOQ2011[indx,"Y 2005.2006"], 

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

### 576 } 

### 577 } 

### 578 ## LOQ is in ug/g creatinine. 

### 579 580 Measured$LOQ < - convert.urine.to.exposure(LOQ) 

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

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

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

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

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

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

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

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

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

### 590 Measured$LCL90 < - convert.ur ine.to.exposure(Measured$LCL90) 

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

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

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

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

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

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

### 597 598 ## --------------------------------------------------------------- 

### 599 ## Build a matrix with rows for parents and columns for metabolite s 

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

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

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

### 603 ## 1 across rows), 

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

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

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

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

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

### 609 ## listed below, 

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

### 611 ## are also in NHANES. 

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

### 613 dimnames=l ist(pred.data$CAS, Measured$CAS)) 

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

### 615 match(mapping$CAS.1,colnames(Xref))) 

### 616 617 Xref[indx] < - NA 

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

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

### 620 ## This is probably crude 

### 621 for (i in seq(along=Nmets)) { 

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

### 623 } 

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

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

### 626 ## A few parent - metabolite pairs are known to be 1 - 1. 

### 627fixed < - cbind(Parent=c("2921 - 88 - 2","68359 - 37 - 5","52918 - 63 - 5","52918 - 63 - 5", 

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

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

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

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

### 632 ## to 1 

### 633 indx < - cbind(match(fixed[,1],rownames(Xref)),match(fixed[,2], c olnames(Xref))) 

### 634 Xref[indx] < - 1 

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

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

### 638 ## parent? 

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

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

### 641 N.metabs < - N.metabs[indx] 

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

### 643 ## assign useful rownames to pred.data: 

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

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

### 646 pred.data < - tmp 

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

### 648 649 ## +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

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

### 651 ## and metabolite (product) concentrations to moles 

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

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

### 654 ## rest. 

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

### 657 Parent MW < - unique(subset(oldmap, select=c("CAS","MW"))) 

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

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

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

### 662 nam es(MetaboliteMW) < - c("CAS","MW") 

### 663 664 ## Combine these two. What are missing? 

### 665 666 MW < - unique(rbind(ParentMW, MetaboliteMW)) 

### 667 rownames(MW) < - as.character(MW$CAS) 

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

### 670 rownames(MW)))]) 

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

### 673 rownames(MW) < - as.character(MW$CAS) 

### 674 675 676 677 ## Need to convert s.U, and s.R to use molar units. 

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

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

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

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

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

### 684 685 Measured$Value < - with(Measured, Value /(MW * 1000)) 

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

### 687 688 Measured$LOQ < - with(Measured, LOQ / (MW * 1000)) 

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

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

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

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

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

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

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

### 696 Measured$LCL50 < - with(Measured, LCL50 /(MW * 10 00)) 

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

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

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

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

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

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

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

### 704 705 ## Measurements are moles / kg bodyweight / day 

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

### 707 ## nmoles / kg bodyweight / day 

### 708 p red.data[,c("s.U","s.R")] < - 1e9 * pred.data[,c("s.U","s.R")] 

### 709 710 Measured[,c("Value","LOD","LOQ","gm","gmlcl","gmucl", 

### 711 "Value50","Value75","Value90","Value95", 

### 712 "LCL50","LCL75","LCL90","LCL95", 

### 713 "UCL50","UCL75","UCL90","UCL9 5")] < - 

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

### 715 "Value50","Value75","Value90","Value95", 

### 716 "LCL50","LCL75","LCL90","LCL95", 

### 717 "UCL50","UCL75","UCL90","UCL95")] 

### 718 ## Now, values in 'Measured' a re nmoles / kg bodyweight / day 

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

### 721 ## Convert cutpoints to kg / day 

### 722 cutpoints < - cutpoints * 0.45359237 / 365.25 

### 723 724 ## Standard errors for Measured. If all values a re above LOD, then sexx is 

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

### 726 727 ## Median 

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

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

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

### 731 732 se2 < - with(Measured, ((log(UCL75) - log(LCL75)) / (2 * 1.96))) 

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

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

### 735 736 se2 < - with(Measured, ((log(UCL90) - log(LCL90)) / (2 * 1.9 6))) 

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

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

### 739 740 se2 < - with(Measured, ((log(UCL95) - log(LCL95)) / (2 * 1.96))) 

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

### 742 Measured$se95 < - ifelse(i s.na(se2), se1,se2) 

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

### 745 ## Build a data frame. 

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

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

### 748 timevar="Quantile") 

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

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

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

### 752 ## Now, get the predictions for all the values in Zd 

### 753anova(out,o ut2) 

### 754 anova(out,out3) 

### 755 756 ## No real differences among quantiles, but significant differences among chems. 

### 757 ## So, set a single se for each chemical, the sqrt of the average of the quant 

### 758 ## specific variances. 

### 759 760 se < - apply(data.matrix(Z),1,function(x) { 

### 761 if (!all(is.na(x))) { 

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

### 763 } else NA 

### 764 }) 

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

### 767 768 se[is.na(se)] < - mean(se, na.rm=TRUE) 

### 769 770 Measured$se < - se 

### 771 attach("UseData - 041013.RData") 

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

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

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

### 775 Uses < - AliciaUseTable[indx,] 

### 776 ##Uses2 < - UseData[indx2,] 

### 777 Uses$NearField < - as.numeric(with(Uses, FRAGRANCE | F OOD.ADDITIVE | CONSUMER.USE | 

### 778 PERSONAL.CARE.PRODUCT | PHARMACEUTICAL)) 

### 779 pred.data$NearField < - Uses$NearField 

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

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

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

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

### 784 file=file.path("..","MCMCdata.RData")) 

### 785 786 ## Counts for uses in the NHANES data: 

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

### 789 ncol=1, 

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

### 791 792 writeLines(c("","Breakdown for Near Field:", 

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

### 794 "Order is: Consumer use, food addi tive, fragrance, personal care product, 

### 795 pharmaceutical", 

### 796 "so 10100 is classified consumer use and fragrance")) 

### 797 tmp < - table(with(Uses, 

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

### 800 ))) 

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

### 802 writeLines(c("", 

### 803 "Breakdown for the whole dataset:", 

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

### 805 product,", 

### 806 "antimicrobial,chemical industrial proc ess, colorant,", 

### 807 "fertilizer, herbicide, pesticide, pharmaceutical, toxin,inert, flame 

### 808 retardant")) 

### 809 tmp < - table(with(Uses, 

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

### 811 ANTIMICROBIAL,CHEMICAL.IN DUSTRIAL.PROCESS,COLORANT, 

### 812 FERTILIZER,HERBICIDE, 

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

### 814 815 matrix(tmp, ncol=1, dimnames=list(names(tmp), "Count")) 

### 816

# OnlyP/ model0.jag: JAGS Code to Estimate Parent Chemical Exposures

from NHANES model 

### 1 { 

### 2 ## +++++++++++++++++++++++++++++++++++++++++++++ 

### 3 ## The heart of the model: 

### 4 ## Model the parent exposures 

### 5 ## +++++++++++++++++++++++++++++++++++++++++++++ 

### 6 for (i in 1:N) { 

### 7 ## log of product of pro duction volume and unit exposure 

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

### 9 P[i]< - exp(lP[i]) 

### 10 } 

### 11 ## hyperparameter for lP 

### 12 lPmu ~ dnorm(0, 0.001) 

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

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

### 16 sd.V < - abs(sd.dum) 

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

### 18 ## ++++++++++++++++++++++++++++++++++++++++++++++ 

### 19 ## Link the unobserved parent exposure to the observed 

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

### 21 ## values. 

### 22 ## 

### 23 ## Model metabolite concentrations as (usually) unkn own fractions 

### 24 ## of the parent. 

### 25 ## 

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

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

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

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

### 30 for (j in 1:M) { 

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

### 32 for (i in 1:4) { 

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

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

### 35 } 

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

### 37 ## lognormal distribution for chemical - specific population sds 

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

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

### 40 } 

### 41 lsd ~ dnorm(0, 0.001) 

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

### 43 lsd.sd < - abs(lsd.dum) 

### 44 tau.lsd < - 1/(l sd.sd * lsd.sd) 

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

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

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

### 48 ## indx[i,1:Nmets[i]] give t he columns of row i that are being estimated, 

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

### 50 ## metabolite matched to parent i. 

### 51 ## 

### 52 ## The elements of phi are essentially unknown. 

### 53 for (k in 1:Ndelta) { 

### 54 delta[k] ~ dunif(0,1) 

### 55 } 

### 56 for (i in 1:N.nontriv) { 

### 57 for (j in 1:Nmets[i]) { 

### 58phi[dstart[i] + j - 1] < - delta[dstart[i] + j - 1] / 

### 59 sum(delta[dstart[i]:dstop[i]]) 

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

### 61 } 

### 62 } 

### 63 } 

### 64

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

NHANES library(doMC) 

### 1 registerDoMC(cores=3) 

### 2 library(rjags) 

### 3 load.module("lecuyer") 

### 4 library(random) 

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

### 6 7 doONEchain < - TRUE 

### 8 9 ## JAGS model is in model0. jag 

### 10 ## It needs the following data: 

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

### 12 ## N = number of parents, 

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

### 14 ## in R.x) 

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

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

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

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

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

### 20 ## M = number of metabolites, 

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

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

### 23 ## metabolite r ates. 

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

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

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

### 27 ## 

### 28 ## The following variables need initiali zing 

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

### 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 if (doONEchain) { 

### 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 quant iles 

### 63 nhanesdata$ly < - 

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

### 65 ## Q, the proportions associated with the quantiles 

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

### 67 ## se - precision of measured metabolite value s, from reported confidence limits on 

### 68 quantiles 

### 69 nhanesdata$se < - Measured$se 

### 70 ## llod - log(Limit of Detection) 

### 71 nhanesdata$llod < - log(Measured$LOQ) 

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

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

### 74 ## Some counts: 

### 75 ## N is the number of parents 

### 76 nhanesdata$N < - nrow(pred.data) 

### 77 ## M is the number of metabolites 

### 78 nhanesdata$M < - nrow(Measured) 

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

### 80 nhan esdata$N.IUR < - length(nhanesdata$IURindx) 

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

### 82 nhanesdata$N.CPRI < - length(nhanesdata$CPRIindx) 

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

### 84 nhanesdata$N.no ntriv < - sum(N.metabs > 1) 

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

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

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

### 88 ## of the ith parent 

### 89 nhanesd ata$indx < - 

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

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

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

### 93 tmp < - cumsum(nhanesdata$Nmets) 

### 94 nhanesdata$dstart < - 1+c(0,tmp[ - nha nesdata$N.nontriv]) 

### 95 nhanesdata$dstop < - tmp 

### 96 97 ## Phi 

### 98 nhanesdata$Phi < - Xref 

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

### 100 101 ## --------------------------------------------------------------------- 

### 102 ## Initialization 

### 103 104 inits < - list(list()) 

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

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

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

### 108 for (j in 1:ncol(ly)) { 

### 109 ## if there are at least 2 values > LOQ, then est imate lmu and sd using least 

### 110 squares 

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

### 112 z < - ly[,j] 

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

### 114 keep < - !is.na(z) 

### 115 x < - Q[keep] 

### 116 y < - z[keep] 

### 117 pstart < - c(min(y) - 1,0) 

### 118 obfun < - function(p) { 

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

### 120 } 

### 121out < - try(optim(pstart, obfun)) 

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

### 123 lU[j] < - out$par[1] 

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

### 125 } else { 

### 126 lU[j] < - min(y) 

### 127 inits[[1]]$sd[j] < - 1 

### 128 } 

### 129 } else { 

### 130 lU[j] < - nhanesdata$llod[j] - 1 

### 131 inits[[1]]$sd[j] < - 1 

### 132 } 

### 133 } 

### 134 XX < - Xref 

### 135 for (i in 1:nrow(XX)) { 

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

### 137 if (nNA > 0) { 

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

### 139 } 

### 140 } 

### 141 XXt < - t(XX) 

### 142 meanlU < - mean(lU) 

### 143 Usc < - exp(lU - meanlU) 

### 144 obfun < - function(p) { 

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

### 146 } 

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

### 148 pstart[pstart< 0] < - .1 

### 149 pstart < - log(pstart) 

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

### 151 152 inits[[1]]$lP < - out$par + meanlU 

### 153 154 inits[[1]]$sd.dum < - sd(inits[[1]]$lP) 

### 155 156 ## for the censored values in ly 

### 157 ily < - is.na(ly) 

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

### 159 lyinit < - ifelse(ily,loqmx - 2,NA) 

### 160 inits[[1]]$ly < - lyinit 

### 161 162 inits[[1]]$lsd < - mean(log(inits[[1]]$sd)) 

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

### 164 165 inits[[1]]$delta < - rep(1, nhanesdata$Ndelta) 

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

### 167 168 doinits < - function(chain) {inits[[chain]]} 

### 169 assign("inits", inits, envir=environment(doinits)) 

### 170 171 ## ------------------------------------------------- 

### 172 ## Run JAGS 

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

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

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

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

### 178 N.Iters < - 1000 

### 179 N.Thin < - 500 

### 180 N.Burnin < - 1000 

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

### 182 n.chains=1,n.adapt=N.Burnin) 

### 183 184out.samples < - coda.samples(model.jags, 

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

### 186 "sd","lsd","lsd.sd"), 

### 187 n.iter=N.Iters * N.Thin, thin=N.Thin) 

### 188 189 nms < - varnames(out.samples) 

### 190 for (nm in nms) { 

### 191 writeLines(c(""," ----------------- ",nm,"")) 

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

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

### 194 } 

### 195 196 save(out.samples, nhanesdata, inits, file="m0parms1.RData") 

### 197 } 

### 198 ## -------------------------------------------------------------------------------- 

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

### 200 ## and run three new chains. 

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

### 203 ## chains are hyperdispersed relative to the p osterior distribution. getInits 

### 204 ## accomplishes that by first picking ninits random dire ctions in parameter space, 

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

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

### 207 ## return value is the row index of the resulting sampl e in the input X. 

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

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

### 210 dta < - scale(X) 

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

### 212 outinits < - numeric(ninits) 

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

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

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

### 216 for (i in 1:ninits) { 

### 217 u < - rnorm(ncol(dta)) 

### 218 u < - u/sqrt(su m(u*u)) 

### 219 lngt < - dta %*% u 

### 220 indx < - seq(along=lngt) 

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

### 222 oind < - order(lngt) 

### 223 outinits[i] < - oind[keep] 

### 224 } 

### 225 outinits 

### 226 } 

### 227 ## Get 3 hyperdispersed points. 

### 228 out.coda < - o ut.samples[[1]] 

### 229 230 Indx < - getInits(out.coda,3) 

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

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

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

### 235 getsamps < - funct ion(nm) { 

### 236 nm < - gsub(" \ \ ."," \ \ \ \ .",nm) 

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

### 238 out.coda[,cindx,drop=FALSE] 

### 239 } 

### 240 241 ## We run three chains. 

### 242 inits2 < - list(list(),list(),list()) 

### 243 for (i in 1:3) { 

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

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

### 246 inits2[[i]]$sd < - getsamps("^sd \ \ [")[Indx[i],] 

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

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

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

### 250 i nits2[[i]]$ly< - inits[[1]]$ly 

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

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

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

### 254 } 

### 255 256 doinits < - function(chain){inits2[[chain]]} 

### 257 env1 < - new.env() 

### 258 assign("inits2", inits2, pos=env1) 

### 259 environment(doinits) < - env1 

### 260 261 NIters < - 1000 

### 262 Thin < - 500 

### 263 NBurnin < - 5000 

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

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

### 267 n.a dapt=NBurnin) 

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

### 269 "sd","lsd","lsd.sd"), 

### 270 n.iter=NIters*Thin, thin=Thin) 

### 271 return(result) 

### 272 } 

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

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

### 276 nms < - varnames(out.coda3R) 

### 277 for (nm in nms) { 

### 278 writeLines(c(""," ----------------- ",nm,"")) 

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

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

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

### 282 } 

### 283

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

Volume library(rjags) 

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

### 2 ## Compare RAIDAR and U SEtox to NHANES urinary data 

### 3 ## JAGS model is in nhanes2.jag 

### 4 ## It needs the following data: 

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

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

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

### 8 ## N = numbe r of parents, 

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

### 10 ## in R.x) 

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

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

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

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

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

### 16 ## M = number of metabolites, 

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

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

### 19 ## metabolite rates. 

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

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

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

### 23 ## 

### 24 ## The following variables need initializing 

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

### 26 ## a 

### 27 ## m.R 

### 28 ## m.U 

### 29 ## X1 

### 30 ## X2 

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

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

### 34 35 attach(file.path("..","MCMCdata.RData")) 

### 36 37 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

### 38 ## JAGS inputs: nhanesdata and inits 

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

### 40 41 nhanesdata < - list() 

### 42 ## The Jags models need these data: 

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

### 44 nhanesdata$IURcuts < - log(cutpoints) 

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

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

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

### 48 nhanesdata$IURvols < - 

### 49 as.numeric(pred.data$V.ca t[nhanesdata$IURindx[1:length(nhanesdata$IURindx)]]) 

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

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

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

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

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

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

### 56 nhanesdata$NearField < - pred.data$NearField 

### 57 ## Mean of ls.R and ls.U, for centering 

### 58nhanesdata$meanls.R < - mean(nh anesdata$ls.R) 

### 59 nhanesdata$meanls.U < - mean(nhanesdata$ls.U) 

### 60 61 ## ly - the log of the observed metabolite quantiles 

### 62 nhanesdata$ly < - 

### 63 ly < - 

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

### 65 ## Q, the proportions associated wit h the quantiles 

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

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

### 68 quantiles 

### 69 nhanesdata$se < - Measured$se 

### 70 ## llod - log(Limit of Detection) 

### 71 nhanesdata$llod < - log(Measured$LOQ) 

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

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

### 74 ## Some counts: 

### 75 ## N is the number of parents 

### 76 nhanesdata$N < - nrow(pred.data) 

### 77 ## M is the number of metabolites 

### 78 nhanesdata$M < - nrow(Measured) 

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

### 80 nhanesdata$N.IUR < - length(nhanesdata$IURindx) 

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

### 82 nhanesdata$N.CPRI < - length(nhanesdata$CPRIindx) 

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

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

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

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

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

### 88 ## of the ith parent 

### 89 nhanesdata$indx < - 

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

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

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

### 93 tmp < - cumsum(nhanesdata$ Nmets) 

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

### 95 nhanesdata$dstop < - tmp 

### 96 97 ## Phi 

### 98 nhanesdata$Phi < - Xref 

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

### 100 101 ## --------------------------------------------------------------------- 

### 102 ## Initialization 

### 103 104 inits < - list() 

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

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

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

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

### 109 ## inits$lP: get estimates for lU, then work backwards to l P 

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

### 111 ## JAGS code 

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

### 114 115 116 ## Now, run for 50000 samples, thinned every 100, and dump lP, lU, and delta 

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

### 118 119 out.samples < - jags.samples(model, 

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

### 121n.iter=100000, thin=200) 

### 122 out.coda < - lapply(out.samp les, as.mcmc.list) 

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

### 124 lapply(out.coda, function(z) summary(z)) 

### 125 126 save(out.samples, out.coda, nhanesdata, inits, file="m0parms1.RData") 

### 127 ## -------------------------------------------------------------------------------- 

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

### 129 ## and run three new chains. 

### 130 131 132 ## getInits finds extreme values 

### 133 getInits < - function(X, nini ts=3) { 

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

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

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

### 137 rows), 

### 138 ## and center and scale the result. 

### 139 dta < - scale(t(do.call("rb ind",tmp1))) 

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

### 141 outinits < - numeric(ninits) 

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

### 143 ## and pick the 5th percentile (since the original directi on u was random u and - u 

### 144 are equally likely). 

### 145 for (i in 1:ninits) { 

### 146 u < - runif(ncol(dta), - 1,1) 

### 147 u < - u/sqrt(sum(u*u)) 

### 148 lngt < - dta %*% u 

### 149 indx < - seq(along=lngt) 

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

### 151 oind < - order(lngt) 

### 152 outinits[i] < - oind[keep] 

### 153 } 

### 154 outinits 

### 155 } 

### 156 ## Get 3 hyperdispersed points. 

### 157 Indx < - getInits(out.samples,3) 

### 158 159 inits2 < - list(list(),list(),list()) 

### 160 for (i in 1:3) { 

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

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

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

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

### 165 } 

### 166 167 168 model < - jags.model("model0.jag",data=nhanesdata,inits=inits2,n.chains=3, 

### 169 n.adapt=45000) 

### 170 171 adaptdone < - FALSE 

### 172 Count < - 0 

### 173 while (!adaptdone && (Count < 10)) { 

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

### 175 Count < - Count + 1 

### 176 } 

### 177 178 out.samps3R < - jags.samples(model, 

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

### 180 n.iter=500000, thin=500) 

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

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

### 183 184for (z in names(out.coda3R)) { 

### 185 writeLines(paste("Variable: ",z)) 

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

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

### 188 print(summary(out.coda3R[[z]])) 

### 189 writeLines(c(" "," ------------------------- ------- "," ")) 

### 190 } 

### 191 192 193 194 195

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

Volume model 

### 1 { 

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

### 3 ## sources, with their own error characteristics. The CPRI dat a are 

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

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

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

### 7 ## scale, that means lot s of normal distributions). IUR data are then 

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

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

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

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

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

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

### 14 ## predictions 

### 15 ## IUR volumes 

### 16 for (i in 1:N.IUR) { 

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

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

### 19 lV[IURindx[i]] < - lV.IUR[i] 

### 20 } 

### 21 lVmu ~ dnorm(0, 0.001) 

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

### 23 sd.lVmu < - abs(sdlVmu.dum) 

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

### 25 ## -------------- ----------------------------- 

### 26 ## CPRI usage 

### 27 for (i in 1:N.CPRI) { 

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

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

### 30 lV[CPRIindx[i]] < - lV.CPRI.mu[i] 

### 31 } 

### 32 } 

### 33

# analysisOmnibus

- F its.R: Evaluating Alternative Models for Predicting

NHANES Chemical Exposures library(compiler) 

### 1 setCompilerOptions(optimize=3) 

### 2 enableJIT(3) 

### 3 library(nlme) 

### 4 library(mvtnorm) 

### 5 library(multicore) 

### 6 source("QFunctions.R") 

### 7 ## This takes a two - stage approach to the problem of validation and calibration 

### 8 ## of RAIDAR and USEtox against NHANES. In the first stage, samples from the 

### 9 ## posterior distribution of parent exposures are drawn using MCMC methods, 

### 10 ## using code in the 'OnlyP' directory, AND samples from the pro duction volume 

### 11 ## posteriors are drawn using code in the 'OnlyV' directory. See those 

### 12 ## models for the details. 

### 13 ## The second stage is to do a conventional analysis, looking for 

### 14 ## relationships between lP and ls.R, ls.U, and lV (and combinations thereof ), 

### 15 ## AND to look for influential observations. This analysis is done for each 

### 16 ## sample in lV and lP computed in the first stage. We use the distributions 

### 17 ## of the various statistics over the samples of lP and lV to make our 

### 18 ## inferences. 

### 19 20 ## FittingM odels controls whether we actually fit models in 

### 21 ## this run, or just print output (mainly for fixing 

### 22 ## bugs in the output code) 

### 23 24 FittingModels < - TRUE 

### 25 26 ## Attach the initial data: 

### 27 28 attach(file.path("..","MCMCdata.RData")) 

### 29 30 ## Get the lP samples 

### 31 32 attach(file.path("..","OnlyP","m3parms3.RData")) 

### 33 34 tmp < - do.call("rbind",out.coda3R) 

### 35 indx < - grep("^lP \ \ [",colnames(tmp)) 

### 36 lP.mx < - tmp[,indx] 

### 37 38 detach(2) 

### 39 40 attach(file.path("..","OnlyV", "m0parms3.RData")) 

### 41 42 lV.mx < - do.call("rbind",out.coda3R$lV) 

### 43 ## We use h alf the sample to estimate correlation and variances, 

### 44 ## the other half to do the testing. 

### 45 ## Get a random sample of half the rows for testing: 

### 46 47 if (FittingModels) { 

### 48 49 CorSampIndx < - sample(3000,1500) 

### 50 RegSampIndx < - (1:3000)[ - CorSampIndx] 

### 51 52 ## Now comp ute the covariance matrix for lV 

### 53 lP.cov < - cov(lP.mx[CorSampIndx,]) 

### 54 lP.cor < - cov2cor(lP.cov) 

### 55 56 corparms < - lP.cor[lower.tri(lP.cor)] 

### 57 lP.vars < - diag(lP.cov) 

### 58

### 59 60 ## This will be easier if everything is in a dataframe. Everything 

### 61 ## but lP and lV will be a constant, so just build it, and leave space 

### 62 ## for those two variables, which will be assigned at each iteration. 

### 63 ## Use the first row of each variable to initialize it. 

### 64 65 66 ## Build a function that carries out the analysis we want to do for 

### 67 ## each pair of lP and lV in lP.mx and lV.mx. The analysis is 

### 68 ## designed to estimate parameters for each of several models, 

### 69 ## corresponding to models fit in the full Bayesian analysis, and 

### 70 ## compute statistics designed to hel p evaluate the importance of 

### 71 ## different factors, as well as the best predictive model for lP. We 

### 72 ## do this for two reasons: the full Bayesian model allows lP and lV 

### 73 ## to move around, because lV is censored or measured with error, and 

### 74 ## lP is really derived from the measurements of metabolites 

### 75 ## compounded with unknown stoichiometry, which gives a slightly 

### 76 ## mysterious feel to evaluating fits; computing, e.g., DIC for this 

### 77 ## model seems to be complicated, and is certainly not immediate ly 

### 78 ## obvious. 

### 79 80 ## It seems natural to use gls() to do the fitting, but I cannot 

### 81 ## figure out how to get it to handle the covariance structure I need 

### 82 ## for this problem. The model is lP = f(predictions, parameters) + 

### 83 ## epsilon + zeta epsilon is normal(0, sigma^2), and zeta is MVnorm(0, 

### 84 ## defined covariance), where the defined covariance is estimated from 

### 85 ## the MCMC samples. To carry this out, I make a new fitting function 

### 86 ## mvtnormfit(), which takes as arguments: 

### 87 ## y 

### 88 ## x1,...,x3 (not all required) 

### 89 ## M (model number, according to the above scheme) 

### 90 ## lPcov: the defined covariance matrix 

### 91 ## Return value is: 

### 92 ## loglik: log likelihood 

### 93 ## AIC: AIC 

### 94 ## parameters: named appropriately 

### 95 ## R2: defined as 1 - va r(y - pred) / var(y) 

### 96 97 98 NMs < - c("intercept_0", 

### 99 "sigma_0", 

### 100 "AIC_0","LL_0", 

### 101 102 "intercept_1","Field_1", 

### 103 "sigma_1", 

### 104 "AIC_1","R2_1","LL_1", 

### 105 106 "intercept_2","Field_2","m.lV_2 ", 

### 107 "sigma_2", 

### 108 "AIC_2","R2_2","LL_2", 

### 109 110 "intercept_3","Field_3","m.lR_3", 

### 111 "sigma_3", 

### 112 "AIC_3","R2_3","LL_3", 

### 113 114 "intercept_4","Field_4","m.luR_4","m.lV_4", 

### 115 "sigma_4", 

### 116 "AIC_4","R2_4","LL_4", 

### 117 118 "intercept_5","Field_5","m.lU_5", 

### 119 "sigma_5", 

### 120 "AIC_5","R2_5","LL_5", 

### 121

### 122 "intercept_6","Field_6","m.luU_6","m.lV_6", 

### 123 "sigma_6" , 

### 124 "AIC_6","R2_6","LL_6", 

### 125 126 "intercept_7","Field_7","m.luR_7","m.luU_7","m.lV_7", 

### 127 "sigma_7", 

### 128 "AIC_7","R2_7","LL_7", 

### 129 130 "intercept_8","Field_8","m.lR_8","m.lU_8", 

### 131 "sigma_8", 

### 132 "AIC_8","R2_8","LL_8", 

### 133 134 "intercept_9","Field_9","m.luR_9","m.luU_9","m.lV_9", 

### 135 "m.luR.N_9", 

### 136 "sigma_9", 

### 137 "AIC_9", "R2_9","LL_9", 

### 138 139 "intercept_10","Field_10","m.luR_10","m. luU_10","m.lV_10", 

### 140 "m.luU.N_10", 

### 141 "sigma_10", 

### 142 "AIC_10", "R2_10","LL_10", 

### 143 144 "intercept_11","Field_11","m.luR_11","m.luU_11","m.lV_11", 

### 145 "m.lV.N_11", 

### 146 "sigma_11", 

### 147 "AIC_11", "R2_11","LL_11", 

### 148 149 "intercept_12","Field_12","m.luR_12","m.luU_12","m.lV_12", 

### 150 "m.luR.N_12","m.luU.N_12","m.lV.N_12", 

### 151 "sigma_12", 

### 152 "AIC_12", "R2_12","LL_12", 

### 153 154 "intercep t_13","Field_13","m.luR_13","m.luU_13","m.lV_13", 

### 155 "m.luR.N_13","m.luU.N_13", 

### 156 "sigma_13", 

### 157 "AIC_13", "R2_13","LL_13", 

### 158 159 "intercept_14","Field_14","m.luR_14","m.luU_14","m.lV_14", 

### 160 "m.luU.N_14"," m.lV.N_14", 

### 161 "sigma_14", 

### 162 "AIC_14", "R2_14","LL_14", 

### 163 164 "intercept_15","Field_15","m.luR_15","m.luU_15","m.lV_15", 

### 165 "m.luR.N_15","m.lV.N_15", 

### 166 "sigma_15", 

### 167 "AIC_15", "R2_15","LL_15", 

### 168 169 "intercept_16","Field_16", 

### 170 "m.luR.N_16","m.luU.N_16","m.lV.N_16", 

### 171 "sigma_16", 

### 172 "AIC_16", "R2_16","LL_16", 

### 173 174 "intercept_17","Field_17", 

### 175 "m.luR.N_17","m.luU.N_17", 

### 176 "sigma_17", 

### 177 "AIC_17", "R2_17","LL_17", 

### 178 179 "intercept_18","Field_18", 

### 180 "m.luR.N_18","m.lV.N_18", 

### 181 "sigma_18", 

### 182 "AIC_18", "R2_18","LL_18", 

### 183 184"intercept_19","Field_19 ", 

### 185 "m.luU.N_19","m.lV.N_19", 

### 186 "sigma_19", 

### 187 "AIC_19", "R2_19","LL_19", 

### 188 189 "intercept_20","Field_20", 

### 190 "m.luR.N_20", 

### 191 "sigma_20", 

### 192 "AIC_20", "R2_20","LL_20", 

### 193 194 "intercept_21", "m.lR_21", 

### 195 "sigma_21", 

### 196 "AIC_21","R2_21","LL_21", 

### 197 198 "intercept_22", "m.lU_22", 

### 199 "sigma_22", 

### 200 "AIC_22","R2_22","LL_22", 

### 201 202 "dR0< - >1","dR0< - >21","dR0< - >22", 

### 203 "dR1< - >2","dR1< - >3","dR1< - >4","dR1< - >5","dR1< - >6","dR1< - >12", 

### 204 "dR2< - >4","dR2< - >6", 

### 205 "dR3< - >4", 

### 206 "dR4< - >7", 

### 207 "dR5< - >6", 

### 208 "dR6< - >7", 

### 209 "dR7< - >9","dR7< - >10","dR7< - >11","dR7< - >12", 

### 210 "dR8< - >7", 

### 211 "dR9< - >12", 

### 212 "dR10< - >12", 

### 213 "dR11< - >12", 

### 214 "dR16< - >12", 

### 215 "dR18< - >12", 

### 216 "dR20< - >12", 

### 217 218 "dLL0< - >1","dLL0< - >21","dLL0< - >22", 

### 219 "dLL1< - >2","dLL1< - >3","dLL1< - >4","dLL1< - >5","dLL1< - >6","dLL1< - >12", 

### 220 "dLL2< - >4","dLL2< - >6", 

### 221 "dLL3< - >4", 

### 222 "dLL4< - >7", 

### 223 "dLL5< - >6", 

### 224 "dLL6< - >7", 

### 225 "dLL7< - >9","dLL7< - >10","dLL7< - >11","dLL7< - >12", 

### 226 "dLL8< - >7", 

### 227 "dLL9< - >12", 

### 228 "dLL10< - >12", 

### 229 "dLL11< - >12", 

### 230 "dLL16< - >12", 

### 231 "dLL18< - >12", 

### 232 "dLL20< - >12" 

### 233 234 ) 

### 235 236 analysis < - function(df,drop=NULL) { 

### 237 if (!is.null(drop)) { 

### 238 df < - df[ - d rop,] 

### 239 } 

### 240 out1 < - rep(as.numeric(NA),times=length(NMs)) 

### 241 names(out1) < - NMs 

### 242 243 ## Model 0 

### 244 md0 < - mvtnormfit(df, 0, lP.cov) 

### 245 if (!is(md0, "try - error")) { 

### 246 out1["intercept_0"] < - unname(coef(md0)["intercept"]) 

### 247out1["sigma_0"] < - exp(unname(coef(md0)["lsigma"])) 

### 248 out1["AIC_0"] < - md0$AIC 

### 249 out1["LL_0"] < - md0$LogLik 

### 250 } 

### 251 252 ## Model 1: 

### 253 md1 < - mvtnormfit(df,1,lP.cov) 

### 254 if (!is(md1, "try - error")) { 

### 255 out1["intercept_1"] < - unname(coef(md1)["intercept"]) 

### 256 out1["Field_1"] < - unname(coef(md1)["NearField"]) 

### 257 out1["sigma_1"] < - exp(unname(coef(md1)["lsigma"])) 

### 258 out1["AIC_1"] < - md1$AIC 

### 259 out1["LL_1"] < - md1$LogLik 

### 260 out1["R2_1"] < - md1$R2 

### 261 } 

### 262 263 ## Model 2 

### 264 md2 < - mvtnormfit(df, 2,lP.cov) 

### 265 if (!is(md2, "try - error")) { 

### 266 out1["intercept_2"] < - unname(coef(md2)["intercept"]) 

### 267 out1["Field_2"] < - unname(coef(md2)["NearField"]) 

### 268 out1["m.lV_2"] < - unname(coef(md2)["m.lV"]) 

### 269 out1["sigma_2"] < - exp(unname(coef(md2)[" lsigma"])) 

### 270 out1["AIC_2"] < - md2$AIC 

### 271 out1["LL_2"] < - md2$LogLik 

### 272 out1["R2_2"] < - md2$R2 

### 273 } 

### 274 275 ## Model 3 

### 276 md3 < - mvtnormfit(df, 3, lP.cov) 

### 277 if (!is(md3, "try - error")){ 

### 278 out1["intercept_3"] < - unname(coef(md3)["intercept"]) 

### 279 out1["Field_3"] < - unname(coef(md3)["NearField"]) 

### 280 out1["m.lR_3"] < - unname(coef(md3)["m.lR"]) 

### 281 out1["sigma_3"] < - exp(unname(coef(md3)["lsigma"])) 

### 282 out1["AIC_3"] < - md3$AIC 

### 283 out1["LL_3"] < - md3$LogLik 

### 284 out1["R2_3"] < - md3$R2 

### 285 } 

### 286 287 ## Model 4 

### 288 md4 < - mvtnormfit(df, 4, lP.cov) 

### 289 if (!is(md4, "try - error")){ 

### 290 out1["intercept_4"] < - unname(coef(md4)["intercept"]) 

### 291 out1["Field_4"] < - unname(coef(md4)["NearField"]) 

### 292 out1["m.luR_4"] < - unname(coef(m d4)["m.luR"]) 

### 293 out1["m.lV_4"] < - unname(coef(md4)["m.lV"]) 

### 294 out1["sigma_4"] < - exp(unname(coef(md4)["lsigma"])) 

### 295 out1["AIC_4"] < - md4$AIC 

### 296 out1["LL_4"] < - md4$LogLik 

### 297 out1["R2_4"] < - md4$R2 

### 298 } 

### 299 ## out2 < - influence.peddling (1, df, md4, lP.cor) 

### 300 301 ## Model 5 

### 302 md5 < - mvtnormfit(df, 5, lP.cov) 

### 303 if (!is(md5, "try - error")) { 

### 304 out1["intercept_5"] < - unname(coef(md5)["intercept"]) 

### 305 out1["Field_5"] < - unname(coef(md5)["NearField"]) 

### 306 out1["m.lU_5"] < - unname( coef(md5)["m.lU"]) 

### 307 out1["sigma_5"] < - exp(unname(coef(md5)["lsigma"])) 

### 308 out1["AIC_5"] < - md5$AIC 

### 309 out1["LL_5"] < - md5$LogLik 

### 310out1["R2_5"] < - md5$R2 

### 311 } 

### 312 313 ## Model 6 

### 314 md6 < - mvtnormfit(df, 6, lP.cov) 

### 315 if (!is(md6, "try - error")){ 

### 316 out1["intercept_6"] < - unname(coef(md6)["intercept"]) 

### 317 out1["Field_6"] < - unname(coef(md6)["NearField"]) 

### 318 out1["m.luU_6"] < - unname(coef(md6)["m.luU"]) 

### 319 out1["m.lV_6"] < - unname(coef(md6)["m.lV"]) 

### 320 out1["sigma_6"] < - exp(unname(coef(md6)["lsigma"])) 

### 321 out1["AIC_6"] < - md6$AIC 

### 322 out1["LL_6"] < - md6$LogLik 

### 323 out1["R2_6"] < - md6$R2 

### 324 } 

### 325 ## out3 < - influence.peddling(2, df, md6, lP.cor) 

### 326 327 ## Model 7 

### 328 md7 < - mvtnormfit(df, 7, lP. cov) 

### 329 if (!is(md7, "try - error")) { 

### 330 out1["intercept_7"] < - unname(coef(md7)["intercept"]) 

### 331 out1["Field_7"] < - unname(coef(md7)["NearField"]) 

### 332 out1["m.luR_7"] < - unname(coef(md7)["m.luR"]) 

### 333 out1["m.luU_7"] < - unname(coef(md7)["m.luU"] ) 

### 334 out1["m.lV_7"] < - unname(coef(md7)["m.lV"]) 

### 335 out1["sigma_7"] < - exp(unname(coef(md7)["lsigma"])) 

### 336 out1["AIC_7"] < - md7$AIC 

### 337 out1["LL_7"] < - md7$LogLik 

### 338 out1["R2_7"] < - md7$R2 

### 339 } 

### 340 341 ## Model 8 

### 342 md8 < - mvtnormfit(df, 8, lP.cov) 

### 343 if (!is(md8, "try - error")) { 

### 344 out1["intercept_8"] < - unname(coef(md8)["intercept"]) 

### 345 out1["Field_8"] < - unname(coef(md8)["NearField"]) 

### 346 out1["m.lR_8"] < - unname(coef(md8)["m.lR"]) 

### 347 out1["m.lU_8"] < - unnam e(coef(md8)["m.lU"]) 

### 348 out1["sigma_8"] < - exp(unname(coef(md8)["lsigma"])) 

### 349 out1["AIC_8"] < - md8$AIC 

### 350 out1["LL_8"] < - md8$LogLik 

### 351 out1["R2_8"] < - md8$R2 

### 352 } 

### 353 354 ## Model 9 

### 355 md9 < - mvtnormfit(df, 9, lP.cov) 

### 356 if (!is(md9, "try - e rror")){ 

### 357 out1["intercept_9"] < - unname(coef(md9)["intercept"]) 

### 358 out1["Field_9"] < - unname(coef(md9)["NearField"]) 

### 359 out1["m.luR_9"] < - unname(coef(md9)["m.luR"]) 

### 360 out1["m.luU_9"] < - unname(coef(md9)["m.luU"]) 

### 361 out1["m.lV_9"] < - unn ame(coef(md9)["m.lV"]) 

### 362 out1["m.luR.N_9"] < - unname(coef(md9)["m.luR.N"]) 

### 363 out1["sigma_9"] < - exp(unname(coef(md9)["lsigma"])) 

### 364 out1["AIC_9"] < - md9$AIC 

### 365 out1["LL_9"] < - md9$LogLik 

### 366 out1["R2_9"] < - md9$R2 

### 367 } 

### 368 369 ## Model 10 

### 370 md10 < - mvtnormfit(df, 10, lP.cov) 

### 371 if (!is(md10, "try - error")){ 

### 372 out1["intercept_10"] < - unname(coef(md10)["intercept"]) 

### 373out1["Field_10"] < - unname(coef(md10)["NearField"]) 

### 374 out1["m.luR_10"] < - unname(coef(md10)["m.luR"]) 

### 375 out1 ["m.luU_10"] < - unname(coef(md10)["m.luU"]) 

### 376 out1["m.lV_10"] < - unname(coef(md10)["m.lV"]) 

### 377 out1["m.luU.N_10"] < - unname(coef(md10)["m.luU.N"]) 

### 378 out1["sigma_10"] < - exp(unname(coef(md10)["lsigma"])) 

### 379 out1["AIC_10"] < - md10$AIC 

### 380 out 1["LL_10"] < - md10$LogLik 

### 381 out1["R2_10"] < - md10$R2 

### 382 } 

### 383 384 ## Model 11 

### 385 md11 < - mvtnormfit(df, 11, lP.cov) 

### 386 if (!is(md11, "try - error")){ 

### 387 out1["intercept_11"] < - unname(coef(md11)["intercept"]) 

### 388 out1["Field_11"] < - unname(coef(md11 )["NearField"]) 

### 389 out1["m.luR_11"] < - unname(coef(md11)["m.luR"]) 

### 390 out1["m.luU_11"] < - unname(coef(md11)["m.luU"]) 

### 391 out1["m.lV_11"] < - unname(coef(md11)["m.lV"]) 

### 392 out1["m.lV.N_11"] < - unname(coef(md11)["m.lV.N"]) 

### 393 out1["sigma_11"] < - exp(unname(coef(md11)["lsigma"])) 

### 394 out1["AIC_11"] < - md11$AIC 

### 395 out1["LL_11"] < - md11$LogLik 

### 396 out1["R2_11"] < - md11$R2 

### 397 } 

### 398 399 ## Model 12 

### 400 md12 < - mvtnormfit(df, 12, lP.cov) 

### 401 if (!is(md12, "try - error")){ 

### 402 out1["intercept_12"] < - unname(coef(md12)["intercept"]) 

### 403 out1["Field_12"] < - unname(coef(md12)["NearField"]) 

### 404 out1["m.luR_12"] < - unname(coef(md12)["m.luR"]) 

### 405 out1["m.luU_12"] < - unname(coef(md12)["m.luU"]) 

### 406 out1["m.lV_12"] < - unn ame(coef(md12)["m.lV"]) 

### 407 out1["m.luR.N_12"] < - unname(coef(md12)["m.luR.N"]) 

### 408 out1["m.luU.N_12"] < - unname(coef(md12)["m.luU.N"]) 

### 409 out1["m.lV.N_12"] < - unname(coef(md12)["m.lV.N"]) 

### 410 out1["sigma_12"] < - exp(unname(coef(md12)["lsigma"])) 

### 411 out1["AIC_12"] < - md12$AIC 

### 412 out1["LL_12"] < - md12$LogLik 

### 413 out1["R2_12"] < - md12$R2 

### 414 } 

### 415 ## Model 13 

### 416 md13 < - mvtnormfit(df, 13, lP.cov) 

### 417 if (!is(md13, "try - error")){ 

### 418 out1["intercept_13"] < - unname(coef(md13)["intercept"]) 

### 419 out1["Field_13"] < - unname(coef(md13)["NearField"]) 

### 420 out1["m.luR_13"] < - unname(coef(md13)["m.luR"]) 

### 421 out1["m.luU_13"] < - unname(coef(md13)["m.luU"]) 

### 422 out1["m.lV_13"] < - unname(coef(md13)["m.lV"]) 

### 423 out1["m.luR.N_13"] < - unname(coef(md13)["m.luR.N"]) 

### 424 out1["m.luU.N_13"] < - unname(coef(md13)["m.luU.N"]) 

### 425 out1["sigma_13"] < - exp(unname(coef(md13)["lsigma"])) 

### 426 out1["AIC_13"] < - md13$AIC 

### 427 out1["LL_13"] < - md13$LogLik 

### 428 out1["R2_13"] < - md13$R2 

### 429 } 

### 430 ## Model 14 

### 431 md14 < - mvtnormfit(df, 14, lP.cov) 

### 432 if (!is(md14, "try - error")){ 

### 433 out1["intercept_14"] < - unname(coef(md14)["intercept"]) 

### 434 out1["Field_14"] < - unname(coef(md14)["NearField"]) 

### 435 out1["m.luR_14"] < - unname(coef(md14)["m.lu R"]) 

### 436out1["m.luU_14"] < - unname(coef(md14)["m.luU"]) 

### 437 out1["m.lV_14"] < - unname(coef(md14)["m.lV"]) 

### 438 out1["m.luU.N_14"] < - unname(coef(md14)["m.luU.N"]) 

### 439 out1["m.lV.N_14"] < - unname(coef(md14)["m.lV.N"]) 

### 440 out1["sigma_14"] < - exp(u nname(coef(md14)["lsigma"])) 

### 441 out1["AIC_14"] < - md14$AIC 

### 442 out1["LL_14"] < - md14$LogLik 

### 443 out1["R2_14"] < - md14$R2 

### 444 } 

### 445 446 ## Model 15 

### 447 md15 < - mvtnormfit(df, 15, lP.cov) 

### 448 if (!is(md15, "try - error")){ 

### 449 out1["intercept_15"] < - unname(coef(md15)["intercept"]) 

### 450 out1["Field_15"] < - unname(coef(md15)["NearField"]) 

### 451 out1["m.luR_15"] < - unname(coef(md15)["m.luR"]) 

### 452 out1["m.luU_15"] < - unname(coef(md15)["m.luU"]) 

### 453 out1["m.lV_15"] < - unname(coef(md15)["m.lV"]) 

### 454 out1["m.luR.N_15"] < - unname(coef(md15)["m.luR.N"]) 

### 455 out1["m.lV.N_15"] < - unname(coef(md15)["m.lV.N"]) 

### 456 out1["sigma_15"] < - exp(unname(coef(md15)["lsigma"])) 

### 457 out1["AIC_15"] < - md15$AIC 

### 458 out1["LL_15"] < - md15$LogLik 

### 459 out1["R2 _15"] < - md15$R2 

### 460 } 

### 461 462 ## Model 16 

### 463 md16 < - mvtnormfit(df, 16, lP.cov) 

### 464 if (!is(md16, "try - error")){ 

### 465 out1["intercept_16"] < - unname(coef(md16)["intercept"]) 

### 466 out1["Field_16"] < - unname(coef(md16)["NearField"]) 

### 467 out1["m.luR.N_ 16"] < - unname(coef(md16)["m.luR.N"]) 

### 468 out1["m.luU.N_16"] < - unname(coef(md16)["m.luU.N"]) 

### 469 out1["m.lV.N_16"] < - unname(coef(md16)["m.lV.N"]) 

### 470 out1["sigma_16"] < - exp(unname(coef(md16)["lsigma"])) 

### 471 out1["AIC_16"] < - md16$AIC 

### 472 out1[ "LL_16"] < - md16$LogLik 

### 473 out1["R2_16"] < - md16$R2 

### 474 } 

### 475 476 ## Model 17 

### 477 md17 < - mvtnormfit(df, 17, lP.cov) 

### 478 if (!is(md17, "try - error")){ 

### 479 out1["intercept_17"] < - unname(coef(md17)["intercept"]) 

### 480 out1["Field_17"] < - unname(coef(md17)["NearField"]) 

### 481 out1["m.luR.N_17"] < - unname(coef(md17)["m.luR.N"]) 

### 482 out1["m.luU.N_17"] < - unname(coef(md17)["m.luU.N"]) 

### 483 out1["sigma_17"] < - exp(unname(coef(md17)["lsigma"])) 

### 484 out1["AIC_17"] < - md17$AIC 

### 485 out1["L L_17"] < - md17$LogLik 

### 486 out1["R2_17"] < - md17$R2 

### 487 } 

### 488 ## Model 18 

### 489 md18 < - mvtnormfit(df, 18, lP.cov) 

### 490 if (!is(md18, "try - error")){ 

### 491 out1["intercept_18"] < - unname(coef(md18)["intercept"]) 

### 492 out1["Field_18"] < - unname(coef(md18)["Ne arField"]) 

### 493 out1["m.luR.N_18"] < - unname(coef(md18)["m.luR.N"]) 

### 494 out1["m.lV.N_18"] < - unname(coef(md18)["m.lV.N"]) 

### 495 out1["sigma_18"] < - exp(unname(coef(md18)["lsigma"])) 

### 496 out1["AIC_18"] < - md18$AIC 

### 497 out1["LL_18"] < - md18$LogLik 

### 498 out1["R2_18"] < - md18$R2 

### 499} 

### 500 501 ## Model 19 

### 502 md19 < - mvtnormfit(df, 19, lP.cov) 

### 503 if (!is(md19, "try - error")){ 

### 504 out1["intercept_19"] < - unname(coef(md19)["intercept"]) 

### 505 out1["Field_19"] < - unname(coef(md19)["NearField"]) 

### 506 out1["m.luU.N_19"] < - unname(coef(md19)["m.luU.N"]) 

### 507 out1["m.lV.N_19"] < - unname(coef(md19)["m.lV.N"]) 

### 508 out1["sigma_19"] < - exp(unname(coef(md19)["lsigma"])) 

### 509 out1["AIC_19"] < - md19$AIC 

### 510 out1["LL_19"] < - md19$LogLik 

### 511 out1["R2 _19"] < - md19$R2 

### 512 } 

### 513 514 ## Model 20 

### 515 md20 < - mvtnormfit(df, 20, lP.cov) 

### 516 if (!is(md20, "try - error")){ 

### 517 out1["intercept_20"] < - unname(coef(md20)["intercept"]) 

### 518 out1["Field_20"] < - unname(coef(md20)["NearField"]) 

### 519 out1["m.luR.N_ 20"] < - unname(coef(md20)["m.luR.N"]) 

### 520 out1["sigma_20"] < - exp(unname(coef(md20)["lsigma"])) 

### 521 out1["AIC_20"] < - md20$AIC 

### 522 out1["LL_20"] < - md20$LogLik 

### 523 out1["R2_20"] < - md20$R2 

### 524 } 

### 525 526 ## Model 21 

### 527 md21 < - mvtnormfit(df, 21, lP.co v) 

### 528 if (!is(md21, "try - error")){ 

### 529 out1["intercept_21"] < - unname(coef(md21)["intercept"]) 

### 530 out1["m.lR_21"] < - unname(coef(md21)["m.lR"]) 

### 531 out1["sigma_21"] < - exp(unname(coef(md21)["lsigma"])) 

### 532 out1["AIC_21"] < - md21$AIC 

### 533 out1["L L_21"] < - md21$LogLik 

### 534 out1["R2_21"] < - md21$R2 

### 535 } 

### 536 537 ## Model 22 

### 538 md22 < - mvtnormfit(df, 22, lP.cov) 

### 539 if (!is(md22, "try - error")) { 

### 540 out1["intercept_22"] < - unname(coef(md22)["intercept"]) 

### 541 out1["m.lU_22"] < - unname(coef(md22)["m.lU"]) 

### 542 out1["sigma_22"] < - exp(unname(coef(md22)["lsigma"])) 

### 543 out1["AIC_22"] < - md22$AIC 

### 544 out1["LL_22"] < - md22$LogLik 

### 545 out1["R2_22"] < - md22$R2 

### 546 } 

### 547 548 ## Model Comparisons: record the null - hypothetical delt a - R2 and 

### 549 ## delta log - likelihood for the corresponding two models. Construct 

### 550 ## a new null - hypothetical dataset for each test. 

### 551 552 ## Try each set of fits 3 times before giving up 

### 553 df0 < - df 

### 554 ## newDataset from md0 

### 555 if (!is(md0, "try - err or")) { 

### 556 muhat < - rep(coef(md0)["intercept"], nrow(df)) 

### 557 Sigma < - lP.cov + 

### 558 diag(exp(2*coef(md0)["lsigma"]), 

### 559 nrow(lP.cov), ncol(lP.cov)) 

### 560 Nrep < - 0 

### 561 dorep < - TRUE 

### 562while(dorep && Nrep < 3){ 

### 563 Nrep < - Nrep + 1 

### 564 lP < - rmvnorm(1,muhat,sigma=Sigma) 

### 565 df0$lP < - lP[1,] 

### 566 md0a < - mvtnormfit(df0, 0, lP.cov) 

### 567 md1a < - mvtnormfit(df0, 1, lP.cov) 

### 568 if (!is(md0a, "try - error") && !is(md1a,"try - error")) dore p < - FALSE 

### 569 } 

### 570 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 571 out1["dR0< - >1"] < - md1a$R2 ## md0a$R2 should be 0 

### 572 out1["dLL0< - >1"] < - md1a$LogLik - md0a$LogLik 

### 573 } 

### 574 575 Nrep < - 0 

### 576 dorep < - TRUE 

### 577 while(dorep & & Nrep < 3){ 

### 578 Nrep < - Nrep + 1 

### 579 lP < - rmvnorm(1,muhat,sigma=Sigma) 

### 580 df0$lP < - lP[1,] 

### 581 md0a < - mvtnormfit(df0, 0, lP.cov) 

### 582 md1a < - mvtnormfit(df0, 21, lP.cov) 

### 583 if (!is(md0a, "try - error") && !is(md1a,"try - error")) do rep < - FALSE 

### 584 } 

### 585 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 586 out1["dR0< - >21"] < - md1a$R2 

### 587 out1["dLL0< - >21"] < - md1a$LogLik - md0a$LogLik 

### 588 } 

### 589 590 Nrep < - 0 

### 591 dorep < - TRUE 

### 592 while(dorep && Nrep < 3){ 

### 593 Nrep < - Nrep + 1 

### 594 lP < - rmvnorm(1,muhat,sigma=Sigma) 

### 595 df0$lP < - lP[1,] 

### 596 md0a < - mvtnormfit(df0, 0, lP.cov) 

### 597 md1a < - mvtnormfit(df0, 22, lP.cov) 

### 598 if (!is(md0a, "try - error") && !is(md1a,"try - error")) dorep < - FALSE 

### 599 } 

### 600 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 601 out1["dR0< - >22"] < - md1a$R2 

### 602 out1["dLL0< - >22"] < - md1a$LogLik - md0a$LogLik 

### 603 } 

### 604 } 

### 605 ## New datasets from md1 

### 606 if (!is(md1, "try - error")) { 

### 607 muhat < - coef( md1)["intercept"] + coef(md1)["NearField"]*df$NearField 

### 608 Sigma < - lP.cov + diag(exp(2*coef(md1)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 609 Nrep < - 0 

### 610 dorep < - TRUE 

### 611 while(dorep && Nrep < 3){ 

### 612 Nrep < - Nrep + 1 

### 613 lP < - rmvnorm (1, muhat, sigma=Sigma) 

### 614 df0$lP < - lP[1,] 

### 615 md0a < - mvtnormfit(df0, 1, lP.cov) 

### 616 md1a < - mvtnormfit(df0, 2, lP.cov) 

### 617 } 

### 618 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 619 out1["dR1< - >2"] < - md1a$R2 - md0a$R2 

### 620 out1["dLL1< - >2"] < - md1a$LogLik - md0a$LogLik 

### 621 } 

### 622 623 Nrep < - 0 

### 624 dorep < - TRUE 

### 625while(dorep && Nrep < 3){ 

### 626 Nrep < - Nrep + 1 

### 627 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 628 df0$lP < - lP[1,] 

### 629 md0a < - mvtnormf it(df0, 1, lP.cov) 

### 630 md1a < - mvtnormfit(df0, 3, lP.cov) 

### 631 } 

### 632 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 633 out1["dR1< - >3"] < - md1a$R2 - md0a$R2 

### 634 out1["dLL1< - >3"] < - md1a$LogLik - md0a$LogLik 

### 635 } 

### 636 637 Nrep < - 0 

### 638 dorep < - TRUE 

### 639 while(dorep && Nrep < 3){ 

### 640 Nrep < - Nrep + 1 

### 641 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 642 df0$lP < - lP[1,] 

### 643 md0a < - mvtnormfit(df0, 1, lP.cov) 

### 644 md1a < - mvtnormfit(df0, 4, lP.cov) 

### 645 } 

### 646 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 647 out1["dR1< - >4"] < - md1a$R2 - md0a$R2 

### 648 out1["dLL1< - >4"] < - md1a$LogLik - md0a$LogLik 

### 649 } 

### 650 651 Nrep < - 0 

### 652 dorep < - TRUE 

### 653 while(dorep && Nrep < 3){ 

### 654 Nrep < - Nrep + 1 

### 655 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 656 df0$lP < - lP[1,] 

### 657 md0a < - mvtnormfit(df0, 1, lP.cov) 

### 658 md1a < - mvtnormfit(df0, 5, lP.cov) 

### 659 } 

### 660 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 661 out1["dR1< - >5"] < - md1a$R2 - md0a$R2 

### 662 out1["dLL1< - >5"] < - md1a$LogLik - md0a$LogLik 

### 663 } 

### 664 665 Nrep < - 0 

### 666 dorep < - TRUE 

### 667 while(dorep && Nrep < 3){ 

### 668 Nrep < - Nrep + 1 

### 669 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 670 df0$lP < - lP[1, ] 

### 671 md0a < - mvtnormfit(df0, 1, lP.cov) 

### 672 md1a < - mvtnormfit(df0, 6, lP.cov) 

### 673 } 

### 674 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 675 out1["dR1< - >6"] < - md1a$R2 - md0a$R2 

### 676 out1["dLL1< - >6"] < - md1a$LogLik - md0a$LogLik 

### 677 } 

### 678 679 Nrep < - 0 

### 680 dorep < - TRUE 

### 681 while(dorep && Nrep < 3){ 

### 682 Nrep < - Nrep + 1 

### 683 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 684 df0$lP < - lP[1,] 

### 685 md0a < - mvtnormfit(df0, 1, lP.cov) 

### 686 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 687 } 

### 688if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 689 out1["dR1< - >12"] < - md1a$R2 - md0a$R2 

### 690 out1["dLL1< - >12"] < - md1a$LogLik - md0a$LogLik 

### 691 } 

### 692 } 

### 693 694 ## New datasets from M2 

### 695 if (!is(md2, "try - error")) { 

### 696 muhat < - coef(md2)["intercept"] + coef(md2)["NearField"]*df$NearField + 

### 697 coef(md2)["m.lV"] * df$lV 

### 698 Sigma < - lP.cov + diag(exp(2*coef(md2)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 699 700 Nrep < - 0 

### 701 dorep < - TRUE 

### 702 while(dorep && Nrep < 3){ 

### 703 Nrep < - Nrep + 1 

### 704 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 705 df0$lP < - lP[1,] 

### 706 md0a < - mvtnormfit(df0, 2, lP.cov) 

### 707 md1a < - mvtnormfit(df0, 4, lP.cov) 

### 708 } 

### 709 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 710 out1["dR2< - >4"] < - md1a$R2 - md0a$R2 

### 711 out1["dLL2< - >4"] < - md1a$LogLik - md0a$LogLik 

### 712 } 

### 713 714 Nrep < - 0 

### 715 dorep < - TRUE 

### 716 while(dorep && Nrep < 3){ 

### 717 Nrep < - Nrep + 1 

### 718 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 719 df0$lP < - lP[1,] 

### 720 md0a < - mvtnormfit(df0, 2, lP.cov) 

### 721 md1a < - mvtnormfit(df0, 6, lP.cov) 

### 722 } 

### 723 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 724 out1["dR2< - >6"] < - md1a$R2 - md0a$R2 

### 725 out1["dLL2< - >6"] < - md1a$LogLik - md0a$LogLik 

### 726 } 

### 727 } 

### 728 729 ## One new dataset from M3 

### 730 if (!is(md3, "try - error")) { 

### 731 muhat < - coef(md3)["intercept"] + coef(md3)["NearField"]*df$NearField + 

### 732 coef(md3)["m.lR"] * df$lR 

### 733 Sigma < - lP.cov + diag(exp(2*coef(md3)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 734 Nrep < - 0 

### 735 dorep < - TRUE 

### 736 while(dorep && Nrep < 3){ 

### 737 Nrep < - Nrep + 1 

### 738 lP < - rmvnorm(1, muhat, sigma= Sigma) 

### 739 df0$lP < - lP[1,] 

### 740 md0a < - mvtn ormfit(df0, 3, lP.cov) 

### 741 md1a < - mvtnormfit(df0, 4, lP.cov) 

### 742 } 

### 743 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 744 out1["dR3< - >4"] < - md1a$R2 - md0a$R2 

### 745 out1["dLL3< - >4"] < - md1a$LogLik - md0a$LogLik 

### 746 } 

### 747 } 

### 748 749 ## One new dataset from M4 

### 750 if (!is(md4, "try - error")){ 

### 751muhat < - coef(md4)["intercept"] + coef(md4)["NearField"]*df$NearField + 

### 752 coef(md4)["m.luR"] * df$ls.R + coef(md4)["m.lV"] * df$lV 

### 753 Sigma < - lP.cov + diag(exp(2*coef(md4)["lsigma" ]), nrow(lP.cov), ncol(lP.cov)) 

### 754 755 Nrep < - 0 

### 756 dorep < - TRUE 

### 757 while(dorep && Nrep < 3){ 

### 758 Nrep < - Nrep + 1 

### 759 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 760 df0$lP < - lP[1,] 

### 761 md0a < - mvtnormfit(df0, 4, lP.cov) 

### 762 md1a < - mvtnormfit(df0, 7, lP.cov) 

### 763 } 

### 764 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 765 out1["dR4< - >7"] < - md1a$R2 - md0a$R2 

### 766 out1["dLL4< - >7"] < - md1a$LogLik - md0a$LogLik 

### 767 } 

### 768 } 

### 769 770 ## One new dataset fro m M5 

### 771 if (!is(md5, "try - error")) { 

### 772 muhat < - coef(md5)["intercept"] + coef(md5)["NearField"]*df$NearField + 

### 773 coef(md5)["m.lU"] * df$lU 

### 774 Sigma < - lP.cov + diag(exp(2*coef(md5)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 775 776 Nrep < - 0 

### 777 dorep < - TRUE 

### 778 while(dorep && Nrep < 3){ 

### 779 Nrep < - Nrep + 1 

### 780 lP < - rmvnorm(1, muhat, sigma = Sigma) 

### 781 782 df0$lP < - lP[1,] 

### 783 md0a < - mvtnormfit(df0, 5, lP.cov) 

### 784 md1a < - mvtnormfit(df0, 6, lP.cov) 

### 785 } 

### 786 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 787 out1["dR5< - >6"] < - md1a$R2 - md0a$R2 

### 788 out1["dLL5< - >6"] < - md1a$LogLik - md0a$LogLik 

### 789 } 

### 790 } 

### 791 792 ## One new dataset from M6 

### 793 if (!is(md6, "try - error")) { 

### 794 muhat < - coef(md6)["intercept"] + coef(md6)["NearField"]*df$NearField + 

### 795 coef(md6)["m.luU"] * df$ls.U + coef(md6)["m.lV"] * df$lV 

### 796 Sigma < - lP.cov + diag(exp(2*coef(md6)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 797 798 Nrep < - 0 

### 799 dorep < - TRU E 

### 800 while(dorep && Nrep < 3){ 

### 801 Nrep < - Nrep + 1 

### 802 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 803 804 df0$lP < - lP[1,] 

### 805 md0a < - mvtnormfit(df0, 6, lP.cov) 

### 806 md1a < - mvtnormfit(df0, 7, lP.cov) 

### 807 } 

### 808 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 809 out1["dR6< - >7"] < - md1a$R2 - md0a$R2 

### 810 out1["dLL6< - >7"] < - md1a$LogLik - md0a$LogLik 

### 811 } 

### 812 } 

### 813 814## New dataset from M7 

### 815 if (!is(md7, "try - error")) { 

### 816 muhat < - coef(md7)["intercept"] + coef(md7)["NearField"]*df$NearField + 

### 817 coef(md7)["m.luR"] * df$ls.R + 

### 818 coef(md7)["m.luU"] * df$ls.U + coef(md7)["m.lV"] * df$lV 

### 819 Sigma < - lP.cov + diag(exp(2*coef(md7)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 820 821 Nrep < - 0 

### 822 dorep < - TRUE 

### 823 while(dorep && Nrep < 3){ 

### 824 Nrep < - Nrep + 1 

### 825 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 826 827 df0$lP < - lP[1,] 

### 828 md0a < - mvtnormfit(df0, 7, lP.cov) 

### 829 md1a < - mvtnormfit(df0, 9, lP. cov) 

### 830 } 

### 831 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 832 out1["dR7< - >9"] < - md1a$R2 - md0a$R2 

### 833 out1["dLL7< - >9"] < - md1a$LogLik - md0a$LogLik 

### 834 } 

### 835 836 Nrep < - 0 

### 837 dorep < - TRUE 

### 838 while(dorep && Nrep < 3){ 

### 839 Nrep < - Nrep + 1 

### 840 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 841 842 df0$lP < - lP[1,] 

### 843 md0a < - mvtnormfit(df0, 7, lP.cov) 

### 844 md1a < - mvtnormfit(df0, 10, lP.cov) 

### 845 } 

### 846 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 847 out1["dR7< - >10"] < - md1a$R2 - md0a$R2 

### 848 out1["dLL7< - >10"] < - md1a$LogLik - md0a$LogLik 

### 849 } 

### 850 851 Nrep < - 0 

### 852 dorep < - TRUE 

### 853 while(dorep && Nrep < 3){ 

### 854 Nrep < - Nrep + 1 

### 855 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 856 857 df0$lP < - lP[1,] 

### 858 md0a < - mvtnormfit(df0, 7, lP.cov) 

### 859 md1a < - mvtnormfit(df0, 11, lP.cov) 

### 860 } 

### 861 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 862 out1["dR7< - >11"] < - md1a$R2 - md0a$R2 

### 863 out1["dLL7< - >11"] < - md1 a$LogLik - md0a$LogLik 

### 864 } 

### 865 866 Nrep < - 0 

### 867 dorep < - TRUE 

### 868 while(dorep && Nrep < 3){ 

### 869 Nrep < - Nrep + 1 

### 870 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 871 872 df0$lP < - lP[1,] 

### 873 md0a < - mvtnormfit(df0, 7, lP.cov) 

### 874 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 875 } 

### 876 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 877out1["dR7< - >12"] < - md1a$R2 - md0a$R2 

### 878 out1["dLL7< - >12"] < - md1a$LogLik - md0a$LogLik 

### 879 } 

### 880 } 

### 881 882 ## One new dataset from M8 

### 883 if (!is(md8, "try - error")) { 

### 884 muhat < - coef(md8)["intercept"] + coef(md8)["NearField"]*df$NearField + 

### 885 coef(md8)["m.lR"] * df$lR + coef(md8)["m.lU"] * df$lU 

### 886 Sigma < - lP.cov + diag(exp(2*coef(md8)["lsigma"]), nrow(lP.cov), ncol (lP.cov)) 

### 887 Nrep < - 0 

### 888 dorep < - TRUE 

### 889 while(dorep && Nrep < 3){ 

### 890 Nrep < - Nrep + 1 

### 891 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 892 893 df0$lP < - lP[1,] 

### 894 md0a < - mvtnormfit(df0, 8, lP.cov) 

### 895 md1a < - mvtnormfit(df0, 7, lP.cov) 

### 896 } 

### 897 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 898 out1["dR8< - >7"] < - md1a$R2 - md0a$R2 

### 899 out1["dLL8< - >7"] < - md1a$LogLik - md0a$LogLik 

### 900 } 

### 901 } 

### 902 ## One new dataset from M9 

### 903 if (!is(md9, "try - error")) { 

### 904 muhat < - coef(md9)["intercept"] + coef(md9)["NearField"]*df$NearField + 

### 905 coef(md9)["m.luR"] * df$ls.R + coef(md9)["m.luU"] * df$ls.U + 

### 906 coef(md9)["m.lV"] * df$lV + coef(md9)["m.luR.N"] * df$NearField * df$ls.R 

### 907 Sigma < - lP.cov + diag(exp(2*coef(md9)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 908 909 Nrep < - 0 

### 910 dorep < - TRUE 

### 911 while(dorep && Nrep < 3){ 

### 912 Nrep < - Nrep + 1 

### 913 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 914 915 df0$lP < - lP[1,] 

### 916 md0a < - mvtn ormfit(df0, 9, lP.cov) 

### 917 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 918 } 

### 919 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 920 out1["dR9< - >12"] < - md1a$R2 - md0a$R2 

### 921 out1["dLL9< - >12"] < - md1a$LogLik - md0a$LogLik 

### 922 } 

### 923 } 

### 924 925 ## One new dataset from M10 

### 926 if (!is(md10, "try - error")) { 

### 927 muhat < - coef(md10)["intercept"] + coef(md10)["NearField"]*df$NearField + 

### 928 coef(md10)["m.luR"] * df$ls.R + coef(md10)["m.luU"] * df$ls.U + 

### 929 coef(md10)["m.lV"] * df$lV + coef(md10)["m.luU.N"] * df$NearField * df$ls.U 

### 930 Sigma < - lP.cov + diag(exp(2*coef(md10)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 931 932 Nrep < - 0 

### 933 dorep < - TRUE 

### 934 while(dorep && Nrep < 3){ 

### 935 N rep < - Nrep + 1 

### 936 lP < - rmvnorm(1, muhat, sigma=lP.cov + diag(exp(2*coef(md10)["lsigma"]), 

### 937 nrow(lP.cov), ncol(lP.cov))) 

### 938 939 df0$lP < - lP[1,] 

### 940md0a < - mvtnormfit(df0, 10, lP.cov) 

### 941 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 942 } 

### 943 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 944 out1["dR10< - >12"] < - md1a$R2 - md0a$R2 

### 945 out1["dLL10< - >12"] < - md1a$LogLik - md0a$LogLik 

### 946 } 

### 947 } 

### 948 949 ## One new dataset from M11 

### 950 if (!is(md11, "try - error")) { 

### 951 muhat < - coef(md11)["intercept"] + coef(md11)["NearField"]*df$NearField + 

### 952 coef(md11)["m.luR"] * df$ls.R + coef(md11)["m.luU"] * df$ls.U + 

### 953 coef(md11)["m.lV"] * df$lV + coef(md11)["m.lV.N"] * df$NearField * df$lV 

### 954 Sigma < - lP.co v + diag(exp(2*coef(md11)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 955 Nrep < - 0 

### 956 dorep < - TRUE 

### 957 while(dorep && Nrep < 3){ 

### 958 Nrep < - Nrep + 1 

### 959 lP < - rmvnorm(1, muhat, sigma=Sigma) 

### 960 961 df0$lP < - lP[1,] 

### 962 md0a < - mvtnormfit(df0, 11, lP.cov) 

### 963 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 964 } 

### 965 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 966 out1["dR11< - >12"] < - md1a$R2 - md0a$R2 

### 967 out1["dLL11< - >12"] < - md1a$LogLik - md0a$LogLik 

### 968 } 

### 969 } 

### 970 971 ## dataset from M16 

### 972 if (!is(md16, "try - error")){ 

### 973 muhat < - coef(md16)["intercept"] + 

### 974 df$NearField * (coef(md16)["NearField"] + 

### 975 coef(md16)["m.luR.N"] * df$ls.R + 

### 976 coef(md16)["m.luU .N"] * df$ls.U + 

### 977 coef(md16)["m.lV.N"] * df$lV) 

### 978 Sigma < - lP.cov + 

### 979 diag(exp(2*coef(md16)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 980 981 Nrep < - 0 

### 982 dorep < - TRUE 

### 983 while(dorep && Nrep < 3){ 

### 984 Nrep < - Nrep + 1 

### 985 lP < - rmvnorm(1, muhat, 

### 986 sigma=Sigma) 

### 987 df0$lP < - lP[1,] 

### 988 md0a < - mvtnormfit(df0, 16, lP.cov) 

### 989 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 990 } 

### 991 if (!is(md0a, "try - error") && !is(md1a,"try - error") ) { 

### 992 out1["dR16< - >12"] < - md1a$R2 - md0a$R2 

### 993 out1["dLL16< - >12"] < - md1a$LogLik - md0a$LogLik 

### 994 } 

### 995 } 

### 996 997 ## dataset from M18 

### 998 if (!is(md18, "try - error")){ 

### 999 muhat < - coef(md18)["intercept"] + 

### 1000 df$NearField * (coef(md1 8)["NearField"] + 

### 1001 coef(md18)["m.luR.N"] * df$ls.R + 

### 1002 coef(md18)["m.lV.N"] * df$lV) 

### 1003Sigma < - lP.cov + 

### 1004 diag(exp(2*coef(md18)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 1005 1006 Nrep < - 0 

### 1007 dor ep < - TRUE 

### 1008 while(dorep && Nrep < 3){ 

### 1009 Nrep < - Nrep + 1 

### 1010 lP < - rmvnorm(1, muhat, 

### 1011 sigma=Sigma) 

### 1012 df0$lP < - lP[1,] 

### 1013 md0a < - mvtnormfit(df0, 18, lP.cov) 

### 1014 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 1015 } 

### 1016 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 1017 out1["dR18< - >12"] < - md1a$R2 - md0a$R2 

### 1018 out1["dLL18< - >12"] < - md1a$LogLik - md0a$LogLik 

### 1019 } 

### 1020 } 

### 1021 1022 ## dataset from M20 

### 1023 if (!is(md20, "try - error")){ 

### 1024 muhat < - c oef(md20)["intercept"] + 

### 1025 df$NearField * (coef(md20)["NearField"] + 

### 1026 coef(md20)["m.luR.N"] * df$ls.R) 

### 1027 Sigma < - lP.cov + 

### 1028 diag(exp(2*coef(md20)["lsigma"]), nrow(lP.cov), ncol(lP.cov)) 

### 1029 1030 Nrep < - 0 

### 1031 d orep < - TRUE 

### 1032 while(dorep && Nrep < 3){ 

### 1033 Nrep < - Nrep + 1 

### 1034 lP < - rmvnorm(1, muhat, 

### 1035 sigma=) 

### 1036 df0$lP < - lP[1,] 

### 1037 md0a < - mvtnormfit(df0, 20, lP.cov) 

### 1038 md1a < - mvtnormfit(df0, 12, lP.cov) 

### 1039 } 

### 1040 if (!is(md0a, "try - error") && !is(md1a,"try - error")) { 

### 1041 out1["dR20< - >12"] < - md1a$R2 - md0a$R2 

### 1042 out1["dLL20< - >12"] < - md1a$LogLik - md0a$LogLik 

### 1043 } 

### 1044 } 

### 1045 1046 out1 

### 1047 } 

### 1048 1049 analysisC < - cmpfun(analysis) 

### 1050 1051 df < - data.frame(NearField = pred.data$NearField, 

### 1052 lR = log(pred.data$s.R) + lV.mx[1,], 

### 1053 lU = log(pred.data$s.U) + lV.mx[1,], 

### 1054 ls.R = log(pred.data$s.R), 

### 1055 ls.U = log(pred.data$s.U), 

### 1056 CAS = pred. data$CAS, 

### 1057 Name = pred.data$Name, 

### 1058 MW = pred.data$MW, 

### 1059 lP = lP.mx[1,], 

### 1060 lV = lV.mx[1,], 

### 1061 vars = lP.vars) 

### 1062 1063 out < - mclapply(RegSampIndx, function(j) { 

### 1064 df$lV < - lV .mx[j,] 

### 1065 df$lP < - lP.mx[j,] 

### 1066analysisC(df) 

### 1067 },mc.preschedule=FALSE, mc.cores=6) 

### 1068 Betas < - do.call("rbind", out) 

### 1069 save(Betas, file="AnalysisOmnibusFitssave.RData") 

### 1070 } 

### 1071 1072 if (!FittingModels) load("AnalysisOmnibusFitssave.RData") 

### 1073 1074 doSummarizeAnalysis < - function(Betas) { 

### 1075 1076 ## 

### 1077 ## Which model works best overall? 

### 1078 ## Truncate Betas to the first 1500 values. The rest are all NA 

### 1079 AICs < - Betas[,grep("AIC",colnames(Betas))] 

### 1080 1081 1082 ## Summary of the Betas. Primarily of interest to track NA's. 

### 1083 writeLines( c("", 

### 1084 1085 "========================================================================", 

### 1086 "Summary of All Model Fits")) 

### 1087 print(summary(Betas)) 

### 1088 1089 ## Mean AIC 

### 1090 writeLines(c("", 

### 1091 1092 "=========================== =============================================", 

### 1093 "Mean AIC")) 

### 1094 z < - apply(AICs,2,mean, na.rm=TRUE) 

### 1095 print(z - z[which.min(z)], digits=3) 

### 1096 1097 writeLines(c("", 

### 1098 1099 "==============================================================================", 

### 1100 "Distribution of which models have best AIC")) 

### 1101 1102 print(table(apply(AICs,1,which.min) - 1)) 

### 1103 1104 1105 ## 

### 1106 ## Another way to look at this is with mean r anks: Rank 

### 1107 ## the AICs across models within replicates, then 

### 1108 ## average the ranks across replicates. This will let 

### 1109 ## us see, for instance, if a given model is often 

### 1110 ## second - best 

### 1111 1112 writeLines(c("", 

### 1113 1114 "============================== ===============================================", 

### 1115 "Mean Rank of each model's AIC")) 

### 1116 1117 tmp < - apply(AICs,1,rank, na.last="keep") 

### 1118 mnrank < - apply(tmp,1,mean, na.rm=TRUE) 

### 1119 print(mnrank) 

### 1120 1121 1122 1123 AICsummaries < - 

### 1124 rbind( 

### 1125 "0 < - > 1" = mysummary(AICs[,"AIC_0"] - AICs[,"AIC_1"]), 

### 1126 "0 < - > 21" = mysummary(AICs[,"AIC_0"] - AICs[,"AIC_21"]), 

### 1127 "0 < - > 22" = mysummary(AICs[,"AIC_0"] - AICs[,"AIC_22"]), 

### 1128 "1 < - > 2" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_2"]), 

### 1129"1 < - > 3" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_3"]), 

### 1130 "1 < - > 4" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_4"]), 

### 1131 "1 < - > 5" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_5"]), 

### 1132 "1 < - > 6" = mysummary(AICs[,"AIC_1"] - AICs[,"AI C_6"]), 

### 1133 "1 < - > 12" = mysummary(AICs[,"AIC_1"] - AICs[,"AIC_12"]), 

### 1134 "2 < - > 4" = mysummary(AICs[,"AIC_2"] - AICs[,"AIC_4"]), 

### 1135 "2 < - > 6" = mysummary(AICs[,"AIC_2"] - AICs[,"AIC_6"]), 

### 1136 "3 < - > 4" = mysummary(AICs[,"AIC_3"] - AICs[,"AIC_4"]), 

### 1137 "4 < - > 7" = mysummary(AICs[,"AIC_4"] - AICs[,"AIC_7"]), 

### 1138 "5 < - > 6" = mysummary(AICs[,"AIC_5"] - AICs[,"AIC_6"]), 

### 1139 "6 < - > 7" = mysummary(AICs[,"AIC_6"] - AICs[,"AIC_7"]), 

### 1140 "7 < - > 12" = mysummary(AICs[," AIC_7"] - AICs[,"AIC_12"]), 

### 1141 "7 < - > 9" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_9"]), 

### 1142 "7 < - > 10" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_10"]), 

### 1143 "7 < - > 11" = mysummary(AICs[,"AIC_7"] - AICs[,"AIC_11"]), 

### 1144 "8 < - > 7" = mysu mmary(AICs[,"AIC_8"] - AICs[,"AIC_7"]), 

### 1145 "9 < - > 12" = mysummary(AICs[,"AIC_9"] - AICs[,"AIC_12"]), 

### 1146 "10 < - > 12" = mysummary(AICs[,"AIC_10"] - AICs[,"AIC_12"]), 

### 1147 "11 < - > 12" = mysummary(AICs[,"AIC_11"] - AICs[,"AIC_12"]), 

### 1148 "16 < - > 12" = mysummary(AICs[,"AIC_16"] - AICs[,"AIC_12"]), 

### 1149 "18 < - > 12" = mysummary(AICs[,"AIC_18"] - AICs[,"AIC_12"]), 

### 1150 "20 < - > 12" = mysummary(AICs[,"AIC_20"] - AICs[,"AIC_12"]) 

### 1151 ) 

### 1152 1153 1154 writeLines(c("", 

### 1155 1156 "====== =======================================================================", 

### 1157 "Distributions of AIC differences for Model Comparisons")) 

### 1158 print(AICsummaries, digits=3) 

### 1159 1160 ## Test Results. Compare log - likelihood differences with 

### 1161 ## differences generated from the simpler model for all these 

### 1162 ## comparisons. Tabulate the fraction of replicates for which 

### 1163 ## LL_j - LL_k > dLLjk 

### 1164 1165 ## Also, What is the change in R2 for e ach treatment? Put them both in 

### 1166 ## the same table 

### 1167 1168 1169 1170 Tests < - 

### 1171 matrix(c( 

### 1172 mycomps(Betas,c(0,1)), 

### 1173 mycomps(Betas,c(0,21)), 

### 1174 mycomps(Betas,c(0,22)), 

### 1175 mycomps(Betas,c(1,2)), 

### 1176 mycomps(Betas,c(1 ,3)), 

### 1177 mycomps(Betas,c(1,4)), 

### 1178 mycomps(Betas,c(1,5)), 

### 1179 mycomps(Betas,c(1,6)), 

### 1180 mycomps(Betas,c(1,12)), 

### 1181 mycomps(Betas,c(2,4)), 

### 1182 mycomps(Betas,c(2,6)), 

### 1183 mycomps(Betas,c(3,4) ), 

### 1184 mycomps(Betas,c(4,7)), 

### 1185 mycomps(Betas,c(5,6)), 

### 1186 mycomps(Betas,c(6,7)), 

### 1187 mycomps(Betas,c(7,9)), 

### 1188 mycomps(Betas,c(7,10)), 

### 1189 mycomps(Betas,c(7,11)), 

### 1190 mycomps(Betas,c(7,12)), 

### 1191 mycomps(Betas,c(8,7)), 

### 1192mycomps(Betas,c(9,12)), 

### 1193 mycomps(Betas,c(10,12)), 

### 1194 mycomps(Betas,c(11,12)), 

### 1195 mycomps(Betas,c(16,12)), 

### 1196 mycomps(Betas,c(18, 12)), 

### 1197 mycomps(Betas,c(20,12)) 

### 1198 ), 

### 1199 ncol=2, 

### 1200 dimnames=list(c( 

### 1201 "0 < - > 1","0 < - > 21","0 < - > 22", 

### 1202 "1 < - > 2","1 < - > 3","1 < - > 4","1 < - > 5","1 < - > 6","1 < - > 12", 

### 1203 "2 < - > 4","2 < - > 6", 

### 1204 "3 < - > 4", 

### 1205 "4 < - > 7", 

### 1206 "5 < - > 6", 

### 1207 "6 < - > 7", 

### 1208 "7 < - > 9","7 < - > 10","7 < - > 11","7 < - > 12", 

### 1209 "8 < - > 7", 

### 1210 "9 < - > 12","10 < - > 12","11 < - > 12", 

### 1211 "16 < - > 12", "18 < - > 12","20 < - > 12" 

### 1212 ), 

### 1213 c("P - value","delta R2")), 

### 1214 byrow=TRUE) 

### 1215 1216 1217 writeLines(c("", 

### 1218 1219 "=============================================================================", 

### 1220 "Significance of Model Comparisons")) 

### 1221 print(Tests) 

### 1222 1223 1224 ## -------------------------------------------------------- ----------------------- 

### 1225 ## Estimates 

### 1226 1227 writeLines(c("", 

### 1228 1229 "=============================================================================", 

### 1230 "Estimates of Model Parameters, with 95% CI","")) 

### 1231 for (i in 1:22) 

### 1232 printests(i, Betas) 

### 1233 1234 } 

### 1235 1236 try(doSummarizeAnalysis(Betas)) 

### 1237

# QFuncti ons.R: Supplemental Functions for Previous Section

## It seems natural to use gls() to do the fitting, but I cannot 

### 1 ## figure out how to get it to handle the covariance structure I need 

### 2 ## for this problem. The model is lP = f(predictions, parameters) + 

### 3 ## epsilon + zeta epsilon is normal(0, sigma^2), and zeta is MVnorm(0, 

### 4 ## defined covariance), where the defined covariance is estimated from 

### 5 ## the MCMC samples. To carry this out, I make a new fitting function 

### 6 ## mvtnormfit(), which takes as arguments: 

### 7 ## y 

### 8 ## x1,...,x3 (not all required) 

### 9 ## M (model number, according to the above scheme) 

### 10 ## lPcov: the defined covariance matrix 

### 11 ## Return value is: 

### 12 ## loglik: log likelihood 

### 13 ## AIC: AIC 

### 14 ## parameters: named appropriately 

### 15 ## R2: defined as 1 - var(y - pred) / var(y) 

### 16 17 ## getparms just returns a named parameter vector for the specified model 

### 18 19 getparms < - function(M) { 

### 20 switch(M+1, 

### 21 ## ------- 0 

### 22 structure(rep(as.numeric(NA),2), names=c("intercept","lsigma")), 

### 23 ## ------- 1 

### 24 s tructure(rep(as.numeric(NA),3), 

### 25 names=c("intercept","NearField","lsigma")), 

### 26 ## ------- 2 

### 27 structure(rep(as.numeric(NA),4), 

### 28 names=c("intercept","NearField","m.lV","lsigma")), 

### 29 ## ------- 3 

### 30 structure(rep(as.numeric(NA),4), 

### 31 names=c("intercept","NearField","m.lR","lsigma")), 

### 32 ## ------- 4 

### 33 structure(rep(as.numeric(NA),5), 

### 34 names=c("intercept","NearField","m.luR","m.lV", 

### 35 "lsigma")), 

### 36 ## ------- 5 

### 37 structure(rep(as.numeric(NA),4), 

### 38 names=c("intercept","NearField","m.lU","lsigma")), 

### 39 ## ------- 6 

### 40 structure(rep(as.numeric(NA),5), 

### 41 name s=c("intercept","NearField","m.luU","m.lV","lsigma")), 

### 42 ## ------- 7 

### 43 structure(rep(as.numeric(NA),6), 

### 44 names=c("intercept","NearField","m.luR","m.luU","m.lV", 

### 45 "lsigma")), 

### 46 ## ------- 8 

### 47 structure(rep(as.numeric(NA),5), 

### 48 names=c("intercept","NearField","m.lR","m.lU","lsigma")), 

### 49 ## ------- 9 

### 50 structure(rep(as.numeric(NA),7), 

### 51 names=c("intercept","NearField", 

### 52 "m.l uR","m.luU","m.lV", 

### 53 "m.luR.N", 

### 54 "lsigma")), 

### 55 ## ------- 10 

### 56 structure(rep(as.numeric(NA),7), 

### 57 names=c("intercept","NearField", 

### 58 "m.luR","m.luU","m.lV", 

### 59 "m.luU.N", 

### 60"lsigma")), 

### 61 ## ------- 11 

### 62 structure(rep(as.numeric(NA),7), 

### 63 names=c("intercept","NearField", 

### 64 "m.luR","m.luU","m.lV", 

### 65 "m.lV.N" , 

### 66 "lsigma")), 

### 67 ## ------- 12 

### 68 structure(rep(as.numeric(NA),9), 

### 69 names=c("intercept","NearField", 

### 70 "m.luR","m.luU","m.lV", 

### 71 "m.luR.N","m.luU.N","m.lV.N", 

### 72 "lsigma")), 

### 73 ## ------- 13 

### 74 structure(rep(as.numeric(NA),8), 

### 75 names=c("intercept","NearField", 

### 76 "m.luR","m.luU","m.lV", 

### 77 "m.luR.N","m.luU.N", 

### 78 "lsigm a")), 

### 79 ## ------- 14 

### 80 structure(rep(as.numeric(NA),8), 

### 81 names=c("intercept","NearField", 

### 82 "m.luR","m.luU","m.lV", 

### 83 "m.luU.N","m.lV.N", 

### 84 "lsigma")), 

### 85 ## -- ----- 15 

### 86 structure(rep(as.numeric(NA),8), 

### 87 names=c("intercept","NearField", 

### 88 "m.luR","m.luU","m.lV", 

### 89 "m.luR.N","m.lV.N", 

### 90 "lsigma")), 

### 91 ## ------- 16 

### 92 structure(rep(as.numeric(NA),6), 

### 93 names=c("intercept","NearField", 

### 94 "m.luR.N","m.luU.N","m.lV.N", 

### 95 "lsigma")), 

### 96 ## ------- 17 

### 97 structure(rep(as.numeric(NA),5), 

### 98 names=c("intercept","NearField", 

### 99 "m.luR.N","m.luU.N", 

### 100 "lsigma")), 

### 101 ## ------- 18 

### 102 structure(rep(as.numeric(NA),5), 

### 103 names=c("intercept","NearField", 

### 104 "m.l uR.N","m.lV.N", 

### 105 "lsigma")), 

### 106 ## ------- 19 

### 107 structure(rep(as.numeric(NA),5), 

### 108 names=c("intercept","NearField", 

### 109 "m.luU.N","m.lV.N", 

### 110 "lsigma")), 

### 111 ## ---- --- 20 

### 112 structure(rep(as.numeric(NA),4), 

### 113 names=c("intercept","NearField", 

### 114 "m.luR.N", 

### 115 "lsigma")), 

### 116 ## ------- 21 

### 117 structure(rep(as.numeric(NA),3), 

### 118 names=c("intercept","m.lR", 

### 119 "lsigma")), 

### 120 ## ------- 22 

### 121 structure(rep(as.numeric(NA),3), 

### 122 names=c("intercept","m.lU", 

### 123"lsigma")) 

### 124 ) 

### 125 } 

### 126 mvtnormfit < - funct ion(df, M, lPcov) { 

### 127 parms < - getparms(M) 

### 128 n < - nrow(df) 

### 129 130 mufun < - with(c(as.list(df),list(n=n)), 

### 131 switch(M+1, 

### 132 ## ------ M0 

### 133 function(p){ 

### 134 rep(p[1],n) 

### 135 }, 

### 136 137 ## ------ M1 

### 138 function(p){ 

### 139 p[1] + p[2] * NearField 

### 140 }, 

### 141 ## ------ M2 

### 142 function(p){ 

### 143 p[1] + p[2] * NearField + p[3] * lV 

### 144 }, 

### 145 ## ------ M3 

### 146 function(p){ 

### 147 p[1] + p[2] * NearField + p[3] * lR 

### 148 }, 

### 149 ## ------ M4 

### 150 function(p){ 

### 151 p[1] + p[2] * NearField + p[3] * ls.R + p[4] * lV 

### 152 }, 

### 153 ## ------ M5 

### 154 function(p){ 

### 155 p[1] + p[2] * NearField + p[3] * lU 

### 156 }, 

### 157 ## ------ M6 

### 158 function(p){ 

### 159 p[1] + p[2] * NearField + p[3] * ls.U + p[4] * lV 

### 160 }, 

### 161 ## ------ M7 

### 162 function(p){ 

### 163 p[1] + p[2] * NearField + p[3] * ls.R + p[4] * ls.U + 

### 164 p[5] * lV 

### 165 }, 

### 166 ## ------ M8 

### 167 function(p){ 

### 168 p[1] + p[2] * NearField + p[3] * lR + p[4] * lU 

### 169 }, 

### 170 ## ------ M9 

### 171 function(p){ 

### 172 p[1] + p[2] * NearField + 

### 173 p[3] * ls.R + 

### 174 p[4] * ls.U + 

### 175 p[5] * lV + 

### 176 p[6] * NearField * ls.R 

### 177 }, 

### 178 ## ------ M10 

### 179 function(p){ 

### 180 p[1] + p[2] * NearField + 

### 181 p[3] * ls.R + 

### 182 p[4] * ls.U + 

### 183 p[5] * lV + 

### 184 p[6] * NearField * ls.U 

### 185 }, 

### 186## ------ M11 

### 187 function(p){ 

### 188 p[1] + p[2] * NearField + 

### 189 p[3] * ls.R + 

### 190 p[4] * ls.U + 

### 191 p[5] * lV + 

### 192 p[6] * NearField * lV 

### 193 }, 

### 194 ## ------ M12 

### 195 function(p){ 

### 196 p[1] + p[2] * NearField + 

### 197 p[3] * ls.R + 

### 198 p[4] * ls.U + 

### 199 p[5] * lV + 

### 200 p[6] * NearField * ls.R + 

### 201 p[7] * NearField * ls.U + 

### 202 p[8] * NearField * lV 

### 203 }, 

### 204 ## ------ M13 

### 205 function(p){ 

### 206 p[1] + p[2] * NearField + 

### 207 p[3] * ls.R + 

### 208 p[4] * ls.U + 

### 209 p[5] * lV + 

### 210 p[6] * NearField * ls.R + 

### 211 p[7] * NearField * ls.U 

### 212 }, 

### 213 ## ------ M14 

### 214 function(p){ 

### 215 p[1] + p[2] * NearField + 

### 216 p[3] * ls.R + 

### 217 p[4] * ls.U + 

### 218 p[5] * lV + 

### 219 p[6] * NearField * ls.U + 

### 220 p[7] * NearField * lV 

### 221 }, 

### 222 ## ------ M15 

### 223 function(p){ 

### 224 p[1] + p[2] * NearField + 

### 225 p[3] * ls.R + 

### 226 p[4] * ls.U + 

### 227 p[5] * lV + 

### 228 p[6] * NearField * ls.R + 

### 229 p[7] * NearField * lV 

### 230 }, 

### 231 ## ------ M16 

### 232 function(p){ 

### 233 p[1] + p[2] * NearField + 

### 234 p[3] * NearField * ls.R + 

### 235 p[4] * NearField * ls.U + 

### 236 p[5] * NearField * lV 

### 237 }, 

### 238 ## ------ M17 

### 239 function(p){ 

### 240 p[1] + p[2] * NearField + 

### 241 p[3] * NearField * ls.R + 

### 242 p[4] * NearField * ls.U 

### 243 }, 

### 244 ## ------ M18 

### 245 function(p){ 

### 246 p[1] + p[2] * NearField + 

### 247 p[3] * NearField * ls.R + 

### 248 p[4] * NearField * lV 

### 249}, 

### 250 ## ------ M19 

### 251 function(p){ 

### 252 p[1] + p[2] * NearField + 

### 253 p[3] * NearField * ls.U + 

### 254 p[4] * NearField * lV 

### 255 }, 

### 256 ## ------ M20 

### 257 function(p){ 

### 258 p[1] + p[2] * NearField + 

### 259 p[3] * NearField * ls.R 

### 260 }, 

### 261 ## ------ M21 

### 262 function(p){ 

### 263 p[1] + p[2] * lR 

### 264 }, 

### 265 ## ------ M22 

### 266 function(p){ 

### 267 p[1] + p[2] * lU 

### 268 }, 

### 269 270 )) 

### 271 of < - function(p) { 

### 272 mu < - mufun(p) 

### 273 dmvnorm(df$lP, mu, lPcov + diag(exp(p[length(p)]*2), nrow(df), nrow(df)), 

### 274 log=TRUE) 

### 275 } 

### 276 277 ## For these models, it works be tter if NearField is a factor 

### 278 df0 < - df 

### 279 df0$NearFieldn < - df0$NearField 

### 280 df0$NearField < - factor(df0$NearField) 

### 281 ## Get initial value of p. Use gls() 

### 282 out < - switch(M+1, 

### 283 ## ------ M0 

### 284 gls(lP ~ 1, 

### 285 data=df0, 

### 286 correlation=corSymm(value=corparms, fixed=TRUE), 

### 287 weights=varFixed(value = ~vars), method="ML"), 

### 288 ## ------ M1 

### 289 gls(lP ~ NearField, 

### 290 data =df0, 

### 291 correlation=corSymm(value=corparms, fixed=TRUE), 

### 292 weights=varFixed(value = ~vars), method="ML"), 

### 293 ## ------ M2 

### 294 gls(lP ~ NearField + lV, 

### 295 data=df0, 

### 296 correlation=corSymm(value=corparms, fixed=TRUE), 

### 297 weights=varFixed(value = ~vars), method="ML"), 

### 298 ## ------ M3 

### 299 gls(lP ~ NearField + lR, 

### 300 data=df0, 

### 301 correlation =corSymm(value=corparms, fixed=TRUE), 

### 302 weights=varFixed(value = ~vars), method="ML"), 

### 303 ## ------ M4 

### 304 gls(lP ~ NearField + ls.R + lV, 

### 305 data=df0, 

### 306 correlation=corSymm(val ue=corparms, fixed=TRUE), 

### 307 weights=varFixed(value = ~vars), method="ML"), 

### 308 ## ------ M5 

### 309 gls(lP ~ NearField + lU, 

### 310 data=df0, 

### 311 correlation=corSymm(value=corparms, fixed=TRUE), 

### 312weights=varFixed(value = ~vars), method="ML"), 

### 313 ## ------ M6 

### 314 gls(lP ~ NearField + ls.U + lV, 

### 315 data=df, 

### 316 correlation=corSymm(value=corparms, fixed=TRUE), 

### 317 weights=varFixed(value = ~vars), method="ML"), 

### 318 ## ------ M7 

### 319 gls(lP ~ NearField + ls.R + ls.U + lV, 

### 320 data=df0, 

### 321 correlation=corSymm(value=corparms, fixed=TRUE), 

### 322 weights=varFixed(value = ~vars), method="ML"), 

### 323 ## ------ M8 

### 324 gls(lP ~ NearField + lR + lU, 

### 325 data=df0, 

### 326 correlation= corSymm(value=corparms, fixed=TRUE), 

### 327 weights=varFixed(value = ~vars), method="ML"), 

### 328 ## ------ M9 

### 329 gls(lP ~ NearField * ls.R + ls.U + lV, 

### 330 data=df0, 

### 331 correlation=corSy mm(value=corparms, fixed=TRUE), 

### 332 weights=varFixed(value = ~vars), method="ML"), 

### 333 ## ------ M10 

### 334 gls(lP ~ NearField + ls.R + NearField * ls.U + lV, 

### 335 data=df0, 

### 336 correlati on=corSymm(value=corparms, fixed=TRUE), 

### 337 weights=varFixed(value = ~vars), method="ML"), 

### 338 ## ------ M11 

### 339 gls(lP ~ NearField + ls.R + ls.U + NearField * lV, 

### 340 data=df0, 

### 341 correlation=corSymm(value=corparms, fixed=TRUE), 

### 342 weights=varFixed(value = ~vars), method="ML"), 

### 343 ## ------ M12 

### 344 gls(lP ~ NearField * (ls.R + ls.U + lV), 

### 345 data=df0, 

### 346 correlation=corSymm(value=corparms, fixed=TRUE), 

### 347 weights=varFixed(value = ~vars), method="ML"), 

### 348 ## ------ M13 

### 349 gls(lP ~ NearField + NearField * (ls.R + ls.U) + lV, 

### 350 dat a=df0, 

### 351 correlation=corSymm(value=corparms, fixed=TRUE), 

### 352 weights=varFixed(value = ~vars), method="ML"), 

### 353 ## ------ M14 

### 354 gls(lP ~ NearField + ls.R + NearField * (ls.U + lV), 

### 355 data=df0, 

### 356 correlation=corSymm(value=corparms, fixed=TRUE), 

### 357 weights=varFixed(value = ~vars), method="ML"), 

### 358 ## ------ M15 

### 359 gls(lP ~ NearField + NearField * ls.R + ls.U + NearField * lV, 

### 360 data=df0, 

### 361 correlation=corSymm(value=corparms, fixed=TRUE), 

### 362 weights=varFixed(value = ~vars), method="ML"), 

### 363 ## ------ M16 

### 364 gls(lP ~ NearField + I(NearFieldn*ls .R) + 

### 365 I(NearFieldn*ls.U) + I(NearFieldn*lV), 

### 366 data=df0, 

### 367 correlation=corSymm(value=corparms, fixed=TRUE), 

### 368 weights=varFixed(value = ~vars), method="ML"), 

### 369 ## ------ M17 

### 370 gls(lP ~ NearField + I(NearFieldn*ls.R) + I(NearFieldn*ls.U), 

### 371 data=df0, 

### 372 correlation=corSymm(value=corparms, fixed=TRUE), 

### 373 weights=varFixed(value = ~va rs), method="ML"), 

### 374 ## ------ M18 

### 375gls(lP ~ NearField + I(NearFieldn*ls.R) + I(NearFieldn*lV), 

### 376 data=df0, 

### 377 correlation=corSymm(value=corparms, fixed=TRUE), 

### 378 weights=var Fixed(value = ~vars), method="ML"), 

### 379 ## ------ M19 

### 380 gls(lP ~ NearField + I(NearFieldn*ls.U) + I(NearFieldn*lV), 

### 381 data=df0, 

### 382 correlation=corSymm(value=corparms, fixed=TRUE), 

### 383 weights=varFixed(value = ~vars), method="ML"), 

### 384 ## ------ M20 

### 385 gls(lP ~ NearField + I(NearFieldn*ls.R), 

### 386 data=df0, 

### 387 correlation=corSymm(value=corparms, fixed=TRUE), 

### 388 weights=varFixed(value = ~vars), method="ML"), 

### 389 ## ------ M21 

### 390 gls(lP ~ lR, 

### 391 data=df0, 

### 392 correlation=corSymm(value=corparms, fixed=TRUE), 

### 393 weights=varFixed(value = ~ vars), method="ML"), 

### 394 ## ------ M22 

### 395 gls(lP ~ lU, 

### 396 data=df0, 

### 397 correlation=corSymm(value=corparms, fixed=TRUE), 

### 398 weights=varFixed(value = ~vars), method="ML"), 

### 399 ) 

### 400 401 SS < - median((out$sigma^2 - 1)*diag(lPcov)) 

### 402 pstart < - c(coef(out),log(SS)/2) 

### 403 out2 < - try(optim(pstart, of, method="BFGS",control=list(fnscale= - 1))) 

### 404 ## If try - error, first try using simplex method 

### 405 if (is(out2, "try - error")) { 

### 406 out2 < - try(optim(pstart, of, control=list(fnscale= - 1, maxit=2000))) 

### 407 } 

### 408 if (!is(out2, "try - error")) { 

### 409 muhat < - with(df, mufun(out2$par)) 

### 410 par < - structure(out2$par, names=names(parms)) 

### 411 list(LogLik = out2$value, 

### 412 ## AIC = - 2* (out2$value) - 2*length(pstart), WRONG! Corrected 7/5/2012 

### 413 ## RWS 

### 414 AIC = - 2*(out2$value) + 2*length(pstart), 

### 415 R2 = 1 - var(df$lP - muhat)/var(df$lP), 

### 416 coefficients = par) 

### 417 } else { 

### 418 out2 

### 419 } 

### 420 } 

### 421 422 423 ## Distributions of AIC values across replicate samples 

### 424 mysummary < - function(x){ 

### 425 qq < - stats::quantile(x, pr=c(0,0.025,0.25,0.5,0.75,0.975,1), na.rm=TRUE) 

### 426 c(qq[1:4],mean=mean(x, na.rm=TRUE),qq[5:7]) 

### 427 } 

### 428 429 ## Changes 7/5/2012: separator for dll is now "< - >" RWS 

### 430 ## Also, com parison between delta LL and null delta LL was wrong. RWS 

### 431 mycomps < - function(z, comp=c(0,1),NEW = TRUE) { 

### 432 dll < - if (NEW) { 

### 433 paste("dLL",comp[1],"< - >",comp[2],sep="") 

### 434 } else { 

### 435 paste("dLL",comp[1],comp[2],sep="") 

### 436 } 

### 437 LLsimp < - paste("LL",comp[ 1],sep="_") 

### 438LLcompl < - paste("LL",comp[2],sep="_") 

### 439 R2simp < - paste("R2",comp[1],sep="_") 

### 440 R2compl < - paste("R2",comp[2],sep="_") 

### 441 c( 

### 442 sum(as.numeric(abs(z[,LLcompl] - z[,LLsimp]) < abs(z[,dll])))/nrow(z), 

### 443 if(comp[1] > 0) { 

### 444 mean(z[,R2comp l] - z[,R2simp], na.rm=TRUE) 

### 445 } else { 

### 446 mean(z[,R2compl]) 

### 447 } 

### 448 ) 

### 449 } 

### 450 451 printests < - function(n, Betas) { 

### 452 pat < - paste("_",n,"$",sep="") 

### 453 indx < - grep(pat, colnames(Betas)) 

### 454 mns < - apply(Betas[,indx],2,mean) 

### 455 lcl < - apply(Betas[,indx],2,quantile,pr=0.025) 

### 456 ucl < - apply(Betas[,indx],2,quantile,pr=0.975) 

### 457 458 out < - cbind(mns, lcl, ucl) 

### 459 rnames < - names(mns) 

### 460 rnames[1] < - "Intercept" 

### 461 rownames(out) < - rnames 

### 462 colnames(out) < - c("Mean","LCL","UCL") 

### 463 writeLines(c("" ,paste("Bootstrap Distributions: Model",n))) 

### 464 print(out) 

### 465 } 

### 466

# Fitmodel16.R: R code to fit Model 16 to NHANES Parent Chemical Data

library(doMC) 

### 1 registerDoMC(cores=3) 

### 2 library(rjags) 

### 3 load.module("lecuyer") 

### 4 library(random) 

### 5 library(R2jags) 

### 6 7 ## JAGS mod el is in model16.jag 

### 8 ## It needs the following data: 

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

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

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

### 12 ## N = number of parents, 

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

### 14 ## in R.x) 

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

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

### 17 ## i contributes. Corresponding values o f Phi[i,] are non - zero. 

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

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

### 20 ## M = number of metabolites, 

### 21 ## tau[1:M] : precision for observed metabolite r ates. 

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

### 23 ## metabolite rates. 

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

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

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

### 27 ## 

### 28 ## The following variables need initializing 

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

### 30 ## a 

### 31 ## m.R 

### 32 ## m.U 

### 33 ## X1 

### 34 ## X2 

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

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

### 38 39 attach(file.path("..","MCMCdata.RData")) 

### 40 z < - runif(1) 

### 41 doONEchain < - TRUE 

### 42 if (doONEchain) { ## this is a restart; use the chains we already have 

### 43 ## to initialize stuff. 

### 44 45 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

### 46 ## JAGS inputs: nhanesdata and inits 

### 47 ## ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 

### 48 49 nhanesdata < - list() 

### 50 ## The Jags models need these data: 

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

### 52 nhanesdata$IURcuts < - log(cutpoints) 

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

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

### 55 ## IURvols - indicator for which inter val an IUR production volume falls 

### 56 nhanesdata$IURvols < - 

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

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

### 59 nhanesdata$CPRIindx < - which(pred.data$PUsource == "CPR I 2002") 

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

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

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

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

### 64 nhanesdata$NearField < - pred.d ata$NearField 

### 65 66 ## ly - the log of the observed metabolite quantiles 

### 67 nhanesdata$ly < - 

### 68 ly < - 

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

### 70 ## Q, the proportions associated with the quantiles 

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

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

### 73 quantiles 

### 74 nhanesdata$se < - Measured$se 

### 75 ## llod - log(Limit of Detection) 

### 76 nhanesdata$llod < - log(Measured$LOQ) 

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

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

### 79 ## Some counts: 

### 80 ## N is the number of parents 

### 81 nhanesdata$N < - nrow(pred.data) 

### 82 ## M is the number of metabolites 

### 83 nhanesdata$M < - nrow(Measured) 

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

### 85 nhanesdata$N.IUR < - length(nhanesdata$IURindx) 

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

### 87 nhanesdata$N.CPRI < - length(nhanesdata$CPRIindx) 

### 88 ## N.nontriv is the number of parents that result in >1 me tabolite 

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

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

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

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

### 93 ## of the ith parent 

### 94 nhanesdata$indx < - 

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

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

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

### 98 tmp < - cumsum(nhanesdata$Nmets) 

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

### 100 nhanesdata$dstop < - tmp 

### 101 102 ## Phi 

### 103 nhanesdata$Phi < - Xref 

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

### 105 106 ## --------------------------------------------------------------------- 

### 107 ## Initialization 

### 108 109 inits < - list(list()) 

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

### 111 inits[[1]]$sdlVmu.dum < - sd(log(pred.data$V)) 

### 112 inits[[1]]$lVmu < - mean(log(pred.data$V)) 

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

### 114 inits[[1]]$lP < - rep(mean(nhanesdata$ly,na.rm=TRUE), nrow(pre d.data)) 

### 115 116 inits[[1]]$a.F < - - 7 ## mean(nhanesdata$ly, na.rm=TRUE) 

### 117 inits[[1]]$a.N < - inits[[1]]$a.F 

### 118 inits[[1]]$m.R.N < - - 1.2 

### 119 inits[[1]]$m.U.N < - 0.60 

### 120 inits[[1]]$m.V.N < - 0.20 

### 121 inits[[1]]$sd.dum < - 3 

### 122 inits[[1]]$sd < - rep(1.0, ncol(ly)) 

### 123## for the censored val ues in ly 

### 124 ily < - is.na(ly) 

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

### 126 lyinit < - ifelse(ily,loqmx - 2,NA) 

### 127 inits[[1]]$ly < - lyinit 

### 128 inits[[1]]$lsd < - 1 

### 129 inits[[1]]$lsd.dum < - 2 

### 130 131 inits[[1]]$delta < - rep(1, nhanesdata$Ndelta) 

### 132 133 134 doinits < - function(chain) {inits[[chain]]} 

### 135 136 assign("inits", inits, envir=environment(doinits)) 

### 137 138 ## JAGS code 

### 139 140 ## model < - jags.model("model3.jag",data=nhanesdata,inits=inits[[1]], 

### 141 ## n.chains=1,n.adapt=1000) 

### 142 143 144 ## Now, run for 300000 samples, thinned every 100, and dump lP, lU, and delta 

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

### 146 147 out.samples < - jags(data=nhanesdata, inits=doinits, 

### 148 parameters.to.save=c("a.F","a.N", 

### 149 "m.R.N","m.U.N","m.V.N", 

### 150 "lV","sd.V","lP","phi","lU","sd.lVmu","lVmu","Z", 

### 151 "sd","lsd","lsd.sd"), 

### 152 model.file="model16.jag", 

### 153 n.chains=1, n.iter=3000, n.burnin=1000, n.thin=2, 

### 154 DIC=FALSE,progress.bar="text") 

### 155 156 157 out.coda < - as.mcmc(out.samples) 

### 158 nms < - varnames(out.coda) 

### 159 for (nm in nms) { 

### 160 writeLines(c(""," ----------------- ",nm,"")) 

### 161 print(try(heidel.diag(out.coda[,nm]))) 

### 162 print(try(summary(out.coda[,nm]))) 

### 163 } 

### 164 165 save(out.samples, out.coda, nhanesdata, inits, file="m3parms1.RData") 

### 166 } 

### 167 168 load("m3parms1.RData") 

### 169 170 ## -------------------------------------------------------------------------------- 

### 171 ## This is the initial, hopefully converged, chain. Now, use these to sta rt with 

### 172 ## and run three new chains. 

### 173 174 175 ## getInits finds extreme values 

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

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

### 178 dta < - scale(X) 

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

### 180 outinits < - numeric(ninits) 

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

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

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

### 184 for (i in 1:ninits) { 

### 185 u < - rnorm(ncol(dta)) 

### 186u < - u/sqrt(sum(u*u)) 

### 187 lngt < - dta %*% u 

### 188 indx < - seq(along=lngt) 

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

### 190 oind < - order(lngt) 

### 191 outinits[i] < - oind[keep] 

### 192 } 

### 193 outinits 

### 194 } 

### 195 ## Get 3 hyperdispersed points. 

### 196 Indx < - getInits(out.coda,3) 

### 197 198 199 200 getsamps < - function(nm) { 

### 201 nm < - gsub(" \ \ ."," \ \ \ \ .",nm) 

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

### 203 out.coda[,cindx,drop=FALSE] 

### 204 } 

### 205 206 inits2 < - list(list(),list(),list()) 

### 207 for (i in 1:3) { 

### 208 inits2[[i]]$l V.IUR < - getsamps("^lV \ \ [")[Indx[i],nhanesdata$IURindx] 

### 209 inits2[[i]]$sdlVmu.dum < - getsamps("^sd.lVmu")[Indx[i],] 

### 210 inits2[[i]]$lVmu < - getsamps("^lVmu")[Indx[i],] 

### 211 inits2[[i]]$lV.CPRI.mu < - getsamps("^lV \ \ [")[Indx[i],nhanesdata$CPRIindx] 

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

### 213 ##inits2[[i]]$lP < - inits[[1]]$lP 

### 214 inits2[[i]]$a.F < - getsamps("^a.F")[Indx[i],] 

### 215 ##inits2[[i]]$a.F < - inits[[1]]$a.F 

### 216 inits2[[i]]$a.N < - getsamps("^a.N")[Indx[i],] 

### 217 inits2[[i]]$m.R.N < - getsamps("^m.R .N")[Indx[i],] 

### 218 inits2[[i]]$m.U.N < - getsamps("^m.U.N")[Indx[i],] 

### 219 inits2[[i]]$m.V.N < - getsamps("^m.V.N")[Indx[i],] 

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

### 221 inits2[[i]]$sd < - getsamps("^sd \ \ [")[Indx[i],] 

### 222 inits2[[i]]$lsd < - getsamps("^lsd$ ")[Indx[i],] 

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

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

### 225 inits2[[i]]$ly< - inits[[1]]$ly 

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

### 227 inits2[[i]]$.RNG.seed=randomNumbers(n=1, min=1, max=1e+06, co l=1) 

### 228 } 

### 229 230 ## quickie tool to compare inits and inits2 visually 

### 231 ## cmpin < - function(nm) { 

### 232 ## cbind(inits[[1]][[nm]], inits2[[1]][[nm]], 

### 233 ## inits2[[2]][[nm]], inits2[[3]][[nm]]) 

### 234 ## } 

### 235 236 ## jags.parallel doesn't export all the global environment, so in its2 

### 237 ## gets lost, and doinits cannot find it in their subprocesses. 

### 238 ## The workaround used here is to create an environment for the 

### 239 ## function doinits and explicitly assign inits2 to it. Then, when 

### 240 ## doinits gets exported, inits2 goes with it. 

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

### 242 env1 < - new.env() 

### 243 assign("inits2", inits2, pos=env1) 

### 244 environment(doinits) < - env1 

### 245 246 NIters < - 1000 

### 247 Thin < - 500 

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

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

### 251 n.adapt=100000) 

### 252 result < - coda.samples(model, variable.names=c("a.F","a.N", 

### 253 "m.R.N","m.U.N","m.V.N", 

### 254 "lV","sd.V","lP","lU","sd.lVmu","lVmu","Z", 

### 255 "sd","lsd","lsd.sd"), 

### 256 n.iter=NIters*Thin, thin=Thin) 

### 257 return(result) 

### 258 } 

### 259 260 261 262 out.coda3R < - do.call("mcmc.list", lapply(out.samps3R, function(z) z [[1]])) 

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

### 264 nms < - varnames(out.coda3R) 

### 265 for (nm in nms) { 

### 266 writeLines(c(""," ----------------- ",nm,"")) 

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

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

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

### 270 } 

### 271 272 273 274 275 276

# model16.jag: JAGS Code for Fitting Model 16 to NHANES Parent

Chemical Data model 

### 1 { 

### 2 ## Production volumes are crudely estimated or r eported. We have two 

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

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

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

### 6 ## of us e/production as lognormal (since everything here is on a log 

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

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

### 9 ## while each CPRI estimate is sampled from a lognorma l distribution 

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

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

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

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

### 14 ## predictions 

### 15 ## ++++++++++++++++++++++++++++++++++++++++++++ 

### 16 ## Model the volumes 

### 17 ## ++++++++++++++++++++++++++++++++++++++++++++ 

### 18 ## ------------------------------------------- 

### 19 ## IUR volumes 

### 20 for (i in 1:N.IUR) { 

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

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

### 23 lV[IURindx[i]] < - lV.IUR[i] 

### 24 } 

### 25 lVmu ~ dnorm(0, 0.001) 

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

### 27 sd.lVmu < - abs(sdlVmu.dum) 

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

### 29 ## -------------- ----------------------------- 

### 30 ## CPRI usage 

### 31 for (i in 1:N.CPRI) { 

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

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

### 34 lV[CPRIindx[i]] < - lV.CPRI.mu[i] 

### 35 } 

### 36 ## +++++++++++++++++++++++++ ++++++++++++++++++++ 

### 37 ## The heart of the model: 

### 38 ## Model the parent exposures. Fit separate 

### 39 ## coefficients for far - and nearfield chemicals 

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

### 41 for (i in 1:N) { 

### 42 Z[i] < - 

### 43 (1 - NearField[i]) * a.F + 

### 44 NearField[i]*(a.N + m.R.N * ls.R[i] + 

### 45 m.U.N * ls.U[i] + 

### 46 m.V.N * lV[i]) 

### 47 lP[i] ~ dnorm(Z[i], tau.V) 

### 48 P[i]< - exp(lP[i]) 

### 49 } 

### 50 ## Priors for a, m.R, m.U, and m.V 

### 51 a.F ~ dno rm(0, 0.001) 

### 52 a.N ~ dnorm(0, 0.001) 

### 53 m.R.N ~ dnorm(0, 0.001) 

### 54 m.U.N ~ dnorm(0, 0.001) 

### 55 m.V.N ~ dnorm(0, 0.001) 

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

### 57 ## Fit a separate standard deviation for far and near 

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

### 59 sd.V < - abs(sd.dum) 

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

### 61 ## ++++++++++++++++++++++++++++++++++++++++++++++ 

### 62 ## Link the unobserved parent exposure to the observed 

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

### 64 ## values. 

### 65 ## 

### 66 ## Model metabolit e concentrations as (usually) unknown fractions 

### 67 ## of the parent. 

### 68 ## 

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

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

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

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

### 73 for (j in 1:M) { 

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

### 75 for (i in 1:4) { 

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

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

### 78 } 

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

### 80 ## lognormal distribution for chemical - specific population sds 

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

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

### 83 } 

### 84 lsd ~ dnorm(0, 0.001) 

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

### 86 lsd.sd < - abs(lsd.dum) 

### 87 tau.lsd < - 1/(lsd.sd * lsd.sd) 

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

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

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

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

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

### 93 ## metabolite matched to parent i. 

### 94 ## 

### 95 ## The elements of ph i are essentially unknown. 

### 96 for (k in 1:Ndelta) { 

### 97 delta[k] ~ dunif(0,1) 

### 98 } 

### 99 for (i in 1:N.nontriv) { 

### 100 for (j in 1:Nmets[i]) { 

### 101 phi[dstart[i] + j - 1] < - delta[dstart[i] + j - 1] / 

### 102 sum(delta[dstart[i]:dstop[i]]) 

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

### 104 } 

### 105 } 

### 106 } 

### 107 108