[R] why this function does not run correctly?

[Peter Ehlers]

On 2010-03-08 14:45, Marco Bressan wrote:
> Hi,
> my name is Marco Bressan i'm working to improve ADati package. I study psicology ad Padua University (Italy). I have this problem: why bartlett.test function running good and my anova.welch function no?
>
> Ciao,
> il mio nome � Marco Bressan e sto lavorando per migliorare il pacchetto ADati. Studio psicologia all'universit� di Padova. Non capisco come mai la funzione che ho fatto mi dia quell'errore
>
> this is anova.welch:
>
> anova.welch<- function(x, ...) UseMethod("anova.welch")
>
> anova.welch.default<-                                                        ## this is the algoritm, I think it's ok (I copy this from Welch() inside ADati)
> function (x, y = NULL, nu = c(0,0) ,...)
> {
>      mx<- tapply(x,y,mean)
>      s2x<- tapply(x,y,var)
>      k<- length(nu)
>      w<- nu/s2x
>      Xp<- sum(w * mx)/sum(w)
>      Fnum<- sum(w * (mx - Xp)^2)/(k - 1)
>      Fden<- 0
>      for (h in 1:k) {
>          a<- 1/(nu[h] - 1)
>          b<- (1 - (w[h]/sum(w)))^2
>          Fden<- Fden + (a * b)
>      }
>      gl2.den<- Fden * 3
>      Fden<- Fden * ((2 * (k - 2))/(k^2 - 1)) + 1
>      Fw<- Fnum/Fden
>      STATISTIC<- Fw
>      gl1<- k - 1
>      gl2.num<- k^2 - 1
>      gl2<- gl2.num/gl2.den
>      PARAMETER<- c(gl1, gl2)
>      PVAL<- pf(Fw, gl1, gl2, lower.tail = FALSE)
>      METHOD<- "Welch ANOVA"
>      DNAME<- NA
>      names(STATISTIC)<- "F"
>      names(PARAMETER)<- c("num df", "denom df")
>      RVAL<- list(statistic = STATISTIC, parameter = PARAMETER,
>          p.value = PVAL, method = METHOD)
>      class(RVAL)<- "htest"
>      return(RVAL)
> }
>
> anova.welch.formula<-                                                            ## I copy this from bartlett.test
> function(formula, data, subset, na.action, ...)
> {
>      if(missing(formula) || (length(formula) != 3L))
>          stop("'formula' mancante o incorretta")
>      m<- match.call(expand.dots = FALSE)
>      if(is.matrix(eval(m$data, parent.frame())))
>          m$data<- as.data.frame(data)
>      m[[1L]]<- as.name("model.frame")
>      mf<- eval(m, parent.frame())
>      DNAME<- paste(names(mf), collapse = " by ")
>      names(mf)<- NULL
>      y<- do.call("anova.welch", as.list(mf))
>      y$data.name<- DNAME
>      y
> }

bartlett.test doesn't have a 'nu' argument. Try this:

anova.welch.formula <-
function(formula, data, subset, na.action, nu = c(0,0), ...)
{
     if(missing(formula) || (length(formula) != 3L))
         stop("'formula' mancante o incorretta")
     m <- match.call(expand.dots = FALSE)
     m[["nu"]] <- NULL  ## needed so that eval(m,....) will work
     if(is.matrix(eval(m$data, parent.frame())))
         m$data <- as.data.frame(data)
     m[[1L]] <- as.name("model.frame")
     mf <- eval(m, parent.frame())
     DNAME <- paste(names(mf), collapse = " by ")
     names(mf) <- NULL
     y <- do.call("anova.welch", c(as.list(mf), list(nu)))
              ## include 'nu' in the parameters passed to
              ##    anova.welch.default
     y$data.name <- DNAME
     y
}

(You might also find the code for lm() instructive.)

I assume that you're aware of the oneway.test() function in
package:stats. So why re-invent the wheel?

  -Peter Ehlers

>
>> n1=10                                                                                  ## this is the test
>> n2=15
>> n3=20
>> y=c(rnorm((n1+n2),5,2),rnorm(n3,7,8))
>> A=factor(c(rep(1,n1),rep(2,n2),rep(3,n3)))
>> anova.welch(y,A,c(n1,n2,n3))
>
>          Welch ANOVA
>
> data:
> F = 2.3025, num df = 2.000, denom df = 27.384, p-value = 0.1191
>
>> anova.welch(y~A,nu=c(n1,n2,n3))
> Errore in model.frame.default(formula = y ~ A, ... = list(nu = c(n1, n2,  :
>    invalid type (pairlist) for variable '(...)'
>
>
> Sorry for my english, tanks you if you can help me :)
> 	[[alternative HTML version deleted]]
>
>
>
>
> ______________________________________________
> 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.

-- 
Peter Ehlers
University of Calgary