Re: [R] p-values

From: Peter Ho <peter_at_estg.ipvc.pt>
Date: Thu 11 Aug 2005 - 23:28:00 EST

Spencer,

Here is an example from rayner and best 2001 and the script sent by Felipe. This can be done as follows using the function durbin.grupos() in the attached file

 > ###Ice cream example from Rayner and Best 2001 . Chapter 7
 > judge <- rep(c(1:7),rep(3,7))
 > variety <- c(1,2,4,2,3,5,3,4,6,4,5,7,1,5,6,2,6,7,1,3,7)
 > cream <- c(2,3,1,3,1,2,2,1,3,1,2,3,3,1,2,3,1,2,3,1,2)
 > durbin.grupos(judge,variety,cream,k=3,r=3,alpha=0.01)

Prueba de Durbin
..............

Chi Cuadrado :  12
Gl.          :  6
P-valor      :  0.0619688

..............
Comparación de tratamientos
Alpha        :  0.01
Gl.          :  8
t-Student    :  3.355387

Diferencia minima
para la diferencia entre suma de rangos = 4.109493

Grupos, Tratamientos y la Suma de sus rangos

a        2       9
ab       1       8
abc      7       7
abc      6       6
abc      5       5
 bc      3       4
  c      4       3

  trat prom M
1 2 9 a
2 1 8 ab
3 7 7 abc
4 6 6 abc
5 5 5 abc
6 3 4 bc
7 4 3 c
 >

You can see that the p-value is the same with

 > pchisq(12, df= 6, lower.tail=F)
[1] 0.0619688

I am hoping that someone, maybe Torsten, might be able to suggest how I can incorporate Monte-Carlo p-values using pperm(). The statistical issues are beyond my comprehension and I assume that Rayner and Best suggestion to use Monte-Carlo p-values instead of Chi-square p-values to be correct. In the above example the Monte-Carlo p-value is 0.02. This is a significant difference, resulting in the rejection of the null hypothesis when using Monte-Carlo p-values.

I hope this example might help. Thanks again for your answer and also to Felipe for sending the function for Durbin's test.

Peter

Spencer Graves wrote:

> pperm seems reasonable, though I have not looked at the details.
>
> We should be careful about terminology, however. So-called
> "exact p-values" are generally p-values computed assuming a
> distribution over a finite set of possible outcomes assuming some
> constraints to make the outcome space finite. For example, Fisher's
> exact test for a 2x2 table assumes the marginals are fixed.
>
> I don't remember the details now, but I believe there is
> literature claiming that this may not be the best thing to do when,
> for example, when it is reasonable to assume that the number in each
> cell is Poisson. In such cases, you may lose statistical power by
> conditioning on the marginals. I hope someone else will enlighten us
> both, because I'm not current on the literature in this area.
>
> The situation with "exact tests" and "exact p-values" is not
> nearly as bad as with so-called ""exact confidence intervals", which
> promise to deliver at least the indicated coverage probability. With
> discrete distributions, it is known that 'Approximate is better than
> "exact' for interval estimation of binomial proportions', as noted in
> an article of this title by A. Agresti and B. A. Coull (1998) American
> Statistician, 52: 119-126. (For more on this particular issue, see
> Brown, Cai and Dasgupta 2003 "Interval Estimation in Exponential
> Families", Statistica Sinica 13: 19-49).
>
> If this does not answer your question adequately, may I suggest
> you try the posting guide. People report having found answers to
> difficult questions in the process of preparing a question following
> that guide, and when they do post a question, they are much more
> likely to get a useful reply.
>
> spencer graves
>
> Peter Ho wrote:
>
>> Spencer,
>>
>>
>> Thank you for referring me to your other email on Exact
>> goodness-of-fit test. However, I'm not entirely sure if what you
>> mentioned is the same for my case. I'm not a statistician and it
>> would help me if you could explain what you meant in a little more
>> detail. Perhaps I need to explain the problem in more detail.
>>
>> I am looking for a way to calculate exaxt p-values by Monte Carlo
>> Simulation for Durbin's test. Durbin's test statistic is similar to
>> Friedman's statistic, but considers the case of Balanced Incomplete
>> block designs. I have found a function written by Felipe de Mendiburu
>> for calculating Durbin's statistic, which gives the chi-squared
>> p-value. I have also been read an article by Torsten Hothorn "On
>> exact rank Tests in R" (R News 1(1), 11–12.) and he has shown how to
>> calculate Monte Carlo p-values using pperm. In the article by Torsten
>> Hothorn he gives:
>>
>> R> pperm(W, ranks, length(x))
>>
>> He compares his method to that of StatXact, which is the program
>> Rayner and Best suggested using. Is there a way to do this for
>> example for the friedman test.
>>
>> A paper by Joachim Rohmel discusses "The permutation distribution for
>> the friendman test" (Computational Statistics & Data Analysis 1997,
>> 26: 83-99). This seems to be on the lines of what I need, although I
>> am not quite sure. Has anyone tried to recode his APL program for R?
>>
>> I have tried a number of things, all unsucessful. Searching through
>> previous postings have not been very successful either. It seems that
>> pperm is the way to go, but I would need help from someone on this.
>>
>> Any hints on how to continue would be much appreciated.
>>
>>
>> Peter
>>
>>
>> Spencer Graves wrote:
>>
>>> Hi, Peter:
>>>
>>> Please see my reply of a few minutes ago subject: exact
>>> goodness-of-fit test. I don't know Rayner and Best, but the same
>>> method, I think, should apply. spencer graves
>>>
>>> Peter Ho wrote:
>>>
>>>
>>>
>>>> HI R-users,
>>>>
>>>> I am trying to repeat an example from Rayner and Best "A
>>>> contingency table approach to nonparametric testing (Chapter 7, Ice
>>>> cream example).
>>>>
>>>> In their book they calculate Durbin's statistic, D1, a dispersion
>>>> statistics, D2, and a residual. P-values for each statistic is

