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

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 )
```