R/fgls.iter.siga.univ.R

In dmm: Dyadic Mixed Model for Pedigree Data

fgls.iter.siga.univ <-
function(am, start.siga,mmat,dyad.explist,evec,dmefglsopt,dmeopt,ctable){
# fgls.iter.siga() - iterate and compute new siga from DME by fgls
# Univariate version
    stopcrit <- 10
    stoptol <- dmefglsopt$stoptol *0.1
    maxiter <- dmefglsopt$maxiter
    bdamp <- dmefglsopt$bdamp
    newsiga <- start.siga
    count <- 0
    degf <- am$n - am$k
    bias <- am$n/degf
    ipkminv <- ginvipk(am$n,am$n)
    while (stopcrit > stoptol && count < maxiter) {
#       cat("Iteration round: ", count, "\n")
        oldsiga <- newsiga
#    compute V matrix
        v <- expect.v(am,oldsiga,dyad.explist) # v is lxl blocks each nxn
#       cat("V matrix:\n")
#       print(v)
        vinv <- ginv(v)
#    solve for siga
        newsiga.list <- fgls.update.univ(am, v, mmat, dyad.explist, evec, ipkminv)
        newsiga <- newsiga.list$siga
#    damp the siga update
        newsiga <- oldsiga + (newsiga - oldsiga) * bdamp
#       cat("Updated siga:\n")
#       print(newsiga)

#  check updated siga posdef
        newsiga <- siga.posdef(newsiga, am, ctable)
#       cat("Updated siga made positive definite:\n")
#       print(newsiga)

#    look at stopcrit
        sumdev <- 0
        for (ll in 1:(am$l*am$l)) {
            for (i in 1:am$v) {
                sumdev <- sumdev + abs(newsiga[i, ll] - oldsiga[i, ll])
                # (newsiga - oldsiga) is always (new-old)*bdamp - ie part of the difference
            }
        }
        stopcrit <- sumdev/(am$l * am$l * am$v)
#       cat("stopcrit = ", stopcrit, "\n")
        count <- count + 1
        cat("Round = ",count," Stopcrit = ",stopcrit,"\n")
    }
#     end of iteration
      cat("Iteration completed - count = ",count,"\n")
#   convergence check
      if(count == maxiter){
        cat("Failed to converge\n")
        return()
      }
      cat("Convergence achieved\n")

#   SE of siga 
    degfd <- am$n * am$n - am$v
#     vsiga <- kronecker(vard, solve(crossprod(qr.R(newemat.qr))), make.dimnames=T)
      vsiga <- newsiga.list$vsiga  # may be univariate ?
      sesiga <- matrix(sqrt(diag(vsiga)), am$v, am$l * am$l, dimnames=dimnames(start.siga))

    parlist <- list(  siga=newsiga,sesiga=sesiga, vsiga=vsiga) #  vard.ols=vard, msr.ols=msr, msrdf=degf, vara.ols=msa)
    return(parlist)
}