REML<-function(y,X,n_comp,conv_crit,n_maxiter,lambda_start,delta, Z1, Z2=0, print_results=FALSE, step=1) { ###################################################################################### ## Responsible programmer: Lars Rönnegård (lars.ronnegard@lcb.uu.se) ## Date 20060925 ## Code published in appendix of Rönnegård & Carlborg 2007, BMC Genetics 8:1 #This is the Fisher Scoring algorithm for REML #See e.g. Johnson & Thompson 1995 J. of Dairy Science 78:449-456. #The function returns the estimated variance components and fixed effects, together with the log-likelihood #and convergence information. #INPUT PARAMTERS #y Response vector #X Design matrix for fixed effects #n_comp Number of different random effects in the model (max.=2 in this version) #conv_crit Value that the change in loglikelihood should be less than #n_maxiter Maximum number of iterations #lambda_start Initial ratio of variance components #delta Minimum possible value of the variance component #WORKING VARIABLES AND OUTPUT PARAMETERS #A A matrix which stores ZZ’ for all variance components #phi_start Staring values for the variance components #M_phi Matrix with VC estimates at each iteration #phi VC estimate from the latest iteration #DL Gradient of the restricted log-likelihood #FS Fisher’s Information matrix #V Variance matrix of y #P The projection matrix #llh Restricted log-likelihood #beta_hat Estimates of fixed effects #conv_test Binary variable equal to 1 if the algorithm converges within n_maxiter iterations ###################################################################################### print("REML iteration number") min.error<-10^-8 n_comp1<-n_comp+1 n_rows<-length(y) A<-matrix(0,(n_comp1*n_rows),n_rows) Aj<-matrix(0,n_rows,n_rows) for (i_comp in 1:n_comp) { if (i_comp==1) Aj<-Z1%*%t(Z1) if (i_comp==2) Aj<-Z2%*%t(Z2) A[((i_comp-1)*n_rows+1):(n_rows*i_comp),1:n_rows]<-Aj } A[(n_comp*n_rows+1):(n_comp1*n_rows),1:n_rows]<-diag(rep(1,n_rows)) res_var<-var(y-X%*%solve(t(X)%*%X)%*%t(X)%*%y) phi_start<-numeric(n_comp1) for (i in 1:(n_comp1-1)) { phi_start[i]<-lambda_start*res_var } phi_start[n_comp1]<-res_var dimIBD<-min(dim(A)) M_phi<-matrix(0,(n_maxiter+1),n_comp1) M_phi[1,1:n_comp1]<-phi_start[1:n_comp1] phi<-numeric(n_comp1) DL<-numeric(n_comp1) DL[1:n_comp1]<-conv_crit+1 FS<-matrix(0,n_comp1,n_comp1) Aj<-matrix(0,dimIBD,dimIBD) Ak<-matrix(0,dimIBD,dimIBD) llh.prev<-1+conv_crit llh<-0 for (i in 1:n_maxiter){ V<-matrix(0,dimIBD,dimIBD) if (abs(llh-llh.prev)>conv_crit){ phi<-M_phi[i,] for (j in 1:n_comp1){ Aj<-A[((j-1)*dimIBD+1):(j*dimIBD),1:dimIBD] V<-V+phi[j]*Aj } invV<-solve(V) temp<-solve(t(X)%*%invV%*%X) P<-invV-invV%*%X%*%(temp)%*%t(X)%*%invV for (j in 1:n_comp1){ Aj<-A[((j-1)*dimIBD+1):(j*dimIBD),1:dimIBD] DL[j]<-sum(diag(Aj%*%P))-t(y)%*%P%*%Aj%*%P%*%y for (k in j:n_comp1) { Ak<-A[((k-1)*dimIBD+1):(k*dimIBD),1:dimIBD] FS[j,k]<-sum(diag(Aj%*%P%*%Ak%*%P)) FS[k,j]<-FS[j,k] } } FS.egen<-eigen(FS, only.values=TRUE) FS.min<-min(FS.egen$values) if (FS.min>min.error) M_phi[(i+1),]<-phi-step*solve(FS)%*%DL if (FS.min<=min.error) { print("Negative Hessian") identitet<-diag(rep((0.5+abs(FS.min)),min(dim(FS)))) M_phi[(i+1),]<-phi-step*solve(FS+identitet)%*%DL } egen<-eigen(V, only.values=TRUE) if (min(egen$values>0))ldV<-sum(log(egen$values)) egen2<-eigen(t(X)%*%invV%*%X, only.values=TRUE) ldXVX<-sum(log(egen2$values)) llh.prev<-llh if (min(egen$values)>0) llh<-(ldV+ldXVX+t(y)%*%P%*%y)*(-0.5) if (min(egen$values)<=0) llh<-llh.prev-1 #Truncation at zero M_phi[(i+1),]<-0.5*(M_phi[(i+1),]+delta+abs(M_phi[(i+1),]-delta)) conv_val<-abs(llh-llh.prev) if (print_results==TRUE) { print("--------------------------------") print("Iteration:") print(i) print("Convergence criteria: Change in log-likelihood") print(conv_val) print("log-likelihood") print(llh) print("REML estimates of variance components") print("Genotype variance [1] and residual variance [2]") print(M_phi[(i+1),]) } if (print_results==FALSE) print(paste(" ",i)) } } conv_test<-1 if (abs(llh-llh.prev)>conv_crit) conv_test<-0 beta_hat<-numeric(min(dim(X))) beta_hat<-solve(t(X)%*%invV%*%X)%*%t(X)%*%invV%*%y if (print_results==TRUE) { print("Estimates of fixed effects") print(beta_hat) } print("log-likelihood") print(llh) print("Estimates of fixed effects") print(beta_hat) print("REML estimates of variance components") print("Estimate of variances for the random effects") print(phi[(1:n_comp)]) print("Estimate of residual variance") print(phi[n_comp1]) list(beta_hat=beta_hat,conv_test=conv_test,conv_val=conv_val,phi=phi,phi_iteration=M_phi,llh=llh) }