# R/est.R In sasLM: 'SAS' Linear Model

#### Documented in est

```est = function(L, X, rx, conf.level=0.95, adj="lsd", paired=FALSE)
{
nL = nrow(L)
PE = L %*% rx\$coefficients

if (rx\$DFr > 0) {
Var = L %*% rx\$g2 %*% t(L) * rx\$SSE/rx\$DFr
SE = sqrt(diag(Var))
Tval = PE/SE
Pval = 2*(1 - pt(abs(Tval), rx\$DFr))
DL = qt(0.5 + conf.level/2, rx\$DFr)*SE
} else if (tolower(adj) == "tukey") {
if (paired) {
nLevel = round(0.5 + sqrt(1 + 8*nL)/2)
if (abs(nLevel*(nLevel - 1)/2 - nL) > 1e-8) stop("Row count of L does not match for Tukey adjustment!")
} else {
nLevel = nL
}
SE2 = SE/sqrt(2)
Pval = 1 - ptukey(abs(PE/SE2), nLevel, rx\$DFr)
DL = qtukey(conf.level, nLevel, rx\$DFr)*SE2
#    } else if (tolower(adj) == "duncan") {
#      if (paired) {
#        nLevel = round(0.5 + sqrt(1 + 8*nL)/2)
#        if (abs(nLevel*(nLevel - 1)/2 - nL) > 1e-8) stop("Row count of L does not match for Duncan adjustment!")
#      } else {
#        nLevel = nL
#      }
#      SE2 = SE/sqrt(2)
#      Pval = 1 - ptukey(abs(PE/SE2), nLevel, rx\$DFr)^(1/(nLevel - 1))
#      DL = qtukey(conf.level^(nLevel - 1), nLevel, rx\$DFr)*SE2
} else if (substr(tolower(adj), 1, 3) == "bon") {
Pval = nL*2*(1 - pt(abs(Tval), rx\$DFr))
Pval[Pval > 1] = 1
DL = qt(1 - (1 - conf.level)/2/nL, rx\$DFr)*SE
} else if (tolower(adj) == "scheffe") {
DFg = round(0.5 + sqrt(1 + 8*nL)/2) - 1 # Df group
Fs = rep(0, nL)
for (j in 1:nL) Fs[j] = cSS(L[j, , drop=F], rx)[["F value"]]
Fval = Fs/DFg
Pval = 1 - pf(Fval, DFg, rx\$DFr)
DL = sqrt(DFg*qf(conf.level, DFg, rx\$DFr))*SE
} else if (tolower(adj) == "dunnett") {
mC = matrix(0.5, nrow=nL, ncol=nL) + diag(0.5, nrow=nL)
Pval = rep(NA, nL)
DL = rep(NA, nL)
for (k in 1:nL) {
Pval[k] = 1 - pmvt(lower=rep(-abs(Tval[k]), nL), upper=rep(abs(Tval[k]), nL), df=rx\$DFr, corr=mC)
DL[k] = qmvt(0.5 + conf.level/2, df=rx\$DFr, corr=mC)\$quantile*SE[k]
}
} else {
}
LL =  PE - DL
UL =  PE + DL
} else {
SE = NA
Tval = NA
Pval = NA
LL = NA
UL = NA
}

if (tolower(adj) %in% c("lsd", "tukey", "bon", "dunnett")) {
Res = cbind(PE, LL, UL, SE, Tval, rx\$DFr, Pval)
colnames(Res) = c("Estimate", "Lower CL", "Upper CL", "Std. Error", "t value", "Df", "Pr(>|t|)")
attr(Res, "Estimability") = estmb(L, X, rx\$g2)
} else {
Res = cbind(PE, LL, UL, SE, Fval, rx\$DFr, Pval)
colnames(Res) = c("Estimate", "Lower CL", "Upper CL", "Std. Error", "F value", "Df2", "Pr(>F)")
attr(Res, "Estimability") = estmb(L, X, rx\$g2)
}

return(Res)
}
```

## Try the sasLM package in your browser

Any scripts or data that you put into this service are public.

sasLM documentation built on Nov. 10, 2022, 5:07 p.m.