kbcon <-
function(x,con=0,tr=.2,alpha=.05,pr=T){
#
# A heteroscedastic test of d linear contrasts using trimmed means.
# based on the Kaiser-Bowden method
#
# The data are assumed to be stored in $x$ in list mode or in a matrix
# Length(x) (or ncol(x))is assumed to correspond to the total number of groups
# It is assumed all groups are independent.
#
# con is a J by d matrix containing the contrast coefficients.
# If con is not specified, all pairwise comparisons are made.
#
# Missing values are automatically removed.
#
if(is.matrix(x))x<-listm(x)
if(!is.list(x))stop("Data must be stored in a matrix or in list mode.")
con<-as.matrix(con)
J<-length(x)
h<-vector("numeric",J)
w<-vector("numeric",J)
xbar<-vector("numeric",J)
for(j in 1:J){
xx<-!is.na(x[[j]])
val<-x[[j]]
x[[j]]<-val[xx] # Remove missing values
h[j]<-length(x[[j]])-2*floor(tr*length(x[[j]]))
# h is the number of observations in the jth group after trimming.
w[j]<-((length(x[[j]])-1)*winvar(x[[j]],tr))/(h[j]*(h[j]-1))
xbar[j]<-mean(x[[j]],tr)
}
if(sum(con^2)==0){
CC<-(J^2-J)/2
v1=J-1
psihat<-matrix(0,CC,6)
dimnames(psihat)<-list(NULL,c("Group","Group","psihat","ci.lower","ci.upper",
"p.value"))
test<-matrix(NA,CC,6)
dimnames(test)<-list(NULL,c("Group","Group","test","crit","se","df"))
jcom<-0
for (j in 1:J){
for (k in 1:J){
if (j < k){
jcom<-jcom+1
df<-(w[j]+w[k])^2/(w[j]^2/(h[j]-1)+w[k]^2/(h[k]-1))
term=sqrt((J-1)*(1+(J-2)/df))
test[jcom,3]<-((xbar[j]-xbar[k])/(term*sqrt(w[j]+w[k])))^2
sejk<-sqrt(w[j]+w[k])
test[jcom,5]<-sejk
psihat[jcom,1]<-j
psihat[jcom,2]<-k
test[jcom,1]<-j
test[jcom,2]<-k
psihat[jcom,3]<-(xbar[j]-xbar[k])
test[jcom,6]<-df
crit=qf(1-alpha,v1,df)
aterm=sqrt(crit)*term
test[jcom,4]<-qf(1-alpha,v1,df)
psihat[jcom,4]<-(xbar[j]-xbar[k])-aterm*sejk
psihat[jcom,5]<-(xbar[j]-xbar[k])+aterm*sejk
psihat[jcom,6]<-1-pf(test[jcom,3],v1,df)
}}}}
if(sum(con^2)>0){
if(nrow(con)!=length(x)){
stop("The number of groups does not match the number of contrast coefficients.")
}
v1=nrow(con)-1
psihat<-matrix(0,ncol(con),5)
dimnames(psihat)<-list(NULL,c("con.num","psihat","ci.lower","ci.upper",
"p.value"))
test<-matrix(0,ncol(con),5)
dimnames(test)<-list(NULL,c("con.num","test","crit","se","df"))
df<-0
L=nrow(con)
for (d in 1:ncol(con)){
psihat[d,1]<-d
psihat[d,2]<-sum(con[,d]*xbar)
sejk<-sqrt(sum(con[,d]^2*w))
test[d,1]<-d
df<-(sum(con[,d]^2*w))^2/sum(con[,d]^4*w^2/(h-1))
A=(L-1)*(1+(L-2)/df)
test[d,2]<-(sum(con[,d]*xbar)/(sqrt(A)*sejk))^2
crit=qf(1-alpha,v1,df)
test[d,3]<-crit
test[d,4]<-sejk
test[d,5]<-df
psihat[d,3]<-psihat[d,2]-sqrt(crit*A)*sejk
psihat[d,4]<-psihat[d,2]+sqrt(crit*A)*sejk
psihat[d,5]<-1-pf(test[d,2],v1,df)
}}
#
if(pr){
print("Note: confidence intervals are adjusted to control FWE")
print("But p-values are not adjusted to control FWE")
}
list(test=test,psihat=psihat)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.