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