[R] Duncan's MRT: limitations to qtukey() function?

[Peter Dalgaard]

Andrea Onofri wrote:
> Dear all,
>
> I'm using R to perform multiple comparison testing on agriculture 
> genotype trials. To perform  the Duncan's MRT, I use the qtukey() 
> function with the following syntax:
>
> qtukey(p = ((1 - 0.05) ^ (pos - 1)), nmeans = pos, df = ni)
>
> I experience a strange behaviour when the number of means in the trial 
> and the number of residual degrees of freedom (ni) becomes high 
> (orientatively pos > 30 and ni > 60), i.e. the function returns NaN 
> together with the following message:
>
> convergence failed in 'qtukey'
>
> Running the code below should reproduce a  statistical table for 
> Duncan's MRT, such those reported on old books of agricultural 
> statistics (ex. Steel and Torrie, 1960).
>
> table <- matrix(0, 28, 14)
> pos <- c(10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 50, 100, 200)
> df <- c(2:20, seq(22, 30, by = 2), 40, 60, 100, 1000)
> for (i in 1:28) table [i, ] <- qtukey(p = ((1 - 0.05) ^ (pos - 1)), 
> nmeans = pos, df = df[i])
> colnames(table) <- pos
> rownames(table) <- df
>
> Indeed, I get several NaN and some odd figures.
>
> Am I doing something wrong or is there any intrinsic limit for the 
> values of p, nmeans and df within the qtukey() function? Just in this 
> latter case, is it possible to raise those limits? Some genotype trials 
> can have more than thirty treatments  under comparison.
>
> I thank you all for your attention.
>
> Andrea Onofri
> Department of Agriculture and Environmental Sciences
> University of Perugia (ITALY)
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>   
Yes, this looks a bit suspicious, especially the fact that it worsens
with increasing df. The core of the matter seems to be this

> qtukey(log(3.689754e-05),200,1000,log=T)
[1] NaN
Warning messages:
1: In qtukey(log(3.689754e-05), 200, 1000, log = T) :
  convergence failed in 'qtukey'
2: In qtukey(p, nranges, nmeans, df, lower.tail, log.p) : NaNs produced

which presumably is from inverting ptukey, but

> f <- Vectorize(function(x)ptukey(x,200,1000,log=T)-log(3.689754e-05))
> curve(f, from=3,to=5)
> uniroot(f, c(3,5))
$root
[1] 3.760707

$f.root
[1] -5.756398e-05

$iter
[1] 6

$estim.prec
[1] 6.103516e-05

is looking perfectly nice.

-- 
   O__  ---- Peter Dalgaard             Øster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907