# Returns MW estimator for P(Y1>Y2)
eblup.auc.REML <- function(Y1,Y2){
Y1.df <- as.data.frame(Y1) # a data frame for lmer
Y2.df <- as.data.frame(Y2) # a data frame for lmer
############## EBLUPS #################
# factor to index random subject effects
Y1.df$subject <- factor(1:nrow(Y1))
Y2.df$subject <- factor(1:nrow(Y2))
# Reshape the data matrix for lmer eblup
Y1.long <- reshape2::melt(Y1.df,id.vars="subject")
Y2.long <- reshape2::melt(Y2.df,id.vars="subject")
# fit a linear mixed model
fit1 <- lmer(value~1|subject,data = Y1.long,REML = TRUE)
fit2 <- lmer(value~1|subject,data = Y2.long,REML = TRUE)
# extract parameter estimates
s1 <- summary(fit1)
mu1.hat <- s1$coefficients[1]
eblups1 <- ((ranef(fit1)$subject[1][,1]))
mod1 <- tidy(fit1)
sigmaX.hat <- mod1$estimate[2]
sigmaEps.hat <- mod1$estimate[3]
s2 <- summary(fit2)
mu2.hat <- s2$coefficients[1]
eblups2 <- ((ranef(fit2)$subject[1][,1]))
mod2 <- tidy(fit2)
sigmaY.hat <- mod2$estimate[2]
sigmaEta.hat <- mod2$estimate[3]
estimates1 <- eblups1+mu1.hat
estimates2 <- eblups2+mu2.hat
# compute Mann-Whitney statistic for estimates2,1
A.hat <- mean(outer(estimates1,estimates2,">"))
return(list(A=A.hat,sigmas=c(sigmaX.hat,sigmaY.hat),taus=c(sigmaEps.hat,sigmaEta.hat)
,blups1=estimates1,blups2=estimates2))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.