in addition, I get this error whenever I want to plot the predicted values
plot(augPred(mod2))Error in sprintf(gettext(fmt, domain = domain), ...) :
invalid type of argument[1]: 'symbol'
From: ommo1@hotmail.com
To: r-help@r-project.org
Subject: [R] multivariate, hierarchical model
Date: Thu, 2 May 2013 10:26:57 +0000
Sorry for the last email, sent too early.
I have a small data set that has a hierarchical structure. It has both temporal
(year, months) and spatial (treatment code and zone code). The following
explains the data:
WSZ_Code the
water supply zone code (1 to 8)
Treatment_Code the
treatment plant which supplies each water supply zone (1 to 4)
Year year
of sampling (1996 - 2000)
Month month
of the year of sampling, 1=January
TTHM the
total trihalomethane concentration (ìg/L) in the tap water sample
CL2_Free concentration
of free chlorine (mg/L) in the water sample - indicates the level of the
disinfecting chlorine dose which has not reacted with organic matter in the
water.
BrO3 concentration
of bromate (?g/L) in the water sample. BrO3 is formed during certain types of
water treatment. This variable contains some missing values due to information
not being recorded for the sample - these are denoted by "NA".
Colour measure
of colour of the water sample - this is one possible indicator of the level of
organic matter in the water. Units = Hazen.
pH pH
of the water sample.
Turbidity measure of "cloudiness" of
the water sample - caused by particles suspended in the water. Units = FTU
The aim of the analysis is to produce exposure estimates for TTHM. I has split
the months into seasons as there appears to be a seasonality trend. I've
also carried out multiple imputation for any missing values. I've done a
correlation analysis on all the variables and observed which ones are well
correlated with TTHM. It appears that CL2_Free is the most significantly
correlated so i've decided to include that in the final model. Positively
skewed
I'm thinking of doing a mixed effects model with random intercepts as the
treatment code and zones within the treatment cose and random slopes as the
seasons.
mod2 <- lme(tthm ~ cl2free, random= ~ seasons| treatcode/loc_code)
but that doesn't work.
these seems to work good:
mod2 <- lme(tthm ~ cl2free, random= ~ 1| loc_code, data=new.data,
method="ML")
mod3 <- lme(tthm ~ cl2free, random= ~ 1| treatcode/loc_code, data=new.data,
method ="ML")
mod2 has a lower AIC, so it appears better.
Should I group (merge) treatment code and location code?
Also, cl2free has a very positively skewed distribution, should I transform it?
I want to incorporate seasons, possibly as a random slope, but R doesn't
seem to like it.
Help greatly appreciated (and sorry for long email).
Omnia
[[alternative HTML version deleted]]