Nothing
m_ng <- function(y, data, FS, medstar = c(0.01,0.0001), numb = 100, burnin = 1, every = 1 ) {
g1 <- which( FS > 0.5 ) ##### Group1: FS > 0.5
g2 <- which( FS == 0.5 ) ##### Group2: FS = 0.5
g3 <- which( FS < 0.5 ) ##### Group3: FS < 0.5
g4 <- which( is.na( FS ) ) ##### Group4: FS = NA
group <- list( g1, g2, g3, g4 ) ##### Create list for the 4 Groups
n.g <- length( group )
len_group <- rep( 0, n.g )
for ( i in 1:n.g ) len_group[i] = length( group[[i]] ) #### number of SNPs in each Group
x <- Rfast::standardise( as.matrix( data ) )
x <- model.matrix(~x)
#mod <- glm( y ~ x, binomial )
#betas <- as.vector( mod$coefficients )
#covs <- vcov( mod )
#solvecovs <- solve( covs, tol = 1e-30 )
dm <- dim(data)
n <- dm[1] ; p <- dm[2]
newlambdastar <- lambdastar <- lambdasd <- lambdaaccept <- lambdacount <- accept <- logaccept <- sha <- gammasq <- newgammasq <- rep( 1, n.g )
newlambda <- list(999)
lambda <- list( rep( lambdastar[1], len_group[1] ), rep( lambdastar[2], len_group[2] ),
rep( lambdastar[3], len_group[3] ), rep( lambdastar[4], len_group[4] ) )
LAMBDA <- GAMMASQ <- rep( 1, p )
w <- 1
h <- 0
sigmasq <- 1
psi <- rep( 2 * lambda[[1]][1] * gammasq[1] * 0.01, p )
numbofits <- burnin + every * numb ## total number of iterations
holdpsi <- matrix( rep( 0, p * numb ), ncol = numb )
holdbeta <- matrix( rep( 0, p * numb ), ncol = numb )
holdgammasq <- matrix( rep( 0, n.g * numb ), ncol = numb )
holdlambda <- matrix( rep( 0, n.g * numb ), ncol = numb )
holdLAMacsept <- matrix( rep( 0, n.g * numb ), ncol = numb )
holdalpha <- numeric( numb )
holdsigmasq <- numeric( numb )
holdW <- numeric( numb )
holdH <- numeric( numb )
lambdasd <- rep( 0.01, n.g )
#const <- solvecovs %*% betas
crosspord.x <- crossprod(x)
crossprod.xy <- crossprod(x,y)
################# Updating beta and alpha #############
for (i in 1:numbofits ){
LAM <- diag( c( 0, 1 / psi ) )
#LAM[LAM>1e+7]=1e+7
#LAM[LAM<1e-7 & LAM!=0]=1e-7
sigLAM <- sigmasq * LAM
#covsLAM[covsLAM <1e-7 & covsLAM!=0]=1e-7
varstar <- sigmasq * solve( crosspord.x + sigLAM, tol = 1e-30 )
expec <- solve( crosspord.x + sigLAM, tol = 1e-30 ) %*% crossprod.xy
cholstar <- chol( varstar, pivot = TRUE )
if ( attr( cholstar, "rank" ) == ncol( cholstar ) ) {
pivot <- attr( cholstar, "pivot" )
cholstar <- cholstar[ ,order( pivot ) ]
randn <- rnorm( p + 1 )
be <- expec + crossprod( cholstar, randn )
alpha <- be[1]
beta <- be[-1]
}
################### Updating sigmasq #######################
c.star <- n / 2
dstard <- y - alpha - data %*% beta
d.star <- crossprod( dstard )
sigmasq <- 1 / rgamma( 1, c.star, scale = 1 / d.star )
################### Updating Psi #######
for (j in 1:p) {
### main check 1
if ( beta[j]^2 < 10^( -5 ) ) {
check <- 0
### subcheck 1 for lambda
if ( LAMBDA[j] < 0.5 ) {
while ( check == 0 ) {
psi[j] <- 1 / rgamma( 1, 0.5 - LAMBDA[j], scale = 1 / ( 0.5 * beta[j]^2 ) )
u <- runif( 1 )
check <- as.numeric( u < exp(- psi[j] / ( 2 * GAMMASQ[j] ) ) ) }
} else {
while ( check == 0) {
psi[j] <- rgamma( 1, LAMBDA[j] - 0.5, scale = 2 * GAMMASQ[j] )
u <- runif( 1 )
check <- as.numeric( u < exp( - 0.5 * beta[j]^2 / psi[j] ) ) } }
} else { psi[j] <- rgig( 1, LAMBDA[j] - 0.5, beta[j]^2, 1 / GAMMASQ[j] ) }
} ### end for (j in 1:p )
#psi[psi>1e+4]=1e+4
psi[ psi < 1e-10 & psi != 0 ] <- 1e-10
########################################## update \lambda for each Group #######
for ( g in 1:n.g ) {
mupsi <- 2 * lambdastar[g] * gammasq[g]
newlambdastar[g] <- lambdastar[g] * exp( lambdasd[g] *rnorm( 1 ) )
newgammasq[g] <- mupsi / ( 2 * newlambdastar[g] )
newlambda[[g]] <- rep( newlambdastar[g], len_group[g] )
logaccept[g] <- log( newlambdastar[g] ) - log( lambdastar[g] ) - 142.85 * ( newlambdastar[g] - lambdastar[g] )
logaccept[g] <- logaccept[g] - len_group[g] * newlambdastar[g] * log( 2 * newgammasq[g] ) - len_group[g] * lgamma( newlambdastar[g] )
logaccept[g] <- logaccept[g] + len_group[g] * lambdastar[g] * log( 2 * gammasq[g] ) + len_group[g] * lgamma( lambdastar[g] )
logaccept[g] <- logaccept[g] + newlambdastar[g] * sum( log( psi[group[[g]]] ) ) - sum( psi[group[[g]]] ) / ( 2 * newgammasq[g] )
logaccept[g] <- logaccept[g] - lambdastar[g] * sum( log( psi[group[[g]]] ) ) +sum( psi[group[[g]]] ) / ( 2 * gammasq[g] )
accept[g] <- 1
if ( logaccept[g] < 0 ) { accept[g] <- exp( logaccept[g] ) }
lambdasd[g] <- lambdasd[g] + ( accept[g] - 0.3) / i
lambdaaccept[g] <- lambdaaccept[g] + accept[g]
lambdacount[g] <- lambdacount[g] + 1
u <- runif( 1 )
if ( u < accept[g] ) {
lambda[[g]] <- newlambda[[g]]
lambdastar[g] <- newlambdastar[g]
gammasq[g] <- newgammasq[g]
}
################### Updating gamma^2 for each Group #######
if (g == 1) {
sha[g] <- 0.5 * sum( psi[group[[g]]] ) + medstar[1] / ( 2 * lambdastar[g] )
gammasq[g] <- 1 / rgamma( 1, sum( unlist( lambda[[g]] ) ) + 2, scale = 1 / sha[g] )
} else if (g == 2) {
A <- 0.5 * sum( psi[group[[g]]] ) + medstar[1] / ( 2 * lambdastar[g] )
B <- 0.5 * sum( psi[group[[g]]] ) + medstar[2] / ( 2 * lambdastar[g] )
delta1 <- ( w * ( medstar[1])^2 ) / ( ( w * ( medstar[1])^2 ) + ( ( A / B )^( len_group[g] * lambdastar[g] + 2 ) * ( 1 - w ) * ( medstar[2])^2 ) )
U <- runif( 1 )
if ( U < delta1 ) {
sha[g] <- 0.5 * sum( psi[group[[g]]] ) + medstar[1] / ( 2 * lambdastar[g] )
gammasq[g] <- 1 / rgamma( 1, sum( unlist( lambda[[g]] ) ) + 2, scale = 1 / sha[g] )
} else {
sha[g] <- 0.5 * sum( psi[group[[g]]] ) + medstar[2] / ( 2 * lambdastar[g] )
gammasq[g] <- 1 / rgamma( 1, sum( unlist( lambda[[g]] ) ) + 2, scale = 1 / sha[g] )
}
} else if (g == 3) {
sha[g] <- 0.5 * sum( psi[group[[g]]] ) + medstar[2] / ( 2 * lambdastar[g] )
gammasq[g] <- 1 / rgamma( 1, sum( unlist( lambda[[g]] ) ) + 2, scale = 1 / sha[g] )
} else {
A <- 0.5 * sum( psi[group[[g]]] ) + medstar[1] / ( 2 * lambdastar[g] )
B <- 0.5 * sum( psi[group[[g]]] ) + medstar[2] / ( 2 * lambdastar[g] )
delta2 <- ( h * ( medstar[1])^2 ) / ( ( h * (medstar[1])^2 ) + ( ( A / B )^(len_group[g] * lambdastar[g] + 2) * ( 1 - h ) * ( medstar[2])^2 ) )
U <- runif( 1 )
if (U < delta2) {
sha[g] <- 0.5 * sum( psi[group[[g]]] ) + medstar[1] / ( 2 * lambdastar[g] )
gammasq[g] <- 1 / rgamma( 1, sum( unlist( lambda[[g]] ) ) + 2, scale = 1 / sha[g] )
} else {
sha[g] <- 0.5 * sum( psi[group[[g]]] ) + medstar[2] / ( 2 * lambdastar[g] )
gammasq[g] <- 1 / rgamma( 1, sum( unlist( lambda[[g]] ) ) + 2, scale = 1 / sha[g] )
}
}
}
LAMBDA <- unlist( lambda )[order( unlist( group ) )]
GAMMASQ <- c( rep( gammasq[1], len_group[1] ),rep( gammasq[2], len_group[2] ), rep(gammasq[3], len_group[3] ), rep( gammasq[4], len_group[4] ) )
GAMMASQ <- GAMMASQ[ order( unlist( group ) ) ]
################### Updating w for each Group #######
U <- runif( 1 )
sha1 <- medstar[1] / ( 2 * lambdastar[2] )
sha2 <- medstar[2] / ( 2 * lambdastar[2] )
A <- dgamma( gammasq[2], 2, scale = 1 / ( sha1 ) )
B <- dgamma( gammasq[2], 2, scale = 1 / ( sha2 ) )
delta3 <- A / ( A + B )
if ( U < delta3 ) {
w <- rbeta( 1, 3, 2 )
} else {
w <- rbeta( 1, 2, 3 )
}
################### Updating h for each Group #######
U <- runif( 1 )
sha1 <- medstar[1] / ( 2 * lambdastar[4] )
sha2 <- medstar[2] / ( 2 * lambdastar[4] )
A <- dgamma( gammasq[4], 2, scale = 1 / ( sha1 ) )
B <- dgamma( gammasq[4], 2, scale = 1 / ( sha2 ) )
delta4 <- A / ( A + ( 4 * B ) )
if ( U < delta4 ) {
h <- rbeta( 1, 2, 4 )
} else {
h <- rbeta( 1, 1, 5 )
}
if ( i > burnin & ( i - burnin ) %% every == 0 ) {
holdlambda[ , ( i - burnin ) / every ] <- lambdastar
holdgammasq[ , ( i - burnin ) / every ] <- gammasq
holdLAMacsept[ , ( i - burnin ) / every ] <- lambdaaccept
holdalpha[ ( i - burnin ) / every ] <- alpha
holdbeta[ , ( i - burnin ) / every ] <- beta
holdsigmasq[ ( i - burnin ) / every ] <- sigmasq
holdpsi[ , ( i - burnin ) / every ] <- psi
holdW[ ( i - burnin ) / every ] <- w
holdH[ ( i - burnin ) / every ] <- h
}
} ## iterations are over
list( alpha = holdalpha, beta = holdbeta, sigmasq = holdsigmasq, psi = holdpsi,
lambda = holdlambda, gammasq = holdgammasq, H = holdH, W = holdW)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.