ds=read.table("otitis-data-R.txt",head=T)
ds

eggs=1

d=ds[ds$z==eggs,2:3] ; d # 1st 2 monozygotic twins

k=length(d[,1]); n=2; k ; n.. = k * n

y=c(d[,1],d[,2])

ybar.. = mean(y) ; ybar..

ss.overall = sum((y - ybar..)^2) ; ss.overall 

   # check: (n..-1)*var(y)

ybar.i=(d[,1]+d[,2])/2 ; hist(ybar.i)

ss.B = n*sum((ybar.i - ybar..)^2) ; ss.B

   # check: 2*(k-1)*var(ybar.i)

ss.W = sum((d[,1]-ybar.i)^2+(d[,2]-ybar.i)^2); ss.W

   # check: (n-1)*sum(apply(d,1,var))
   
   # v.i=((d[,1]-ybar.i)^2+(d[,2]-ybar.i)^2)/(2-1) ;
   # hist(v.i)

ms.B=ss.B/(k-1) ; ms.W=ss.W/(n..-k) ; c(ms.B,ms.W)

sigma.sq.W.hat = ms.W ; sigma.sq.W.hat
sigma.sq.B.hat = (ms.B - ms.W)/2 ; sigma.sq.B.hat