RAZAW<-function(Y,X,Z1,Z2) { ############################################# # This function calculates Henderson 3 estimates and modified Hendrson 3 estimates for a # mixed linear model with two random effects and iid residuals. # # Code written by Lars Ronnegard 2008-07-29 # # Input: # Y = response vector # X = design matrix for fixed effects # Z1 = incidence matrix for first random effect # Z2 = incidence matrix for second random effect ############################################# library(MASS) n<-length(Y) ############################################# X1<-cbind(X,Z1) X2<-cbind(X1,Z2) PX<-X%*%ginv(t(X)%*%X)%*%t(X) PX1<-X1%*%ginv(t(X1)%*%X1)%*%t(X1) PX2<-X2%*%ginv(t(X2)%*%X2)%*%t(X2) I<-diag(1,length(Y)) A<-PX1-PX B<-PX2-PX1 C<-I-PX2 V1<-Z1%*%t(Z1) V2<-Z2%*%t(Z2) a<-sum(diag((PX1-PX)%*%V1)) b<-sum(diag((PX2-PX1)%*%V2)) c<-sum(diag(I-PX2)) d<-sum(diag((PX1-PX)%*%V2)) e<-sum(diag(PX2-PX1)) #Added f to equations f<-sum(diag(PX2-PX)) #Used lower case k since it is a scalar k<-d*e-f*b ############ c1<-1/(2/(a^2)*sum(diag(A%*%V1%*%A%*%V1))+1) e1<-1/(2/(b^2)*sum(diag(B%*%V2%*%B%*%V2))+1) e2<-(d/b*e1*sum(diag(B))-sum(diag(A)))/((k/b)*(2/c+1)) # Modified.Henderson3<-numeric(3) Original.Henderson3<-numeric(3) Modified.Henderson3[1]<-c1/a*(t(Y)%*%A%*%Y-d/b*e1*(t(Y)%*%B%*%Y)+k/(b*c)*e2*(t(Y)%*%C%*%Y)) Original.Henderson3[1]<-1/a*(t(Y)%*%A%*%Y-d/b*(t(Y)%*%B%*%Y)+k/(b*c)*(t(Y)%*%C%*%Y)) Modified.Henderson3[2]<-c1*e1/b*(t(Y)%*%B%*%Y)-c1*e2/(b*c)*(t(Y)%*%C%*%Y) Original.Henderson3[2]<-1/b*(t(Y)%*%B%*%Y)-1/(b*c)*(t(Y)%*%C%*%Y) Modified.Henderson3[3]<-c1*e2/c*(t(Y)%*%C%*%Y) Original.Henderson3[3]<-1/c*(t(Y)%*%C%*%Y) print(" Modified VC estimates are given as $Modified.Henderson3") print("Non-modified VC estimates are given as $Original.Henderson3") list(Original.Henderson3=Original.Henderson3,Modified.Henderson3=Modified.Henderson3,c1=c1,e1=e1,e2=e2,k=k) }