>>>> calculated from a chi-square distribution and also Monte Carlo
>>>> p-values.
>>>>
>>>> I have found similar p-values based on the chi-square distribution
>>>> by using:
>>>>
>>>> > pchisq(12, df= 6, lower.tail=F)
>>>> [1] 0.0619688
>>>> > pchisq(5.1, df= 6, lower.tail=F)
>>>> [1] 0.5310529
>>>>
>>>> Is there a way to calculate the equivalent Monte Carlo p-values?
>>>>
>>>> The values were 0.02 and 0.138 respectively.
>>>>
>>>> The use of the approximate chi-square probabilities for Durbin's
>>>> test are considered not good enough according to Van der Laan (The
>>>> American Statistician 1988,42,165-166).
>>>>
>>>>
>>>> Peter
>>>> --------------------------------
>>>> ESTG-IPVC
>>>>
>>>> ______________________________________________
>>>> R-help@stat.math.ethz.ch mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>>> PLEASE do read the posting guide!
>>>> http://www.R-project.org/posting-guide.html
>>>>
>>>
>>>
>>>
>>>
>>>
>>
>

#---------------------------------
# Prueba de Durbin para bloques incompletos
# k = nro de unidades experimentales por bloque
# r = numero de veces que aparece el tratamiento # Lambda = numero de bloques en el que dos tratamientos estan juntos (diseño) durbin.comp<-function(juez,trat,eval,k,r,alpha=0.05,main=" ") {
x<-data.frame(juez,trat,eval)
t<-length(unique(x$trat))
b<-length(unique(x$juez))

# Determina el rango dentro de cada juez z <- by(x,x$juez,function(x) rank(x$eval))
y<-data.frame(c(z))
m<-dim(y)
n<-m[1]*m[2]

rango <- 1:n
for (i in 1:m[1]) {
for (j in 1:m[2]) {
kk=i+m[1]*(j-1)
rango[kk]<-y[i,j]
}
}
x<-data.frame(x,rango)
z <-by(x,x$trat,function(x) sum(x$rango)) y<-as.vector(c(z))
s <- (y-r*(k+1)/2)^2
s1 <- sum(s)
# determina el valor de Durbin
gl1<-t-1 ;gl2<-b*k-t-b+1
s <- 12*(t-1)*s1/(r*t*(k-1)*(k+1))
prob<-1-pchisq(s,gl1); Tprob<-qt(1-alpha/2,gl2) Mc <-Tprob*sqrt( r*k*(k+1)*(b*(k-1)-s)/(6*gl2))

# s,prob,Tprob,Mc,gl1,gl2)
# Impresion de resultados

cat("\nEstudio:",main,"\n")
cat("Prueba de Durbin\n\n")
cat("Valor calculado    = ",s,"\n")
cat("Grados de libertad = ",gl1,"\n")
cat("P-valor            = ",prob,"\n\n")
cat("t-Student          = ",Tprob,"\n")
cat("Grados de libertad = ",gl2,"\n\n")
cat("Diferencia minima\npara la diferencia entre suma de rangos = ",Mc,"\n\n")

