library(nlme) source("rmultnorm.R") effect <- .3 # effect size corr <- .4 # intrafamilial correlation numsim <- 1000 n1fam <- 50 # families with 2 patients and 1 relative n2fam <- 50 # families with 1 patient and 2 relatives vmat <- matrix(c(1, corr,corr, corr,1, corr, corr,corr, 1),3,3) x <- c(rep(1,n1fam),rep(1,n1fam),rep(0,n1fam), rep(1,n2fam),rep(0,n2fam),rep(0,n2fam)) id <- c(1:n1fam,1:n1fam,1:n1fam, (n1fam+1:n2fam),(n1fam+1:n2fam),(n1fam+1:n2fam)) power <- rep(0,numsim) for (i in 1:numsim) { cat(i," ") grp1 <- rmultnorm(n1fam, c(effect,effect, 0), vmat) grp2 <- rmultnorm(n2fam, c(effect,0,0), vmat) y <- c(grp1[,1],grp1[,2],grp1[,3], grp2[,1],grp2[,2],grp2[,3]) group <- groupedData(y~x|id) res <- lme(group,random = ~ 1) pval <- summary(res)$tTable[2,5] power[i] <- pval<=0.05 } cat("\nempirical power for effect size of ",effect, " is ",round(sum(power)/numsim,3),".\n",sep="")