Sunday, December 1, 2013

In my previous post I ignored the fact that some data are below the detection level. I would not know how to handle those in a mixed model from lme4 or nlme. However, JAGS can handle these values. Next to that I kept the usual independent variables, such as random effects for location, machine etc.

Data

Data has been acquired and processed similar to last week. For script at the bottom of the page.

detection level

In the spreadsheet a '4' is used to indicate a value below detection limit. (* een " 4 "-teken in kolom C geeft aan dat de analyseresultaten beneden de detectielimiet liggen). There is something strange there. The detection limit seems to be 0.4, but there is data under 0.4 which is not flagged. In addition, some of the data is below 0, which seems improbable. In the end I flagged all data under 0.4 as NA.

Fe <- rf2[rf2$Compound==levels(rf2$Compound)[9],]

library(lattice)

densityplot( ~ Amount | machine + Status,

groups=Location,data=Fe[!is.na(Fe$Status)

& Fe$Amount<.8,],

panel=function(...) {

panel.densityplot(...)

panel.abline(v=c(0,0.4))

},adjust=.5

)

Analysis

It is a reasonable standard JAGS model. The only tricky thing is initialization. When the parameters are too far away an error is thrown. Hence 'intercept' is initialized as -1. In text the model and inits are given, at the end the full code.

FeModel <- function() {

for (i in 1:nobs) {

BLZ[i] ~ dinterval( LAmount[i] , -0.4034029 )

ExpLAmount[i] <- intercept +

rate*time[i] +

rateL[Location[i]]*time[i]+

FMach[Machine[i]] +

FLoc[Location[i]]

LAmount[i] ~ dnorm(ExpLAmount[i],tau[Machine[i]])

}

for (i in 1:2) {

tau[i] <- pow(sd[i],-2)

sd[i] ~ dunif(0,1000)

}

intercept ~ dnorm(0,.0001)

FMach[1] ~dnorm(0,tauMachShift)

FMach[2] <- -FMach[1]

tauMachShift <- pow(sdMachShift,-2)

sdMachShift ~ dunif(0,10)

for (iL in 1:nloc) {

FLoc[iL] ~ dnorm(0,tauLoc)

rateL[iL] ~dnorm(0,tauRate)

yearlyL[iL] <- pow(pow(10,rateL[iL]+rate),365)

}

tauLoc <- pow(sdLoc,-2)

sdLoc ~ dunif(0,100)

tauRate <- pow(sdRate,-2)

sdRate ~ dunif(0,100)

rate ~ dnorm(0,.0001)

dailydec <- pow(10,rate)

yearlydec <- pow(dailydec,365)

}

inits <- function() {

list(intercept = -1,#rnorm(1,.45,.04),

rate=rnorm(0,.01),

sd=runif(2,0,50),

sdMachShift=runif(1,0,1),

sdLoc=runif(1,0,1))

}

Results

The model shows the decrease of 5% per year which was expected.

Machine 1 (Eigenbrodt NSA 181 KHT) is slightly less accurate (was also observed last week) and has slightly lower values than 2 (Van Essen/ECN vanger). Variation between locations is less than variation between machines.

Inference for Bugs model at "C:/.../AppData/Local/Temp/Rtmpma5UEx/modelc40761211ba.txt", fit using jags,

No comments:

Post a Comment

Wiekvoet

Wiekvoet is about R, JAGS, STAN, and any data I have interest in. Topics range from sensometrics, statistics, chemometrics and biostatistics. For comments or suggestions please email me at wiekvoet at xs4all dot nl.