Как обновить код, чтобы создать функцию для расчета Уэлча для полиномиальных трендов? - PullRequest
Новая
       35

Как обновить код, чтобы создать функцию для расчета Уэлча для полиномиальных трендов?

3 голосов
/ 19 июня 2019

Я пытаюсь воспроизвести вывод SPSS для значимости линейного тренда среди средних значений, когда равные отклонения не предполагаются. Я с благодарностью использовал код из http://www -personal.umich.edu / ~ gonzo / coursenotes / file3.pdf , чтобы создать функцию для вычисления отдельных дисперсий, которую на основе моего поиска я понимаю как «равный отклонения не принимаются »выводится в SPSS.

Моя проблема / цель: Я оцениваю только полиномиальные ортогональные тренды (в основном линейные). Я хочу адаптировать код, создающий функцию, чтобы аргумент контраста мог принимать заранее созданные матрицы контраста, а не указывать коэффициенты каждый раз вручную (место для опечаток!).

… Я попробовал эти точные команды, но получил Error in contrast %*% means : non-conformable arguments. Я поиграл с кодом, но не могу заставить его работать.

Код для создания функции из заметок: 

sepvarcontrast <- function(dv, group, contrast) {
  means <- c(by(dv, group, mean))
  vars <- c(by(dv, group, var))
  ns <- c(by(dv, group, length))
  ihat <- contrast %*% means
  t.denominator <- sqrt(contrast^2 %*% (vars/ns))
  t.welch <- ihat/ t.denominator
  num.contrast <- ifelse(is.null(dim(contrast)),1,dim(contrast)[1])
  df.welch <- rep(0, num.contrast)
  if (is.null(dim(contrast))) contrast <- t(as.matrix(contrast))
  for (i in 1:num.contrast) {
    num <- (contrast[i,]^2 %*% (vars))^2
    den <- sum((contrast[i,]^2 * vars)^2 / (ns-1))
    df.welch[i] <- num/den
    }
  p.welch <- 2*(1- pt(abs(t.welch), df.welch))
  result <- list(ihat = ihat, se.ihat = t.denominator, t.welch = t.welch,
                 df.welch = df.welch, p.welch = p.welch)
  return(result)
}

Я бы хотел использовать такую ​​функцию: 

# Create a polynomial contrast matrix for 5 groups, then save 

contr.mat5 <- contr.poly(5)


# Calculate separate variance

sepvarcontrast(dv, group, contrast = contr.mat5)

Я попробовал эти точные команды, чтобы увидеть, будут ли они работать, но получил Error in contrast %*% means : non-conformable arguments.

Все предложения приветствуются! Я все еще учусь создавать репер ...

...