2012-02-03 1 views
0

3つのグループのカテゴリ予測子xを持つ正規分布変数yを持ち、直交コントラストc1とc2を持つとします。私はRで、x、c1、およびc2が与えられたときに、c1とc2がユーザによって指定された効果サイズr1とr2を持つようなyを作成するプログラムを作成しようとしています。例えば与えられた人口効果の大きさを持つANOVAのようなデザインからランダムを引き出します

は、のは、その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)に等しく、c2によって占められるyの分散はr2(0)に等しい。

私はこれについてどのように考えているのですか?私はrnorm関数を使用すべきであることを知っていますが、サンプリングを行うときにどの人口が/ sds rnormが使用すべきであるかを知りません。あなたの助けを事前に

おかげで、

パトリック

答えて

0

同僚からいくつかの寛大なアドバイスの礼儀は、私は今のグループの指定された数の特定のシミュレートされたデータを作成する一つの機能、コントラストのセットを持って、回帰係数のセット、セルごとに指定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) 
} 

Iはまた、セット回帰係数のセット、セル当たりN、グループの数、与えられた、機能を書い直交コントラスト次のようにいくつかのテストを行った後

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) 
} 

、機能はコントラストC1に対する所望のデルタR^2を生成するように見える:STSは、必要と群内分散を返し、関心のコントラストのデルタ-R^2を所望しました。

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() 

パトリック

関連する問題