use*_*085 5 simulation statistics r genetics
我想生成2个连续随机变量Q1,Q2(数量性状,各自是正常的)和2个二进制随机变量Z1,Z2(二进制性状)与所有可能的对它们的之间给出成对的相关性.说
(Q1,Q2):0.23
(Q1,Z1):0.55
(Q1,Z2):0.45
(Q2,Z1):0.4
(Q2,Z2):0.5
(Z1,Z2):0.47
Run Code Online (Sandbox Code Playgroud)
请帮我在R中生成这样的数据
这很粗糙,但可能会让您朝着正确的方向开始。
library(copula)
options(digits=3)
probs <- c(0.5,0.5)
corrs <- c(0.23,0.55,0.45,0.4,0.5,0.47) ## lower triangle
Run Code Online (Sandbox Code Playgroud)
模拟相关值(前两个定量,后两个转换为二进制)
sim <- function(n,probs,corrs) {
tmp <- normalCopula( corrs, dim=4 , "un")
getSigma(tmp) ## test
x <- rCopula(1000, tmp)
x2 <- x
x2[,3:4] <- qbinom(x[,3:4],size=1,prob=rep(probs,each=nrow(x)))
x2
}
Run Code Online (Sandbox Code Playgroud)
测试观察到的相关性和目标相关性之间的 SSQ 距离:
objfun <- function(corrs,targetcorrs,probs,n=1000) {
cc <- try(cor(sim(n,probs,corrs)),silent=TRUE)
if (is(cc,"try-error")) return(NA)
sum((cc[lower.tri(cc)]-targetcorrs)^2)
}
Run Code Online (Sandbox Code Playgroud)
看看当输入 corrs=target 时情况有多糟糕:
cc0 <- cor(sim(1000,probs=probs,corrs=corrs))
cc0[lower.tri(cc0)]
corrs
objfun(corrs,corrs,probs=probs) ## 0.112
Run Code Online (Sandbox Code Playgroud)
现在尝试优化。
opt1 <- optim(fn=objfun,
par=corrs,
targetcorrs=corrs,probs=c(0.5,0.5))
opt1$value ## 0.0208
Run Code Online (Sandbox Code Playgroud)
501 次迭代后停止,并显示“超出最大迭代次数”。这永远不会很好地工作,因为我们试图在随机目标函数上使用确定性爬山算法......
cc1 <- cor(sim(1000,probs=c(0.5,0.5),corrs=opt1$par))
cc1[lower.tri(cc1)]
corrs
Run Code Online (Sandbox Code Playgroud)
也许尝试模拟退火?
opt2 <- optim(fn=objfun,
par=corrs,
targetcorrs=corrs,probs=c(0.5,0.5),
method="SANN")
Run Code Online (Sandbox Code Playgroud)
看起来并没有比之前的值好多少。两个可能的问题(留给读者作为练习)(1)我们指定了一组与我们选择的边际分布不可行的相关性,或者(2)目标函数曲面中的误差正在进入方式——为了做得更好,我们必须对更多的重复进行平均(即增加n)。