Случайные отрисовки из ANOVA-подобного дизайна с заданными размерами популяции

Question

Допустим, у вас есть нормально распределенная переменная y с категориальным предиктором * 3 *, который имеет ортогональные контрасты c1 и c2. Я пытаюсь создать программу на R, которая с учетом x, c1 и c2 создает y так, что c1 и c2 имеют размеры эффектов r1 и r2, указанные пользователем.

Например, предположим, что x, c1, c2, r1 и r2 были созданы следующим образом:

x <- factor(rep(c(1, 2, 3), 100))
contrasts(x) <- matrix(c(0, -.5, .5, -2/3, 1/3, 1/3), 
  nrow = 3, ncol = 2, dimnames = list(c("1", "2", "3"), c("c1", "c2")))

contrasts(x)
    c1         c2
1  0.0 -0.6666667
2 -0.5  0.3333333
3  0.5  0.3333333

r1 <- .09
r2 <- 0

Я бы хотел, чтобы программа создала y так, чтобы дисперсия y, учитываемая c1, равнялась r1 (.09), а дисперсия y, учитываемая c2, равнялась r2 (0).

Кто-нибудь знает, как я могу поступить об этом? Я знаю, что мне следует использовать функцию rnorm, но я застрял на том, какое количество населения должно использовать / sds rnorm, когда оно делает выборку.

Answer 1

Благодаря некоторым щедрым советам моих коллег, теперь у меня есть одна функция, которая создает смоделированные данные с заданным числом групп, набором контрастов, набором коэффициентов регрессии, указанным N на ячейку и указанным внутри- дисперсия группы

sim.factor <- function(levels, contr, beta, perCell, errorVar){
  # Build design matrix X
  X <- cbind(rep(1,levels*perCell), kronecker(contr, rep(1,perCell)))
  # Generate y
  y <- X %*% beta + rnorm(levels*perCell, sd=sqrt(errorVar))
  # Build and return data frame
  dat <- cbind.data.frame(y, X[,-1])
  names(dat)[-1] <- colnames(contr)
  return(dat)
}

Я также написал функцию, которая, учитывая набор коэффициентов регрессии, N на ячейку, количество групп, набор ортогональных контрастов, желаемую дельта-R ^ 2 для контраста интереса, возвращает требуемую дисперсию внутри группы:

ws.var <- function(levels, contr, beta, perCell, dc){
  # Build design matrix X
  X <- cbind(rep(1,levels), contr)
  # Generate the expected means
  means <- X %*% beta
  # Find the sum of squares due to each contrast 
  var <- (t(means) %*% contr)^2 / apply(contr^2 / perCell, 2, sum)
  # Calculate the within-conditions sum of squares
  wvar <- var[1] / dc - sum(var)
  # Convert the sum of squares to variance
  errorVar <- wvar / (3 * (perCell - 1))
  return(errorVar)
}

После некоторого тестирования следующим образом, кажется, что функции генерируют желаемую дельту R ^ 2 для контраста c1.

contr <- contr.helmert(3)
colnames(contr) <- c("c1","c2")
beta <- c(0, 1, 0)
perCell <- 50
levels = 3
dc <- .08
N <- 1000

# Calculate the error variance
errorVar <- ws.var(levels, contr, beta, perCell, dc)

# To store delta R^2 values
d1 <- vector("numeric", length = N)

# Use the functions
for(i in 1:N)
{
   d <- sim.factor(levels=3,
                   contr=contr,
                   beta=beta,
                   perCell=perCell,
                   errorVar=errorVar)
   d1[i] <- lm.sumSquares(lm(y ~ c1 + c2, data = d))[1, 2] # From the lmSupport package
}

m <- round(mean(d1), digits = 3)

bmp("Testing simulation functions.bmp")
hist(d1, xlab = "Percentage of variance due to c1", main = "")
text(.18, 180, labels = paste("Mean =", m))
dev.off()

Patrick