#install.packages("nlme")
library("nlme")
library(lattice)
#install.packages("lme4")
#library("lme4")

RIKZ = read.table("http://www.math.pku.edu.cn/teachers/xirb/Courses/biostatistics/Biostatistics2014/RIKZ.txt",header=T)


Beta<-vector(length=9)
for (i in 1:9){
 tmpout<-summary(lm(Richness~NAP,subset = (Beach==i),data=RIKZ))
 Beta[i]<-tmpout$coefficients[2,1]
}

aggregate(1:nrow(RIKZ),by=list(RIKZ$Beach),FUN=function(v){RIKZ$Exposure[v[1]]})

###put 8 and 10 as "a", 11 as "b"
ExposureBeach <- c("a","a","b","b","a",
                     "b","b","a","a")
tmp2 <- lm(Beta ~ factor(ExposureBeach),data=RIKZ)


library(nlme)
Mlme1 <- lme(Richness ~ NAP, random = ~1 | Beach,data=RIKZ)
summary(Mlme1)
RIKZ$fBeach <- factor(RIKZ$Beach)
Mlme1 <- lme(Richness ~ NAP, random = ~1 | fBeach,data=RIKZ)
summary(Mlme1)

xyplot(fitted(Mlme1)~NAP|fBeach,data=RIKZ,type="l")

F0<-fitted(Mlme1,level=0)
F1<-fitted(Mlme1,level=1)
I<-order(RIKZ$NAP)
NAPs<-sort(RIKZ$NAP)
plot(NAPs,F0[I],lwd=4,type="l",ylim=c(0,22),
     ylab="Richness",xlab="NAP")
for (i in 1:9){
   x1<-RIKZ$NAP[RIKZ$Beach==i]
   y1<-F1[RIKZ$Beach==i]
   K<-order(x1)
   lines(sort(x1),y1[K])
}
text(RIKZ$NAP,RIKZ$Richness,RIKZ$Beach,cex=0.9)

###random slope but fixed effect


Mlme1.1 <- lme(Richness ~ NAP,
             random = ~ 0 + NAP | fBeach, data = RIKZ)
summary(Mlme1.1)

F0<-fitted(Mlme1.1,level=0)
F1<-fitted(Mlme1.1,level=1)
I<-order(RIKZ$NAP)
NAPs<-sort(RIKZ$NAP)
plot(NAPs,F0[I],lwd=4,type="l",ylim=c(0,22),
     ylab="Richness",xlab="NAP")
for (i in 1:9){
   x1<-RIKZ$NAP[RIKZ$Beach==i]
   y1<-F1[RIKZ$Beach==i]
   K<-order(x1)
   lines(sort(x1),y1[K])
}
text(RIKZ$NAP,RIKZ$Richness,RIKZ$Beach,cex=0.9)

xyplot(fitted(Mlme1.1)~NAP|fBeach,data=RIKZ,type="l")

###random slope and intercept
Mlme2 <- lme(Richness ~ NAP,
             random = ~1 + NAP | fBeach, data = RIKZ)
summary(Mlme2)
## or equavalently 
Mlme2 <- lme(Richness ~ NAP,
             random = ~NAP | fBeach, data = RIKZ)
summary(Mlme2)

F0<-fitted(Mlme2,level=0)
F1<-fitted(Mlme2,level=1)
I<-order(RIKZ$NAP)
NAPs<-sort(RIKZ$NAP)
plot(NAPs,F0[I],lwd=4,type="l",ylim=c(0,22),
     ylab="Richness",xlab="NAP")
for (i in 1:9){
   x1<-RIKZ$NAP[RIKZ$Beach==i]
   y1<-F1[RIKZ$Beach==i]
   K<-order(x1)
   lines(sort(x1),y1[K])
}
text(RIKZ$NAP,RIKZ$Richness,RIKZ$Beach,cex=0.9)


xyplot(fitted(Mlme2)~NAP|fBeach,data=RIKZ,type="l")




#### without the NAP effect
Mlme3 <- lme(Richness ~ 1, random = ~1 | fBeach,
           data = RIKZ)
