resd_ function(model,plotit=T,
               labs=rep("d",model$rows[1])   )
{
    X_ model$X; G_ model$Ghe
    rX1_ model$rows[1]
    yhat_ as.vector(X %*% model$mean.theta)
    rX _ dim(X)[1]
    E_ model$Y - yhat
    svX_ svd(X * G)
    V  _ svX$u %*% t(svX$u)
    sE _ G * E
    stud.res_ sE/sqrt(diag(diag(rX)-V))
    stud.res1_ stud.res[1:rX1];  stud.res2_ stud.res[(rX1+1):rX]
    if (plotit) {
        par(mfrow=c(1,2),pty="s")
        qqnorm(stud.res[1:rX1], main="the data cases")
        abline(1,lty=3)  
        qqnorm(stud.res[(rX1+1):sum(model$ptab.len)], 
                            main="the constraint cases")
        abline(1,lty=3)  
        par(mfrow=c(1,1),pty="m",ask=T)
        plot(yhat[1:rX1],stud.res1, xlab="fitted value",
        ylab="scaled data-case residual",type="n")
        text(yhat,stud.res, labs, cex = 1)
        fit.sry_ lsfit(yhat,stud.res)
        abline(fit.sry)
    }
    z_ list(yhat=as.vector(yhat),stud.res=stud.res,
            slope=1-fit.sry$coef[2],
            std.err=ls.diag(fit.sry)$std.err[2])
    z
}