require(bindata) # Parameters of joint distribution. N <- 10 p1 <- 0.5 p2 <- 0.5 rho<- .6 hwe_diag = (2*freq_af$af * ( 1- freq_af$af) * freq_af$n) #one_hundred = diag(100) varx = 2*p1*(1-p1) vary = 2*p2*(1-p2) # Create one pair of correlated binomial values X = c() N=100 for(i in 1:N){ h1 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) h2 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) g1 = (h1[,1] + h2[,1]) g2 = (h1[,2] + h2[,2]) X = cbind(X, g1,g2) } X2 = t(apply(X,1,function(x){scale(x)})) h1 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) h2 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) g1 = (h1[,1] + h2[,1]) g2 = (h1[,2] + h2[,2]) h1 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) h2 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) g3 = (h1[,1] + h2[,1]) g4= (h1[,2] + h2[,2]) X = cbind(g1,g2,g3,g4) X = X -1 cor(g1,g2) beta = rep(c(0.1,0),100) pca =svd(cor(X2)) W = X2 %*% pca$v[,1:10] W2 = (pca$v[1:10,]) %*% X2 solve(t(W)%*%W) %*% t(W) %*% y # PCA regression from summary statistics, probably can't do this from SS. solve(t(W)%*%W) %*% (diag(diag((t(W)%*%W)))) %*% (diag(diag(ginv(t(W)%*%W))) %*% t(W) %*% y) solve(t(W)%*%W) %*% t(W) %*% y h_sq = .6 eta = X2 %*% beta y = eta + rnorm(N, 0, sd=sqrt((1-h_sq))) X_small = X[1:1e5,] y_small = y[1:1e5] summary(glmn(y ~ X)) summary(lm(y ~ X[,1])) summary(lm(y ~ X[,2])) summary(lm(y ~ X[,3])) summary(lm(y ~ X[,4])) summary(lm(y ~ X[,1:2])) summary(lm(residuals(lm(y ~ X[,1])) ~ -1 + X[,2])) N_Case = 20000 N_Control = 80000 hse <- X[y>qnorm(0.9,mean=0,sd=sd(y)),] Xcase <- Xcase[1:N_Case,] Xcontrol <- X[y<qnorm(0.9,mean=0,sd=sd(y)),] Xcontrol <- Xcontrol[1:N_Control,] X <- rbind(Xcase,Xcontrol) Ycase <- y[y>qnorm(0.9,mean=0,sd=sd(y))] Ycase <- Ycase[1:N_Case] Ycontrol <- y[y<qnorm(0.9,mean=0,sd=sd(y))] Ycontrol <- Ycontrol[1:N_Control] y = c(Ycase,Ycontrol) ycc<- c(rep(T,20000),rep(F,80000)) ld_matrix = cor(X_small) for(i in 1:ncol(X_small)){ g = (summary(lm(y ~ X_small[,i]))) beta = }
N=100 for(i in 1:100){ h1 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) h2 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) g1 = (h1[,1] + h2[,1]) g2 = (h1[,2] + h2[,2]) X = cbind(X, g1,g2) } h1 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) h2 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) g1 = (h1[,1] + h2[,1]) g2 = (h1[,2] + h2[,2]) h1 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) h2 <- rmvbin(N, c(p1,p2), bincorr=(1-rho)*diag(2)+rho) g3 = (h1[,1] + h2[,1]) g4= (h1[,2] + h2[,2]) X = cbind(g1,g2,g3,g4) X = X -1 cor(g1,g2) beta = rep(c(0.1,0),100) h_sq = 12 eta = X %*% beta y = eta + rnorm(N, 0, sd=sqrt((20-h_sq))) X_small = X[1:1e5,] y_small = y[1:1e5] summary(glmn(y ~ X)) #install.packages("glmnet") library(glmnet) library(MASS) (lm.ridge(y ~ X,lambda = 10000)) fit = glmnet(X, y, alpha = 0, weights = c(rep(1,50),rep(2,50)), nlambda = 20) cvfit = cv.glmnet(X, y, type.measure = "mse", nfolds = 20)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.