MATH 6643 Summer 2012 Applications of Mixed Models/Students/smithce/SomethingUseful

From Wiki1

< MATH 6643 Summer 2012 Applications of Mixed Models | Students | smithce
Revision as of 23:42, 19 August 2012 by Smithce (Talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
plotSimulateLME <-
function(x, form = y ~ x | df * method, df = attr(x, "df"), weights,
xlab = "Empirical p-value",
ylab = "Nominal p-value", 
pobs = NULL,
alpha = 0.05, dalpha = 0.005, digits = 4,
xlim = NULL,
ylim = NULL, ...)
{
ML <- !is.null(x$null$ML)
if(ML) {
if (is.null(x$alt$ML))
stop("Plot method only implemented for comparing models")
okML <- x$null$ML[, "info"] < 8 & x$alt$ML[, "info"] < 8
}
REML <- !is.null(x$null$REML)
if(REML) {
if (is.null(x$alt$REML))
stop("Plot method only implemented for comparing models")
okREML <- x$null$REML[, "info"] < 8 & x$alt$REML[, "info"] < 8
}
if (is.null(df)) {
stop("No degrees of freedom specified")
}
if ((ldf <- length(df)) > 1) {
df <- sort(unique(df))
if (missing(weights)) {
weights <- rep.int(1/ldf, ldf)
} else {
if (!identical(weights,FALSE) && length(weights) != ldf)
stop("Degrees of freedom and weights must have the same length")
}
} else {
weights <- FALSE
}
useWgts <- (length(weights) != 1)
if (any(df < 0)) {
stop("Negative degrees of freedom not allowed")
} else {
if ((ldf == 1) && (df == 0)) {
stop("More than one degree of freedom is needed when one them is zero.")
}
}
if (ML) {
MLstat <-
rev(sort(2 * pmax(0, x$alt$ML[okML, "logLik"] - x$null$ML[okML,"logLik"])))
MLy <- lapply(df,
function(df, x) {
if (df > 0) 1 - pchisq(x, df) else 1*(x == 0)
}, x = MLstat)
dfC <- paste("df",df,sep="=")
if (useWgts) {                      # has weights
if (ldf == 2) {                   # will interpolate
MLy <-
c(MLy1, weights[1] * MLy1 + weights[2] * MLy2, MLy2)
MLdf <- rep(c(dfC[1], paste("Mix(",df[1],",",df[2],")",sep=""),
dfC[2]), rep(length(MLstat), ldf + 1))
} else {
aux <- weights[1] * MLy1
auxNam <- paste("Mix(",df[1],sep="")
for(i in 2:ldf) {
aux <- aux + weights[i] * MLyi
auxNam <- paste(auxNam, ",", df[i],sep="")
}
auxNam <- paste(auxNam, ")",sep="")
MLy <- c(unlist(MLy), aux)
MLdf <- rep(c(dfC, auxNam), rep(length(MLstat), ldf + 1))
}
MLx <- rep((1:length(MLstat) - 0.5)/length(MLstat), ldf + 1)
} else {
MLy <- unlist(MLy)
MLdf <- rep(dfC, rep(length(MLstat), ldf))
MLx <- rep((1:length(MLstat) - 0.5)/length(MLstat), ldf)
}
auxInd <- MLdf != "df=0"
meth <- rep("ML", length(MLy))
Mdf <- MLdf
} else {
MLy <- MLdf <- MLx <- auxInd <- meth <- Mdf <- NULL
}
if (REML) {
REMLstat <- rev(sort(2 * pmax(0, x$alt$REML[okREML, "logLik"] -
x$null$REML[okREML, "logLik"])))
REMLy <- lapply(df,
function(df, x) {
if (df > 0) {
1 - pchisq(x, df)
} else {
val <- rep(0, length(x))
val[x == 0] <- 1
val
}
}, x = REMLstat)
dfC <- paste("df",df,sep="=")
if (useWgts) {                      # has weights
if (ldf == 2) {                   # will interpolate
REMLy <-
c(REMLy1, weights[1] * REMLy1 + weights[2] * REMLy2, REMLy2)
REMLdf <- rep(c(dfC[1], paste("Mix(",df[1],",",df[2],")",sep=""),
dfC[2]), rep(length(REMLstat), ldf + 1))
} else {
aux <- weights[1] * REMLy1
auxNam <- paste("Mix(",df[1],sep="")
for(i in 2:ldf) {
aux <- aux + weights[i] * REMLyi
auxNam <- paste(auxNam, ",", df[i],sep="")
}
auxNam <- paste(auxNam, ")",sep="")
REMLy <- c(unlist(REMLy), aux)
REMLdf <- rep(c(dfC, auxNam), rep(length(REMLstat), ldf + 1))
}
REMLx <- rep((1:length(REMLstat) - 0.5)/length(REMLstat), ldf + 1)
} else {
REMLy <- unlist(REMLy)
REMLdf <- rep(dfC, rep(length(REMLstat), ldf))
REMLx <- rep((1:length(REMLstat) - 0.5)/length(REMLstat), ldf)
}
auxInd <- c(auxInd, REMLdf != "df=0")
meth <- c(meth, rep("REML", length(REMLy)))
Mdf <- c(Mdf, REMLdf)
} else {
REMLy <- REMLdf <- REMLx <- NULL
}
meth <- meth[auxInd]
Mdf <- Mdf[auxInd]
Mdf <- ordered(Mdf, levels = unique(Mdf))
frm <- data.frame(x = c(MLx, REMLx)[auxInd], y = c(MLy, REMLy)[auxInd],
df = Mdf,  method = meth) 
MLresult <- frm[frm$method=="ML",]
REMLresult <- frm[frm$method=="REML",]
if( nrow(MLresult)>0 ){
message("Maximum Likelihood \n")
pML.Emp.to.Nom <- MLresult[MLresult$x>(alpha-dalpha)&MLresult$x<(alpha+dalpha),1:2]
colnames(pML.Emp.to.Nom) <- c("Empirical p-value", "Nominal p-value") 
if (nrow(pML.Emp.to.Nom) == 0) {
message ( "p-value not obtained at desired resolution.  Increase 'dalpha' or re-run 
simulation with more iterations" )
} else print(pML.Emp.to.Nom[,2:1], row.names=F, digits=digits)
if(is.numeric(pobs)){
message(paste("Nominal p-value", pobs, 
"approximately equal to empirical p-value", MLresult$x[which.min(abs(MLresult$y - pobs))]))
}
message("")
pML.Nom.to.Emp <- MLresult[MLresult$y>(alpha-dalpha)&MLresult$y<(alpha+dalpha),1:2]
colnames(pML.Nom.to.Emp) <- c("Empirical p-value", "Nominal p-value") 
}
if( nrow(REMLresult)>0 ){
message("REML \n")
pREML.Emp.to.Nom <- REMLresult[REMLresult$x>(alpha-dalpha)&REMLresult$x<(alpha+dalpha),1:2]
colnames(pREML.Emp.to.Nom) <- c("Empirical p-value", "Nominal p-value") 
if (nrow(pREML.Emp.to.Nom) == 0) {
message ( "p-value not obtained at desired resolution.  Increase 'dalpha' or re-run 
simulation with more iterations" )
} else print(pREML.Emp.to.Nom[,2:1], row.names=F, digits=digits)
if(is.numeric(pobs)){
message(paste("Nominal p-value", pobs, 
"approximately equal to empirical p-value", REMLresult$x[which.min(abs(REMLresult$y - pobs))]))
}
}
if(nrow(MLresult)>0 & nrow(REMLresult)>0) {
if(is.null(xlim)) {
xlim <- c(0,
max(c(MLresult$x[which.min(abs(MLresult$x - alpha))],
MLresult$x[which.min(abs(MLresult$y - alpha))],
MLresult$x[which.min(abs(MLresult$y - pobs))],
REMLresult$x[which.min(abs(REMLresult$x - alpha))],
REMLresult$x[which.min(abs(REMLresult$y - alpha))],
REMLresult$x[which.min(abs(REMLresult$y - pobs))])) + dalpha) }
if(is.null(ylim)) {
ylim <- c(0,
max(c(MLresult$y[which.min(abs(MLresult$x - alpha))],
MLresult$y[which.min(abs(MLresult$y - alpha))],
MLresult$y[which.min(abs(MLresult$y - pobs))],
REMLresult$y[which.min(abs(REMLresult$x - alpha))],
REMLresult$y[which.min(abs(REMLresult$y - alpha))],
REMLresult$y[which.min(abs(REMLresult$y - pobs))])) + dalpha)}
op <- par(mfrow = c(1,2))
plot(MLresult$x, MLresult$y, type='l',
main = "Maximum Likelihood",
ylab = "Nominal p-value",
xlab = "Empirical p-value",
xlim = xlim,
ylim = ylim)
# Empirical to Nominal cut-off alpha
lines(x=c(xlim[1],MLresult$x[which.min(abs(MLresult$x - alpha))],MLresult$x[which.min(abs(MLresult$x - alpha))]),y=c(MLresult$y[which.min(abs(MLresult$x   - alpha))],MLresult$y[which.min(abs(MLresult$x - alpha))],0),lty="dashed",col="red")
# Nominal to Empirical cut-off alpha (indicative of degree to which conservative)
lines(x=c(xlim[1],MLresult$x[which.min(abs(MLresult$y - alpha))],MLresult$x[which.min(abs(MLresult$y - alpha))]),y=c(MLresult$y[which.min(abs(MLresult$y - alpha))],MLresult$y[which.min (abs(MLresult$y - alpha))],0),lty="dashed",col="pink")
# Observed nominal to empirical (closest approximation)
if(!is.null(pobs)) {  
lines(x=c(xlim[1],MLresult$x[which.min(abs(MLresult$y - pobs))],MLresult$x[which.min(abs(MLresult$y - pobs))]),y=c(MLresult$y[which.min(abs(MLresult$y - pobs))],MLresult$y[which.min(abs (MLresult$y - pobs))],0),lty="dashed",col="blue")
}
 
plot(REMLresult$x, REMLresult$y, type='l',
main = "REML",
ylab = "Nominal p-value",
xlab = "Empirical p-value",
xlim = xlim,
ylim = ylim)
# Empirical to Nominal cut-off alpha
lines(x=c(xlim[1],REMLresult$x[which.min(abs(REMLresult$x - alpha))],REMLresult$x[which.min(abs(REMLresult$x - alpha))]),y=c(REMLresult$y[which.min(abs(REMLresult$x - alpha))],REMLresult$y[which.min(abs(REMLresult$x - alpha))],0),lty="dashed",col="red")
# Nominal to Empirical cut-off alpha (indicative of degree to which conservative)
lines(x=c(xlim[1],REMLresult$x[which.min(abs(REMLresult$y - alpha))],REMLresult$x[which.min(abs(REMLresult$y - alpha))]),y=c(REMLresult$y[which.min(abs(REMLresult$y - alpha))],REMLresult$y[which.min(abs(REMLresult$y - alpha))],0),lty="dashed",col="pink")
# Observed nominal to empirical (closest approximation)
if(!is.null(pobs)) { 
lines(x=c(xlim[1],REMLresult$x[which.min(abs(REMLresult$y - pobs))],REMLresult$x[which.min(abs(REMLresult$y - pobs))]),
y=c(REMLresult$y[which.min(abs(REMLresult$y - pobs))],REMLresult$y[which.min(abs(REMLresult$y - pobs))],0),
lty="dashed",col="blue")}
par(op)
}else if(nrow(MLresult)>0) {
if(is.null(xlim)) {
xlim <- c(0,
max(c(MLresult$x[which.min(abs(MLresult$x - alpha))],
MLresult$x[which.min(abs(MLresult$y - alpha))],
MLresult$x[which.min(abs(MLresult$y - pobs))])) + dalpha) }
if(is.null(ylim)) {
ylim <- c(0,
max(c(MLresult$y[which.min(abs(MLresult$x - alpha))],
MLresult$y[which.min(abs(MLresult$y - alpha))],
MLresult$y[which.min(abs(MLresult$y - pobs))])) + dalpha) }
plot(MLresult$x, MLresult$y, type='l',
main = "Maximum Likelihood",
ylab = "Nominal p-value",
xlab = "Empirical p-value",
xlim = xlim,
ylim = ylim)
# Empirical to Nominal cut-off alpha
lines(x=c(xlim[1],MLresult$x[which.min(abs(MLresult$x - alpha))],MLresult$x[which.min(abs(MLresult$x - alpha))]),y=c(MLresult$y[which.min(abs(MLresult$x - alpha))],MLresult$y[which.min (abs(MLresult$x - alpha))],0),lty="dashed",col="red")
# Nominal to Empirical cut-off alpha (indicative of degree to which conservative)
lines(x=c(xlim[1],MLresult$x[which.min(abs(MLresult$y - alpha))],MLresult$x[which.min(abs(MLresult$y - alpha))]),y=c(MLresult$y[which.min(abs(MLresult$y - alpha))],MLresult$y[which.min (abs(MLresult$y - alpha))],0),lty="dashed",col="pink")
# Observed nominal to empirical (closest approximation)
if(!is.null(pobs)) {
lines(x=c(xlim[1],MLresult$x[which.min(abs(MLresult$y - pobs))],MLresult$x[which.min(abs(MLresult$y - pobs))]),y=c(MLresult$y[which.min(abs(MLresult$y - pobs))],MLresult$y[which.min(abs (MLresult$y - pobs))],0),lty="dashed",col="blue")}
}else if(nrow(REMLresult)>0){
if(is.null(xlim)) {
xlim <- c(0,
max(c(REMLresult$x[which.min(abs(REMLresult$x - alpha))],
REMLresult$x[which.min(abs(REMLresult$y - alpha))],
REMLresult$x[which.min(abs(REMLresult$y - pobs))])) + dalpha) }
if(is.null(ylim)) {
ylim <- c(0,
max(c(REMLresult$y[which.min(abs(REMLresult$x - alpha))],
REMLresult$y[which.min(abs(REMLresult$y - alpha))],
REMLresult$y[which.min(abs(REMLresult$y - pobs))])) + dalpha)}
plot(REMLresult$x, REMLresult$y, type='l',
main = "REML",
ylab = "Nominal p-value",
xlab = "Empirical p-value",
xlim = xlim,
ylim = ylim)
# Empirical to Nominal cut-off alpha
lines(x=c(xlim[1],REMLresult$x[which.min(abs(REMLresult$x - alpha))],REMLresult$x[which.min(abs(REMLresult$x - alpha))]),y=c(REMLresult$y[which.min(abs(REMLresult$x - alpha))],REMLresult$y[which.min(abs(REMLresult$x - alpha))],0),lty="dashed",col="red")
# Nominal to Empirical cut-off alpha (indicative of degree to which conservative)
lines(x=c(xlim[1],REMLresult$x[which.min(abs(REMLresult$y - alpha))],REMLresult$x[which.min(abs(REMLresult$y - alpha))]),y=c(REMLresult$y[which.min(abs(REMLresult$y - alpha))],REMLresult$y[which.min(abs(REMLresult$y - alpha))],0),lty="dashed",col="pink")
# Observed nominal to empirical (closest approximation)
if(!is.null(pobs)) {
lines(x=c(xlim[1],REMLresult$x[which.min(abs(REMLresult$y - pobs))],REMLresult$x[which.min(abs(REMLresult$y - pobs))]),y=c(REMLresult$y[which.min(abs(REMLresult$y - pobs))],REMLresult$y [which.min(abs(REMLresult$y - pobs))],0),lty="dashed",col="blue")}
}
}
## Example ##
library(spidadev)
library(nlme)
dd <- hs1
str(dd)
fit1 <- lme( mathach ~ ses, data=dd,
random = ~ 1|school, 
na.action = na.exclude )
fit2 <- lme( mathach ~ ses, data=dd,
random = ~ 1 + ses|school, # Model 2 has more complex random part
na.action = na.exclude )
anova( fit1, fit2 )  # Test significance of improvement of random slope
# However this test is conservative
# Simulate to adjust p-values
system.time( sim.out <- simulate( fit1, m2 = fit2, nsim = 1000) )
plotSimulateLME( sim.out , dalpha = .001 )
plotSimulateLME( sim.out , dalpha = .005, pobs = 0.1021 )
system.time( sim.out <- simulate( fit1, m2 = fit2, nsim = 1000, method = "ML") )
plotSimulateLME( sim.out , dalpha = .0025, pobs = 0.1021 )
system.time( sim.out <- simulate( fit1, m2 = fit2, nsim = 1000, method = "REML") )
plotSimulateLME( sim.out , dalpha = .0025, pobs = 0.1021 )
Personal tools