Fisher scoring algorithm for REML estimates in R code
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
#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 (Table 5)
#WORKING VARIABLES AND OUTPUT PARAMETERS (Table 6)
######################################################################################
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)
}
list(beta_hat=beta_hat,conv_test=conv_test,conv_val=conv_val,phi=phi,phi_iteration=M_phi,llh=llh)
}
Algorithm for construction of an F2 pedigree with a fully informative marker in R code
simple_F2<-function(n0.males,n0.females,n1.males,n1.females,n2.males,n2.females,QTLvar,RESvar) {
######################################################################################
## Responsible programmer: Lars Rönnegård (lars.ronnegard@lcb.uu.se)
## Date 20060925
## This function constructs an F2 pedigree for a fully informative marker.
## The function returns the simulated phenotypes (y) and the incidence matrix (Z)
## and also the true simulated QTL alleles of the base individuals.
## The phenotypes are simulated assuming a biallelic QTL at the marker position
#INPUT PARAMTERS (Table 9)
#WORKING VARIABLES AND OUTPUT PARAMETERS (Table 10)
######################################################################################
n0=n0.males+n0.females
n1=n1.males+n1.females
n2=n2.males+n2.females
n=n0+n1+n2
Z<-matrix(0,n2,(2*n0))
index.mat<-matrix(0,n1,2)
for (i in 1:n1) {
sire<-i-n0.males*floor(i/n0.males)+1
sire.allele<-sample(1:2,1)
dam<-i-n0.females*floor(i/n0.females)+1
dam=dam+n0.males
dam.allele<-sample(1:2,1)
index.mat[i,]<-c(((sire-1)*2+sire.allele),((dam-1)*2+dam.allele))
}
for (i in 1:n2) {
sire<-sample(1:n1.males,1)
sire.allele<-sample(1:2,1)
dam<-sample(1:n1.females,1)
dam=dam+n1.males
dam.allele<-sample(1:2,1)
Z[i,index.mat[sire,sire.allele]]=Z[i,index.mat[sire,sire.allele]]+1
Z[i,index.mat[dam,dam.allele]]=Z[i,index.mat[dam,dam.allele]]+1
}
base.alleles <-sample(0:1,(2*n0),replace=TRUE)
v<-base.alleles *sqrt(2*QTLvar)
e<-rnorm(n2,0,sqrt(RESvar))
y<-Z%*%v+e
list(y=y,Z=Z,base.alleles= base.alleles)
}