Я в некотором роде новичок в R. Я пытаюсь понять формат функции basenonpartest, написанной в Bathke et al 2008. Я понимаю, что фрейм должен быть тем фреймом данных, который я хочу ввести, в моем случае это sberry'набор данных, который является общедоступным.Я также понимаю, что, поскольку я хочу, чтобы он выполнял тест ANOVA, мне нужны тесты = 1.Однако я не уверен, что нужно вводить для группы и переменных.Спасибо за помощь.
basenonpartest <-
function(frame,group,vars,tests=c(1,1,1,1)){
if(!is.factor(frame[,group]))
{
frame[,group] <- factor(frame[,group])
}
# Orders data by group
o <- order(frame[,group])
frame <- frame[o,]
# Compute sample sizes per group
p <- length(vars)
a <- length(levels(frame[,group]))
N <- length(frame[,1])
ssize <- array(NA,a)
lims <- matrix(NA,2,a)
for(i in 1:a){
ssize[i] <- length(frame[frame[,group]==levels(frame[,group])[i],1])
lims[1,i] <- min(which(frame[,group]==levels(frame[,group])[i]))
lims[2,i] <- max(which(frame[,group]==levels(frame[,group])[i]))
}
# Sets up R matrix
Rmat <- matrix(NA,N,p)
for(j in 1:p){
Rmat[,j] <- rank(frame[,vars[j]],ties.method="average")
}
# Manipulating R
Rbars <- matrix(NA,a,p)
for(i in 1:a){
for(j in 1:p){
Rbars[i,j] <- mean(Rmat[(lims[1,i]:lims[2,i]),j])
}
}
Rtilda <- (1/a)*colSums(Rbars)
Rbarovr <- (1/N)*colSums(Rmat)
H1 <- matrix(0,p,p)
H2 <- matrix(0,p,p)
G1 <- matrix(0,p,p)
G2 <- matrix(0,p,p)
G3 <- matrix(0,p,p)
for(i in 1:a){
H1 <- H1 + ssize[i]*(Rbars[i,] - Rbarovr)%*%t(Rbars[i,] -Rbarovr)
H2 <- H2 + (Rbars[i,] - Rtilda)%*%t(Rbars[i,] - Rtilda)
for(j in 1:ssize[i]){
G1 <- G1 + (((Rmat[(lims[1,i]+j-1),(1:p)])-
Rbars[i,])%*%t((Rmat[(lims[1,i]+j-1),(1:p)])-Rbars[i,]))
G2 <- G2 + (1-(ssize[i]/N))*(1/(ssize[i]-1))*(((Rmat[(lims[1,i]+j-1),
(1:p)])-Rbars[i,])%*%t((Rmat[(lims[1,i]+j-1),(1:p)])-Rbars[i,]))
G3 <- G3 + (1/(ssize[i]*(ssize[i]-1)))*(((Rmat[(lims[1,i]+j-1),(1:p)])-
Rbars[i,])%*%t((Rmat[(lims[1,i]+j-1),(1:p)])-Rbars[i,]))
}
}
H1 <- (1/(a-1))*H1
H2 <- (1/(a-1))*H2
G1 <- (1/(N-a))*G1
G2 <- (1/(a-1))*G2
G3 <- (1/a)*G3
# ANOVA TYPE TEST
fhat <- ((a-1)*sum(diag(G3))^2)/sum(diag(G3%*%G3))
fhat0 <- ((a^2)/((a-1)*sum(1/(ssize-1))))*fhat
Fanova <- sum(diag(H2))/sum(diag(G3))
pvalanova <- 1 - pf(Fanova,fhat,fhat0)
df1anova<-fhat
df2anova<-fhat