# comparacion de tratamientos.
comb <-combn(t,2)
n<-ncol(comb)
diff<-array(0,n)
stat<-array(" ns ",n)
for (kk in 1:n) {
diff[kk]<-abs(y[comb[1,kk]]-y[comb[2,kk]]) if (diff[kk] >= Mc) stat[kk]<-" * "
}
tcomb<-t(comb)
tstat<-cbind(stat)
tr.comp<-cbind(diff)
cat("Suma de rangos\n\Trat Suma\n----------\n") for(i in 1:t){
a <-c(i)
b <-c(" = ",y[i])
cat(a,b,"\n")
}
compara <- data.frame(tcomb,tr.comp,tstat) return(compara)
}

#---------------------------------
# Prueba de Durbin para bloques incompletos
# k = nro de unidades experimentales por bloque
# r = numero de veces que aparece el tratamiento # Lambda = numero de bloques en el que dos tratamientos estan juntos (diseño) durbin.grupos<-function(juez,trat,eval,k,r,alpha=0.05, main=" ") { name.y <- paste(deparse(substitute(eval)))
x<-data.frame(juez,trat,eval)
t<-length(unique(x$trat))
b<-length(unique(x$juez))

# Determina el rango dentro de cada juez z <- by(x,x$juez,function(x) rank(x$eval))
y<-data.frame(c(z))
m<-dim(y)
n<-m[1]*m[2]

rango <- 1:n
for (i in 1:m[1]) {
for (j in 1:m[2]) {
kk=i+m[1]*(j-1)
rango[kk]<-y[i,j]
}
}
x<-data.frame(x,rango)
z <-by(x,x$trat,function(x) sum(x$rango)) y<-as.vector(c(z))
nombre<-as.character(dimnames(z)$"x$trat") s <- (y-r*(k+1)/2)^2
s1 <- sum(s)
# determina el valor de Durbin
gl1<-t-1 ;gl2<-b*k-t-b+1
s <- 12*(t-1)*s1/(r*t*(k-1)*(k+1))
prob<-1-pchisq(s,gl1); Tprob<-qt(1-alpha/2,gl2) Mc <-Tprob*sqrt( r*k*(k+1)*(b*(k-1)-s)/(6*gl2))

# s,prob,Tprob,Mc,gl1,gl2)
# Impresion de resultados

cat("\nEstudio:",main)
cat("\nPrueba de Durbin")
cat("\n..............")
cat("\nChi Cuadrado : ",s)
cat("\nGl.          : ",gl1)
cat("\nP-valor      : ",prob)
cat("\n..............")
cat("\nComparación de tratamientos\n")
cat("\nAlpha        : ",alpha)
cat("\nGl.          : ",gl2)
cat("\nt-Student    : ",Tprob)
cat("\nDiferencia minima\npara la diferencia entre suma de rangos = ",Mc)
# comparacion de tratamientos.
cat("\n\nGrupos, Tratamientos y la Suma de sus rangos\n") y<-as.numeric(y)
resultado<-orden.est(nombre,y,Mc)
#
return(resultado)
}
#-----------------
# establece el orden estadistico en la comparacion de tratamientos
# Nombre de tratamientos : trat

# promedios : prom
# valor minimo de comparacion

orden.est<-function(tratamiento,promedio,minimo) { n<-length(promedio)
z<-data.frame(tratamiento,promedio)
# ordena tratamientos

w<-z[order(z[,2],decreasing=T),]
M<-rep("",n)
k<-1

k1<-0
j<-1
i<-1
cambio<-n
cambio1<-0
chequeo=0
M[1]<-letters[k]
while(j<n) {
chequeo<-chequeo+1
if (chequeo > n) break
for(i in j:n) {
s<-abs(w[i,2]-w[j,2])<=minimo
if(s) {
if(ultimo.c(M[i]) != letters[k])M[i]<-paste(M[i],letters[k],sep="") }
else {
k<-k+1
cambio<-i
cambio1<-0
ja<-j
for(jj in cambio:n) M[jj]<-paste(M[jj]," ",sep="") M[cambio]<-paste(M[cambio],letters[k],sep="") for( v in ja:cambio) {
if(abs(w[v,2]-w[cambio,2])>minimo) {j<-j+1 cambio1<-1
}
else break
}
break
}
}
if (cambio1 ==0 )j<-j+1
}
#-----------
w<-data.frame(w,stat=M)
trat<-as.character(w$tratamiento)
prom<-as.numeric(w$promedio)
for(i in 1:n){
cat(M[i],"\t",trat[i],"\t",prom[i],"\n") }
resultado<-data.frame(trat,prom,M)
return(resultado)
}
ultimo.c<-function(x) {
y<-sub(" +$", "",x)
p1<-nchar(y)
cc<-substr(y,p1,p1)
return(cc)
}
##



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Aug 11 23:31:57 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 15:14:42 EST