Вопрос о функции basenonparest в пакете npmv? - PullRequest
0 голосов
/ 08 декабря 2018

Я в некотором роде новичок в 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
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...