Nothing
spl.X<-function(x,df=10){
rk<-function(x,z){
((z-0.5)^2-1/12)*((x-0.5)^2-1/12)/4-
((abs(x-z)-0.5)^4-(abs(x-z)-0.5)^2/2+7/240)/24
}
knotnum<-df+2
knots <- stats::quantile(x, probs = seq(0, 1, length =knotnum ))
knots <- unique(knots[!knots %in% knots[c(1, knotnum)]])
n<-length(x)
X<-outer(x,knots,FUN=rk)
return(X)
}
simulation_var<-function(n=n,cov=cov, E=E,G=G,EG=EG, res=res, model=model, pis=pis){
nZ<-ncol(G)
nX<-ncol(cov)+1
if (model=='binomial') Vinv<-diag(pis*(1-pis))
if (model=='gaussian') Vinv<-diag(rep(1,n))
Xtu <- as.matrix(cbind(rep(1,n), cov, G, EG ))
A <- t(Xtu) %*% Vinv %*% Xtu
D <- matrix(0,nrow=(nZ+nZ+nX),ncol=(nZ+nZ+nX))
for (i in 1:n){
Ui2<-as.vector(Xtu[i,] * res[i])
D <- D+(Ui2)%*%t(Ui2)
}
A<-A/n
A11<-A[1:(nX),1:(nX)]
A12<-A[1:(nX),(nX+1):(nX+nZ+nZ)]
A21<-A[(nX+1):(nX+nZ+nZ),1:(nX)]
A22<-A[(nX+1):(nX+nZ+nZ),(nX+1):(nX+nZ+nZ)]
D<-D/(n-1-nZ-nX)
B11<-D[1:(nX),1:(nX)]
B12<-D[1:(nX),(nX+1):(nX+nZ+nZ)]
B21<-t(B12)
B22<-D[(nX+1):(nX+nZ+nZ),(nX+1):(nX+nZ+nZ)]
A11.edecomp <- eigen(A11,symmetric=T) #general inverse
A11.ev.inv <- ifelse(abs(A11.edecomp$values) >1e-8, 1/A11.edecomp$values, 0)
A11inv <- t(A11.edecomp$vectors %*% (t(A11.edecomp$vectors) * A11.ev.inv))
CovS <- (A22 - A21 %*% A11inv %*% A12 )
CovS.edecomp <- eigen(CovS,symmetric=T)
CovS.ev <- ifelse(abs(CovS.edecomp$values) >1e-8, CovS.edecomp$values, 0)
V11 <- t(CovS.edecomp$vectors %*% (t(CovS.edecomp$vectors) * sqrt(CovS.ev)))
return(V11)
}
calculation<-function(method=method, n.perm=n.perm, res=res, XdE=XdE, nZ=nZ, n=n, V11=V11, pow=pow, Ts_main=Ts_main,Ts_int=Ts_int,Ts=Ts, G=G, EG=EG, pis=pis, cov1=cov1, model=model){
pPerm0_main = rep(NA,length(pow))
pPerm0=rep(NA,length(pow))
T0s_main = matrix(NA,n.perm,length(pow))
P0s_main <-matrix(NA,n.perm,length(pow))
pPerm0_int = rep(NA,length(pow))
T0s_int = matrix(NA,n.perm,length(pow))
P0s_int <-matrix(NA,n.perm,length(pow))
T0s = matrix(NA,n.perm,length(pow))
P0s <-matrix(NA,n.perm,length(pow))
Y0 <- numeric(n)
for (b in 1:n.perm){
if (method=='ResidPerm') {res.p<-sample(res,n,replace=T)
U0_main <- as.vector(t(G)%*% t(t(res.p)))
U0_int <- as.vector(t(EG) %*% t(t(res.p)))
} else if (method=='Bootstrap') {
for(i in 1:n) Y0[i] <- sample(c(1,0), 1, prob=c(pis[i], 1-pis[i]))
tdat4=data.frame(trait=Y0,cov1)
temp<-stats::glm(trait~.,data=tdat4,family=model)
pis.p<-temp$fitted.value
res.p<-tdat4$trait-pis.p
U0_main <- as.vector(t(G)%*% t(t(res.p)))
U0_int <-as.vector(t(EG) %*% t(t(res.p)))
} else if (method=='Simulation'){
nX<-ncol(XdE)+2
Us <- stats::rnorm(nZ+nZ)
U0 <- ( V11 %*% Us) * sqrt(n)
U0_main<-U0[1:nZ]
U0_int<-U0[(nZ+1):(nZ+nZ)]
}
for (j in 1:length(pow)){
if (pow[j] < Inf){
T0s_main[b,j] = round( sum( U0_main^pow[j]), digits = 8)
T0s_int[b,j] = round( sum( U0_int^pow[j]), digits = 8)}
if (pow[j] == Inf) {
T0s_main[b,j] = round( max(abs(U0_main)), digits = 8)
T0s_int[b,j] = round( max(abs(U0_int)), digits = 8)}
} }
T0s=T0s_main+T0s_int #unweighted version
for (j in 1:length(pow)){
pPerm0_main[j] = round((sum(abs(Ts_main[j])<=abs(T0s_main[1:(n.perm-1),j]))+1)/(n.perm), digits = 8) #why from B-1 and add 1-> to make p-value (0,1]
P0s_main[,j] = (n.perm-rank(abs(T0s_main[,j]))+1)/(n.perm)
pPerm0_int[j] = round((sum(abs(Ts_int[j])<=abs(T0s_int[1:(n.perm-1),j]))+1)/(n.perm), digits = 8)
P0s_int[,j] = (n.perm-rank(abs(T0s_int[,j]))+1)/(n.perm)
pPerm0[j] = round((sum(abs(Ts[j])<=abs(T0s[1:(n.perm-1),j]))+1)/(n.perm), digits = 8)
P0s[,j] = (n.perm-rank(abs(T0s[,j]))+1)/(n.perm)
}
minp0_main<-apply(P0s_main,1,min)
minp0_int<-apply(P0s_int,1,min)
minp0<-apply(P0s,1,min)
twomin<-cbind(minp0_main,minp0_int)
min_min<-apply(twomin,1,min)
fisher_min<-minp0_main*minp0_int
Paspu_min<-(sum(min_min<=min(min(pPerm0_main),min(pPerm0_int)))+1)/(n.perm+1)
Paspu_fisher<-(sum(log(fisher_min)<=log(min(pPerm0_main))+log(min(pPerm0_int)))+1)/(n.perm+1) #update log here for fisher apr12
Paspu_unweighted<-(sum(minp0<=min(pPerm0))+1)/(n.perm+1)
Ts <- c(Ts,min(pPerm0),min(min(pPerm0_main),min(pPerm0_int)),min(pPerm0_main)*min(pPerm0_int) )
pvs <- c(pPerm0,Paspu_unweighted, Paspu_min, Paspu_fisher)
names(Ts) <- c(paste("aGEj", pow, sep = ""), "aGEj","aGEj_minp","aGEj_fisher")
names(pvs) = names(Ts)
return(pvs)
}
RobustSim<-function(n=n, cov=cov, G=G, XUs=XUs, res=res, n.perm=n.perm, pow=pow, Vinv=Vinv, model=model, Y=Y, phi=phi){
nZ <-ncol(G)
nX <- ncol(cov) +1
Xtu <- as.matrix(cbind( rep(1,nrow(XUs)), cov, XUs )) #add intercept
A <- t(Xtu) %*% Vinv %*% Xtu
B <- matrix(0,nrow=(nZ+nX),ncol=(nZ+nX))
U <- as.vector(t(XUs)%*%res) / phi
for (i in 1:n){
B <- B+1/(phi^2)*tcrossprod(Xtu[i,])*(res[i]^2) }
B <- B / n
A <- A / n
A11inv <- solve(A[1:nX, 1:nX])
A12 <- A[1:nX, (nX+1):(nZ+nX)]
A21 <- A[(nX+1):(nZ+nX),1:nX]
B22 <- B[(nX+1):(nZ+ nX), (nX+1):(nZ+nX)]
B12 <- B[1:nX, (nX+1):(nZ+nX)]
B21 <- B[(nX+1):(nZ+nX),1:nX]
B11 <- B[1:nX, 1:nX]
A22 <- A[(nX+1):(nZ+ nX), (nX+1):(nZ+nX)]
CovS <- (A22 - A21 %*% A11inv %*% A12 )
CovS.edecomp <- eigen(CovS,symmetric=T)
CovS.ev <- zapsmall(CovS.edecomp$values)
V22 <- t(CovS.edecomp$vectors %*% (t(CovS.edecomp$vectors) * sqrt(CovS.ev))) #V22 is the standard deviation
Ts <- rep(NA, length(pow))
for (j in 1:length(pow)) {
if (pow[j] < Inf) {
Ts[j] <- sum((U)^pow[j])
}
else {
Ts[j] <- max(abs(U))
}
}
var <- sum(CovS)
asyT1 <- (1-stats::pnorm(abs(Ts[1]/sqrt(n)), mean=0, sd=sqrt(var)))*2 #2-sided test
CovS.ev1 <- CovS.edecomp$values[zapsmall(CovS.edecomp$values)>0]
asyT2 <- survey::pchisqsum(Ts[2]/n, rep(1, length(CovS.ev1)), CovS.ev1, lower.tail = F, method = "sad")
pPerm0 <- rep(NA,length(pow))
T0s <- matrix(NA,n.perm,length(pow))
for (b in 1:n.perm){
Us <- stats::rnorm(nZ, 0, 1)
U0 <- (V22 * sqrt(n)) %*% Us
for (j in 1:length(pow)){
if (pow[j] < Inf){ T0s[b,j] = round( sum( (U0)^pow[j]), digits = 8) }
if (pow[j] == Inf) {T0s[b,j] = round( max(abs(U0)), digits = 8) }
} }
T0s <- T0s[stats::complete.cases(T0s),]
n.perm <- nrow(T0s)
P0s <- matrix(NA,n.perm,length(pow))
for (j in 1:length(pow)){
pPerm0[j] <- round((sum(abs(Ts[j])<=abs(T0s[1:(n.perm-1),j]))+1)/(n.perm), digits = 8)
P0s[,j] <- (n.perm-rank(abs(T0s[,j]))+1)/(n.perm)
}
minp0 <- apply(P0s,1,min)
fisherp0 <- apply(P0s,1,function(x) sum(log(x)))
Paspu <- (sum(minp0<=min(pPerm0))+1)/(n.perm+1)
Pfisher <- (sum(fisherp0<=sum(log(pPerm0)))+1)/(n.perm+1)
Ts <- c(Ts, min(pPerm0),sum(log(pPerm0)))
pvs <- c(pPerm0, Paspu,Pfisher)
names(pvs) <- c(paste("aGEsm", pow, sep = ""), "aGEsm","aGEsm_fisher")
return(list( pvs = pvs, nperm=n.perm))
}
GLMMaSPU <- function(tdat1=tdat1,model=model,G=G, cov=cov, nc=nc, n=n){
nZ <- ncol(G)
if (nc == 1) {
fixedmod <- "trait ~ cov"
} else {
fixedmod <- paste("trait ~", paste(names(tdat1)[2:(1 + nc)], collapse = " + "))}
if (nZ == 1) {
randommod <- "~ 0 + G"} else {
randommod <- paste("~ 0 +", paste(names(tdat1)[(nc + 2):(1 + nc + nZ)], collapse = " + "))
}
if (model == 'binomial'){
fit1 <- suppressWarnings(MASS::glmmPQL(fixed = stats::as.formula(fixedmod), data = tdat1, random = list(id = nlme::pdIdent(stats::as.formula(randommod))), family = model, verbose = F ))
pis<-tryCatch(as.numeric(stats::predict(fit1,data=tdat1,type='response',level=1)),error=function(c) NA, warning = function(c) NA)
varc <- tryCatch(nlme::VarCorr(fit1), error=function(c) NA, warning = function(c) NA)
tau <- ifelse(all(is.na(varc)), NA, as.numeric(varc[1,1])) #because G follows the same distribution
phi <- ifelse(all(is.na(varc)), NA, as.numeric(varc[nrow(varc),1]))
W <- diag(pis*(1-pis))/phi
} else if (model=='gaussian'){
fit1 <- suppressWarnings(MASS::glmmPQL(fixed = stats::as.formula(fixedmod), data = tdat1, random = list(id = nlme::pdIdent(stats::as.formula(randommod))), family = model, verbose = F ))
pis <- tryCatch(as.numeric(stats::predict(fit1,data=tdat1,type='response',level=1)),error=function(c) NA, warning = function(c) NA)
varc <- tryCatch(nlme::VarCorr(fit1), error=function(c) NA, warning = function(c) NA)
phi <- as.numeric(varc[nrow(varc),1])
tau <- as.numeric(varc[1,1])
W <- diag(n)/phi
}
return(list(pis=pis, phi=phi, tau=tau, W=W))
}
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.