# R/racov.r In npsm: Package for Nonparametric Statistical Methods using R

#### Documented in getxactgetxact2kancovakancovarownonecovaonecovaheteronecovahomog

```onecovahomog = function(levs,data,xcov,print.table=TRUE) {
#
#   Input:
#     levs = vector of levels corresponding to the factors A, B, C, etc.
#     data = matrix whose first column is the response variable then the succeeding columns
#            are the level of the first factor, second factor, etc.
#            For example, the Data
#              for a three way is of the form:
#              Y_ijkl i j k
#            NOTE: The data must be presorted so that k runs the fastest, than j, than i
#     CONTRASTS = T if confidence intervals are to be obtained.
#     CONTRASTSCONF is the confidence coefficient for the confidence intervals for
#                   the contrasts.
#     CONTRASTSMAT is the matrix of contrasts, one row for each contrast.
#

xcov = as.matrix(xcov)
nf = length(levs)
y = data[,1]
levsind = data[,2:(nf+1)]
n = length(y)
xfull = cellx(data[,2:(nf+1)])
#    fitcell = rfit(y~xfull,intercept=F)
fitcell = rfit(y~xfull-1)
dercell = disp(fitcell\$betahat,xfull,fitcell\$y,fitcell\$scores)
pcell = length(xfull[1,])
tab2 = matrix(rep(0,10),ncol=5)

xfull2 = cbind(xfull,xcov)
p = length(xfull2[1,])
dffull = n - p
#    fitF = rfit(y~xfull2,intercept=FALSE)
fitF = rfit(y~xfull2-1)
tau = fitF\$tauhat
derfull = disp(fitF\$betahat,xfull2,fitF\$y,fitF\$scores)

pcov = length(xcov[1,])
fitcov = rfit(y~xcov)
betacov = fitcov\$betahat[2:(pcov+1)]
dercov = disp(betacov,xcov,fitcov\$y,fitcov\$scores)
dropdf = p - (pcov+1)
rdcell = dercov - derfull
rdcellq = rdcell/dropdf
testcell = rdcellq/(tau/2)
pvcell = 1 - pf(testcell,dropdf,dffull)
tab2[1,] = c(dropdf,rdcell,rdcellq,testcell,pvcell)

dropdf = p - pcell
rdcov = dercell - derfull
rdcovq = rdcov/dropdf
testcov = rdcovq/(tau/2)
pvcov = 1 - pf(testcov,dropdf,dffull)
tab2[2,] = c(dropdf,rdcov,rdcovq,testcov,pvcov)
rownames(tab2) = c("Groups","Covariates")
colnames(tab2) = c("df","RD","MRD","F","p-value")
if(print.table){
cat("\n")
cat("Robust ANCOVA (Assuming Same Slopes) Table","\n")
print(tab2)
cat("\n")
}

list(tab=tab2,fit=fitF)
#         list(tab2=tab2)
}

onecovaheter = function(levs,data,xcov,print.table=TRUE) {
#
#   Input:
#     levs = vector of levels corresponding to the factors A, B, C, etc.
#     data = matrix whose first column is the response variable then the succeeding columns
#            are the level of the first factor, second factor, etc.
#            For example, the Data
#              for a three way is of the form:
#              Y_ijkl i j k
#            NOTE: The data must be presorted so that k runs the fastest, than j, than i
#     CONTRASTS = T if confidence intervals are to be obtained.
#     CONTRASTSCONF is the confidence coefficient for the confidence intervals for
#                   the contrasts.
#     CONTRASTSMAT is the matrix of contrasts, one row for each contrast.
#

xcov = as.matrix(xcov)
nf = length(levs)
y = data[,1]
levsind = data[,2:(nf+1)]
n = length(y)
xfull = cellx(data[,2:(nf+1)])
#    fitcell = rfit(y~xfull,intercept=F)
fitcell = rfit(y~xfull-1)
dercell = disp(fitcell\$betahat,xfull,fitcell\$y,fitcell\$scores)
pcell = length(xfull[1,])
tab2 = matrix(rep(0,10),ncol=5)

#   reduced [cell xcov]
xfull2 = cbind(xfull,xcov)
p = length(xfull2[1,])
#    fitF = rfit(y~xfull2,intercept=FALSE)
fitF = rfit(y~xfull2-1)
tau = fitF\$tauhat
derfull = disp(fitF\$betahat,xfull2,fitF\$y,fitF\$scores)

#   full model   [cell xcov Interact]
xint = getxact(xfull,xcov)
pint = length(xint[1,])
#    fint = rfit(y~xint,intercept=F)
fint = rfit(y~xint-1)
dint = disp(fint\$betahat,xint,fint\$y,fint\$scores)
tau = fint\$tauhat
fulldf = n - pint

#   Test homog slopes
dropdf = pint - p
rdhomg = derfull - dint
rdhomgq = rdhomg/dropdf
fhomg = rdhomgq/(tau/2)
pval = 1 -pf(fhomg,dropdf,fulldf)
tab2[2,] = c(dropdf,rdhomg,rdhomgq,fhomg,pval)

#  Test for groups  reduced [xco
xuse = xint[,(pcell+1):pint]
puse = length(xuse[1,])
fuse = rfit(y~xuse)
beta = fuse\$betahat[2:(puse+1)]
duse = disp(beta,xuse,fuse\$y,fuse\$scores)
rdcell = duse - dint
numdf = pcell - 1
rdcellq = rdcell/numdf
fcell = rdcellq/(tau/2)
pval = 1 - pf(fcell,numdf,fulldf)
tab2[1,] = c(numdf,rdcell,rdcellq,fcell,pval)

rownames(tab2) = c("Groups","Homog Slopes")
colnames(tab2) = c("df","RD","MRD","F","p-value")
if(print.table){
cat("\n")
cat("Robust ANCOVA (Assuming Heterogeneous Slopes) Table","\n")
print(tab2)
cat("\n")
}

list(tab=tab2,fit=fint)
#         list(tab2=tab2)
}
getxact = function(amat,bmat){
#
#   Note that the first group is referenced.
#
ca = length(amat[1,])
cb = length(bmat[1,])

cmat = cbind(amat,bmat)
for(i in 2:ca){
for(j in 1:cb){
cmat = cbind(cmat,amat[,i]*bmat[,j])
}
}
cmat
}
getxact2 = function(amat,bmat){
#
#
ca = length(amat[1,])
cb = length(bmat[1,])

cmat = amat
for(i in 1:ca){
for(j in 1:cb){
cmat = cbind(cmat,amat[,i]*bmat[,j])
}
}
nm = rep(0,ca)
for(i in 1:ca){nm[i] = paste("Int ",i)}
for(i in 1:cb){
for(j in 1:ca){
nm =c(nm,paste("rg ",j,i))
}
}
colnames(cmat) = nm
cmat
}
kancova = function(levs,data,xcov,print.table=TRUE) {
#
#   Input:
#     levs = vector of levels corresponding to the factors A, B, C, etc.
#     data = matrix whose first column is the response variable then the succeeding columns
#            are the level of the first factor, second factor, etc.
#            For example, the Data
#              for a three way is of the form:
#              Y_ijkl i j k
#            NOTE: The data must be presorted so that k runs the fastest, than j, than i
#     CONTRASTS = T if confidence intervals are to be obtained.
#     CONTRASTSCONF is the confidence coefficient for the confidence intervals for
#                   the contrasts.
#     CONTRASTSMAT is the matrix of contrasts, one row for each contrast.
#

xcov = as.matrix(xcov)
nf = length(levs)
y = data[,1]
levsind = data[,2:(nf+1)]
n = length(y)
xfull = cellx(data[,2:(nf+1)])
#    fitcell = rfit(y~xfull,intercept=F)
fitcell = rfit(y~xfull-1)
dercell = disp(fitcell\$betahat,xfull,fitcell\$y,fitcell\$scores)
pcell = length(xfull[1,])
xfull2 = cbind(xfull,xcov)
listtests = subsets(nf)
p = length(xfull2[1,])
pcov = length(xcov[1,])
#     fitF = rfit(y~xfull2,intercept=FALSE)
fitF = rfit(y~xfull2-1)
iflagq = 0

drf = disp(fitF\$betahat, xfull2, fitF\$y, fitF\$scores)
ehatr = fitF\$resid
dfull = n - p
nt = length(listtests[,1])
tab2 = matrix(rep(0,(nt+2)*5),ncol=5)
tab3 = matrix(rep(0,4),ncol=2)
rnm = c("0")

for(i in 1:nt){
permh = listtests[i,]
vecl = listtests[i,]
nm = kancovarown(vecl)
rnm = c(rnm,nm)
hmat = khmat(levsind,permh)
q = length(hmat[,1])
mat0 = matrix(rep(0,q*pcov),ncol=pcov)
hmat = cbind(hmat,mat0)
#      see = cbind(rep(i,q),hmat)
#      write(t(see),ncol=13,append=T,file="allh.dat")
xred = redmod(xfull2,hmat)
#      fitr = rfit(y~xred,intercept=FALSE)
fitr = rfit(y~xred-1)
drr = disp(fitr\$betahat, xred, fitr\$y, fitr\$scores)
rd = drr - drf
ft = (rd/q)/(fitF\$tauhat/2)
pv = 1 - pf(ft,q,dfull)
tab2[i,] = c(q,rd,(rd/q),ft,pv)
#      tab3[i,] = c(dfull,0,(fitF\$tauhat/2),drr,drf)
}
tab3[1,] = c(dfull,(fitF\$tauhat/2))
rdcell = dercell - drf
fcell = (rdcell/pcov)/(fitF\$tauhat/2)
pcell = 1 - pf(fcell,pcov,dfull)
tab2[(nt+1),] = c(pcov,rdcell,(rdcell/pcov),fcell,pcell)
rnm = c(rnm,"Covariate")
xint = getxact(xfull,xcov)
pint = length(xint[1,])
#      fint = rfit(y~xint,intercept=F)
fint = rfit(y~xint-1)
dint = disp(fint\$betahat,xint,fint\$y,fint\$scores)
rdint = drf - dint
ardint = rdint/(pint-p)
finttest = ardint/(fint\$tauhat/2)
pvint = 1 - pf(finttest,(pint-p),(n-pint))
tab2[(nt+2),] = c((pint-p),rdint,ardint,finttest,pvint)
rnm = c(rnm,"Hetrog regr")
rnml = length(rnm)
rnm = rnm[2:rnml]
rownames(tab2) = rnm
colnames(tab2) = c("df","RD","MRD","F","p-value")
if(print.table){
cat("\n")
cat("Robust ANCOVA Table","\n")
cat("All tests except last row is with homogeneous slopes","\n")
cat("as the full model.   For the last row the full model is","\n")
cat("with heterscedastic slopes.","\n")
print(tab2)
cat("\n")
}
tab3[2,] = c((n-pint),(fint\$tauhat/2))
#  xint contains the augmented matrix of cell mean model and interaction columns
#  However the first col of int part is full vector(s) of covariates and the
#  rest are increments.

#list(tab=tab2,fit=fitF)
list(listtests=listtests,tab2=tab2,tab3=tab3,xint=xint,fithomog=fitF,fint=fint)

}
kancovarown = function(vec){
num = length(vec)
nm = c(round(vec[1]))
for(j in 2:num){
nm = paste(nm,",",round(vec[j]))
}
nm
}

onecova = function(levs,data,xcov,print.table=TRUE) {
#
#   Input:
#     levs = vector of levels corresponding to the factors A, B, C, etc.
#     data = matrix whose first column is the response variable then the succeeding columns
#            are the level of the first factor, second factor, etc.
#            For example, the Data
#              for a three way is of the form:
#              Y_ijkl i j k
#            NOTE: The data must be presorted so that k runs the fastest, than j, than i
#     CONTRASTS = T if confidence intervals are to be obtained.
#     CONTRASTSCONF is the confidence coefficient for the confidence intervals for
#                   the contrasts.
#     CONTRASTSMAT is the matrix of contrasts, one row for each contrast.
#

xcov = as.matrix(xcov)
nf = length(levs)
y = data[,1]
levsind = data[,2:(nf+1)]
n = length(y)
xfull = cellx(data[,2:(nf+1)])
#    fitcell = rfit(y~xfull,intercept=F)
fitcell = rfit(y~xfull-1)
dercell = disp(fitcell\$betahat,xfull,fitcell\$y,fitcell\$scores)
pcell = length(xfull[1,])
xfull2 = cbind(xfull,xcov)
p = length(xfull2[1,])
pcov = length(xcov[1,])
#     fitF = rfit(y~xfull2,intercept=FALSE)
fitF = rfit(y~xfull2-1)
iflagq = 0

drf = disp(fitF\$betahat, xfull2, fitF\$y, fitF\$scores)
ehatr = fitF\$resid
dfull = n - p
nt = 1
tab2 = matrix(rep(0,(nt+2)*5),ncol=5)
tab3 = matrix(rep(0,4),ncol=2)

hmat = cbind(rep(-1,(pcell-1)),diag(rep(1,(pcell-1))))
q = length(hmat[,1])
mat0 = matrix(rep(0,q*pcov),ncol=pcov)
hmat = cbind(hmat,mat0)
#      see = cbind(rep(i,q),hmat)
#      write(t(see),ncol=13,append=T,file="allh.dat")
xred = redmod(xfull2,hmat)
#      fitr = rfit(y~xred,intercept=FALSE)
fitr = rfit(y~xred-1)
drr = disp(fitr\$betahat, xred, fitr\$y, fitr\$scores)
rd = drr - drf
ft = (rd/q)/(fitF\$tauhat/2)
pv = 1 - pf(ft,q,dfull)
tab2[1,] = c(q,rd,(rd/q),ft,pv)
tab3[1,] = c(dfull,(fitF\$tauhat/2))

rdcell = dercell - drf
fcell = (rdcell/pcov)/(fitF\$tauhat/2)
pcell = 1 - pf(fcell,pcov,dfull)
tab2[(nt+1),] = c(pcov,rdcell,(rdcell/pcov),fcell,pcell)
xint = getxact(xfull,xcov)
pint = length(xint[1,])
#      fint = rfit(y~xint,intercept=F)
fint = rfit(y~xint-1)
dint = disp(fint\$betahat,xint,fint\$y,fint\$scores)
rdint = drf - dint
ardint = rdint/(pint-p)
finttest = ardint/(fint\$tauhat/2)
pvint = 1 - pf(finttest,(pint-p),(n-pint))
tab2[(nt+2),] = c((pint-p),rdint,ardint,finttest,pvint)
tab3[2,] = c((n-pint),(fint\$tauhat/2))
x4 = getxact2(xfull,xcov)
f4 = rfit(y~x4-1)
rownames(tab2) = c("Groups","Covariates","G:C")
colnames(tab2) = c("df","RD","MRD","F","p-value")
if(print.table){
cat("\n")
cat("Robust ANCOVA Table","\n")
print(tab2)
cat("\n")
}

#list(tab=tab2,fit=fitF)
list(tab2=tab2,tab3=tab3,fitF1=fitF,fitF2=f4,xint=xint)

}
```

## Try the npsm package in your browser

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

npsm documentation built on May 29, 2017, 12:12 p.m.