### 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 F_{2} 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)

}