summary(Mlme3)
xyplot(fitted(Mlme3)~NAP|fBeach,data=RIKZ,type="l")



###generalized least square with correlated data
M.mixed <- lme(Richness ~ NAP, random = ~1 | fBeach,
                 method = "REML", data = RIKZ)
M.gls <- gls(Richness ~ NAP, method = "REML",
          correlation = corCompSymm(form =~1 | fBeach),
          data = RIKZ)




###compare different methods and models
RIKZ$fExp<-RIKZ$Exposure
RIKZ$fExp[RIKZ$fExp==8]<-10
RIKZ$fExp<-factor(RIKZ$fExp,levels=c(10,11))
M0.ML <- lme(Richness ~ NAP, data = RIKZ,
              random = ~1| fBeach, method = "ML")
M0.REML <-lme(Richness ~ NAP, random = ~1|fBeach,
              method = "REML", data = RIKZ)
M1.ML <- lme(Richness ~ NAP+fExp, data = RIKZ,
              random = ~1| fBeach, method = "ML")
M1.REML <- lme(Richness ~NAP+fExp, data = RIKZ,
              random = ~1|fBeach, method = "REML")



#####################################
#####nested and crossed model########
#####################################



library("lme4")
library("nlme")
data("Oxide")
Oxide <- as.data.frame(Oxide)
###this is from http://www.r-bloggers.com/nested-vs-non-nested-crossed-random-effects-in-r/

## In the Oxide data, Site is nested in Wafer, which
## is nested in Lot. But, the latter appear crossed (alghough not really crossed):
#### lot -> Wafter -> site
xtabs(~ Lot + Wafer, Oxide)



## Create a variable that identifies unique Wafers
Oxide$LotWafer <- with(Oxide, interaction(Lot,Wafer))



## For continuous response 'Thickness',
## fit nested model E[y_{ijk}|lot_i, wafer_j] = a + b_i + g_{ij}
## for Lot i = 1:8 and Wafer j = 1:3 and Site k = 1:3
## where b_i    ~ N(0, sigma_1)
##       g_{ij} ~ N(0, sigma_2)
## and b_i is independent of g_{ij}
## The following four models are identical:

## lme4
lmer(Thickness ~ (1 | Lot/Wafer), data=Oxide)
lmer(Thickness ~ (1 | Lot) + (1 | LotWafer), data=Oxide)
## Note: the equivalence of the above formulations makes
## clear that the intercept indexed by Wafer within Lot
## has the same variance across Lots.

## nlme
lme(Thickness ~ 1, random= ~1 | Lot/Wafer, data=Oxide)
lme(Thickness ~ 1, random=list(~1|Lot, ~1|Wafer), data=Oxide)
## Note: the second formulation illustrates that lme assumes
## nesting in the order that grouping factors are listed. I
## think that this was a poor implementation decision, and
## that the latter should indicate non-nested grouping.




## Fit non-nested model E[y_{ijk}] = a + b_i + g_j
## for Lot i = 1:8 and Wafer j = 1:3 and Site k = 1:3
## where b_i ~ N(0, sigma_1)
##       g_j ~ N(0, sigma_2)
## and b_i is independent of g_j

## lme4
lmer(Thickness ~ (1 | Lot) + (1 | Wafer), data=Oxide)
lmer(Thickness ~ (1 | Wafer) + (1 | Lot), data=Oxide)


## nlme: There is no 'easy' way to do this with nlme,
## and I couldn't get this to work with nlme. This is a
## trick that gives a random slope for each level of the
## grouping variables, which are indexed by the levels of
## a dummy grouping variable with only one group. We also
## specify, for each grouping factor, that covariance
## matrix is proportional to the identity matrix.
Oxide$Dummy <- factor(1)
Oxide <- groupedData(Thickness ~ 1 | Dummy, Oxide)
lme(Thickness ~ 1, data=Oxide,
    random=pdBlocked(list(pdIdent(~ 0 + Lot),
                          pdIdent(~ 0 + Wafer))))



####there are some useful material at
####http://www.geos.ed.ac.uk/homes/lsmallma/r_course.html