simex.auc.naive.RM <- function(Y2,Y1,
K=1000,LAM1=seq(0.1,1.9,0.2),LAM2=seq(0.1,1.9,0.2),
tau1.reml=NA,tau2.reml=NA){
# Y1, Y2 are matrices of RM (Y2 is the diseased pop.)
# K is the resampling size for the SIMEX procedure
# LAM1,2 are the vectors of Lambda for simex resampling
# can accept REML estimates of variance components
# otherwise, uses ANOVA estimates
n <- nrow(Y1)
p <- ncol(Y1)
# ANOVA estimates in case REML not supplied
if(is.na(tau1.reml)){
y1_bar <- rowMeans(Y1)
y2_bar <- rowMeans(Y2)
# sample estimates
mu1 <- mean(Y1)
mu2 <- mean(Y2)
SSB1 <- p * sum( (y1_bar-mu1)^2)
SSB2 <- p * sum( (y2_bar-mu2)^2)
tau1_hat <-sqrt(( 1 / ((p-1)*n)) * sum((Y1-do.call(cbind,replicate(p,y1_bar,simplify = FALSE)))^2))
tau2_hat <-sqrt(( 1 / ((p-1)*n)) * sum((Y2-do.call(cbind,replicate(p,y2_bar,simplify = FALSE)))^2))
}
# in case REML estimates are supplied
else{
tau1_hat <- tau1.reml
tau2_hat <- tau2.reml
}
TC <- matrix(data = NA,nrow = 0,ncol = 3)
for(lam1 in LAM1){
for(lam2 in LAM2){
w1 <- matrix(rnorm(K*p*n,0,sqrt(lam1)*tau1_hat), n, p*K )
w2 <- matrix(rnorm(K*p*n,0,sqrt(lam2)*tau2_hat), n, p*K )
Y1.rep <- do.call(cbind,replicate(K,Y1,simplify = FALSE))
W1 <- Y1.rep + w1
Y2.rep <- do.call(cbind,replicate(K,Y2,simplify = FALSE))
W2 <- Y2.rep + w2
TC <- rbind(TC,c(lam1,lam2, mw.outer.matrix(W1,W2)))
}
}
TC <- as.data.frame(TC)
names(TC) <- c("L1","L2","E_theta_N")
# Now we obtain the intercept from the linear regression model
TC$L1 <- 1+TC$L1
TC$L2 <- 1+TC$L2
TC$L1sq <- TC$L1^2
TC$L2sq <- TC$L2^2
fit1 <- lm(formula = E_theta_N ~ L1+L2,data = TC)
fit2 <- lm(formula = E_theta_N ~ L1+L2+L1sq+L2sq,data = TC)
s1 <- as.numeric(fit1$coefficients[1])
s2 <- as.numeric(fit2$coefficients[1])
# trimming of values >1
s1 <- min(c(s1,1))
s2 <- min(c(s2,1))
v <- c(s1,s2)
names(v) <- c("Linear","Quadratic")
return(v)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.