library(nlme)      # {nlme} for lme() & invervals() 
library(ape)       # {ape} for varcomp()
#library(help=ape)  # prototype for finding package index

#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#CALCULATING VARIANCE COMPONENTS

#READ STRUCTURED DATA TABLE WITH NUMERIC CODED NESTED FACTORS
K=read.table("c:/2008LinearModelsData/SRBox10.6NEWR.txt")

K
attach(K)

Y=PH           #response variable
A=factor(Dam)  #OUTER FACTOR A
B=factor(Sire) #INNER FACTOR B inside A

r=dim(K)[[1]]  #TOTAL NUMBER OF INDIVIDUALS
r              #quantity 1 in SR

BB=summary(B)
BBsq=BB^2
q2=sum(BBsq)   
q2             #quantity 2 in SR

a=length(levels(A))
AA=summary(A)
AAsq=AA^2
q3=sum(AAsq)   
q3             #quantity 3 in SR

b=length(levels(B))
L=table(A,B)
BBsum=rowSums(L^2)
M=BBsum/AA
q4=sum(M) 
q4             #quantity 4 in SR

k1= (1/(a*((b/a)-1)))*(r-q4)
k1             #quantity k1 in LW

k2= (1/(a-1))*(q4-(q2/r))
k2             #quantity k2 in LW

k3= (1/(a-1))*(r-(q3/r))
k3             #quantity k3 in LW

FM=lm(Y~A+B%in%A)
anova(FM)  #standard Least Squares NESTED ANOVA
MSA = anova(FM)[1,3]
MSB = anova(FM)[2,3]
MSE = anova(FM)[3,3]
MS=c(MSA,MSB,MSE)

#LEAST SQUARES VARIANCE COMPONENTS:
sigA = (MSA-MSE-(k2/k1)*(MSB-MSE))/k3
sigA  #variance component for outside factor A

sigB = (MSB-MSE)/k1
sigB  #variance component for nested factor B inside A

sigE = MSE
sigE  #variance component for error
LEASTSQvar=c(sigA,sigB,sigE)

#MAXIMUM LIKLIHOOD VARIANCE COMPONENTS using lme():
FMe=lme(Y~1,random=~1|A/B,data=K)
FMe
summary(FMe)
varcomp(FMe)    # {ape} Variance Components calculated
MAXLIKvar=c(varcomp(FMe)[1],varcomp(FMe)[2],varcomp(FMe)[3])

intervals(FMe)  # {nlme} PB Confidence Intervals
PBINTlower=c(intervals(FMe)$reStruct$A[1,1]^2,
              intervals(FMe)$reStruct$B[1,1]^2,
              intervals(FMe)$sigma[1]^2)
PBINTest  =c(intervals(FMe)$reStruct$A[1,2]^2,
              intervals(FMe)$reStruct$B[1,2]^2,
              intervals(FMe)$sigma[2]^2)
PBINTupper = c(intervals(FMe)$reStruct$A[1,3]^2,
              intervals(FMe)$reStruct$B[1,3]^2,
              intervals(FMe)$sigma[3]^2)

results=cbind(MS,LEASTSQvar,MAXLIKvar,PBINTlower,PBINTest,PBINTupper)


#LEAST SQUARES AND MAXIMUM LIKLIHOOD ESTIMATES OF VARIANCE COMPONENTS:
results
detach(K)