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@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
______________________________________________
R-help@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.
