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]]

