addvar_ function(model, addv1 = rep(0,model$rows[1]), 
                    addv2 = rep(0,sum(model$ptab.len) - model$rows[1]), 
                    plotit = F, 
                    labs = c(rep("d",model$rows[1]),
                     rep("c",sum(model$ptab.len)-model$rows[1]),
                     rep("p",dim(model$X)[1]-sum(model$ptab.len))) )

{
    rX1_ model$rows[1]; X _ model$X; rX_ dim(X)[1]
    if (length(addv1) > rX1 || length(c(addv1,addv2)) > rX)
        stop("\nadded variable doesn't have the right length.")
    Y _ model$Y
    G _ model$Ghe * model$g.diag
    E _ Y - X %*% as.matrix(model$mean.theta)	
    sE_ as.vector(G * E)
    temp_ svd(X * G)
    V _ temp$u %*% t(temp$u)			
    addv2_ c(rep(0,rX1),addv2,rep(0,rX-length(addv2)-rX1))
    B_ as.matrix(c(addv1,rep(0,rX-rX1)))+ as.matrix(addv2)
    sB _  B * G					
    avB _ as.vector((diag(rX) - V) %*% sB)
    av.ls_ lsfit(avB,sE)
    if (plotit) {
        plot(avB, sE, xlab="Candidate added variable", 
                ylab="Scaled residuals", type="n")
        if (all(addv1==0) && all(addv2==0))
            warning("added variable has all zero values")
        text(avB, sE, labs, cex = 1)
        abline(av.ls)
    }
    av.lsp_ ls.print(av.ls,print.it=F)
    z_ list(avB = as.vector(avB), sE = as.vector(sE),
            slope = av.lsp$coef.table[2,1], 
            std.err = av.lsp$coef.table[2,2], 
            p.value = av.lsp$coef.table[2,4],)
    z
}