**Wiekvoet**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

After reading Blog About Stats’ Open Data Index Blog Post I decided to browse a bit in the Open Data Index. Choosing Netherlands and following Emission of Pollutants I ended on a page from National Institute for Public Health. The page http://www.lml.rivm.nl/data/gevalideerd/index.html is in Dutch and chose composition of rain water 1992-2005.There is also 2006-2008, 2009, 2010 and 2011. In 2006 machines have changed, ‘a report describing changes is in preparation’. I am only using 1992-2005 in this post.

### Data

The data is in a spreadsheet, one page per location. There is quite a big header on each page and each data column has an associated column marking data below detection limit. This being real world data there are also the odd missing data.

Following Read Excel Data from R I chose to read the data via XLConnect. The whole process is wrapped in a function, So I can read all pages (all locations).

library(XLConnect)

wb = loadWorkbook(“lmre_1992-2005.xls”)

sheets <- getSheets(wb)[-1]

readsheet <- function(cursheet) {

df = readWorksheet(wb, sheet = cursheet , header = TRUE)

topline <- 8

test <- which(df[6,]==’C’)+1

numin <- df[topline:nrow(df),test]

names(numin) <- make.names( gsub(‘ ‘,”,df[6,test]))

for(i in 1:ncol(numin)) numin[,i] <- as.numeric(gsub(‘,’,’.’,numin[,i]))

status = df[topline:nrow(df),test-1]

names(status) <- paste(‘C’,make.names( gsub(‘ ‘,”,df[6,test])),sep=”)

numin$StartDate <- as.Date(substr(df[topline:nrow(df),1],1,10))

numin$EndDate <- as.Date(substr(df[topline:nrow(df),2],1,10))

numin <- cbind(numin,status)

tall <- reshape(numin,

varying=list(make.names( gsub(‘ ‘,”,df[6,test])),

paste(‘C’,make.names( gsub(‘ ‘,”,df[6,test])),sep=”)),

v.names=c(‘Amount’,’Status’),

timevar=’Compound’,

idvar=c(‘StartDate’,’EndDate’),

times=paste(df[6,test],'(‘,df[5,test],’)’,sep=”),

direction=’long’)

tall$Compound <- factor(gsub(‘ )’,’)’,gsub(‘ +’,’ ‘,tall$Compound)))

row.names(tall)<- NULL

tall$Location <- cursheet

tall

}

tt <- lapply(sheets,readsheet)

tt2 <- do.call(rbind,tt)

tt2$Location <- factor(tt2$Location)

### Plots

[4] “Cl (umol/l)” “Co (umol/l)” “Cr (umol/l)”

[7] “Cu (umol/l)” “F (umol/l)” “Fe (umol/l)”

[10] “K (umol/l)” “K25 (uS/cm)” “massa_hc (g)”

[16] “Na (umol/l)” “NH4 (umol/l)” “Ni (umol/l)”

[19] “NO3 (umol/l)” “nsl (mm)” “Pb (umol/l)”

[22] “pH (pH-eenheid)” “PO4 (umol/l)” “SO4 (umol/l)”

#### pH

#### Iron

### Mixed model for Iron

Just to get some quantification on the reduction I ran a model. The obvious choice is again to use logarithmic transformed data. However, some zeros unleashed havoc on the subsequent estimates. Hence these data are made NA prior to computations.

library(lme4)

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

# plot(density(Fe$Amount[!is.na(Fe$Status) & Fe$Status!=0],adjust=.5))

Fe$LAmount <- log10(Fe$Amount)

Fe$LAmount[Fe$Amount<=0] <- NA

ll <- lmer(LAmount ~ StartDate + (1|Location) +(StartDate|Location ) ,

data=Fe)

summary(ll)

Linear mixed model fit by REML [‘lmerMod’]

Formula: LAmount ~ StartDate + (1 | Location) +

(StartDate | Location)

Data: Fe

REML criterion at convergence: 2170.498

Random effects:

Groups Name Variance Std.Dev. Corr

Location (Intercept) 4.727e-11 6.875e-06

Location.1 (Intercept) 1.456e-04 1.207e-02

StartDate 1.348e-10 1.161e-05 1.00

Residual 1.436e-01 3.789e-01

Number of obs: 2338, groups: Location, 17

Fixed effects:

Estimate Std. Error t value

(Intercept) 5.001e-01 5.971e-02 8.376

StartDate -6.975e-05 6.338e-06 -11.005

Correlation of Fixed Effects:

(Intr)

StartDate -0.862

**leave a comment**for the author, please follow the link and comment on their blog:

**Wiekvoet**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.