R/PlotGLMERContinuous.R
In timnewbold/StatisticalModels: What the package does (short line)
Defines functions
PlotGLMERContinuous
Documented in
PlotGLMERContinuous
PlotGLMERContinuous<-function(model,data,effects,byContEffect=NULL,
otherContEffects=character(0),
otherFactors=character(0),xlab,ylab,
byFactor=NULL,zlab=NULL,outDir=NULL,
logLink="n",plotRug=FALSE,seMultiplier=1.96,
params=list(),add=FALSE,xlim=NULL,ylim=NULL,
zlim=NULL,line.cols=NULL,line.types=NULL,
plotUncertainty=TRUE,nPanels=1,main=NULL,
yDiv=1,transformX=FALSE){
if ((length(effects)>1) && (!is.null(byFactor))){
stop("If plotting by a factor is specified, only one effect can be plotted")
}
if ((!is.null(byFactor)) && (!is.null(byContEffect))){
stop("You cannot specify both a factor and a continuous effect to divide the results")
}
if ((length(effects)==1) && (is.null(byFactor)) && (is.null(byContEffect))){
eval(substitute(newdat<-data.frame(seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length.out=100)),list(x=effects)))
names(newdat)<-effects
for (o in otherContEffects){
newdat[,o] <- median(data[,o])
}
# Add a column to newdat for each of otherFactors
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
newdat[,names(model@frame)[1]] <- 0
mm<-model.matrix(terms(model),newdat)
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
newdat$y<-mm %*% fixef(model)
if (logLink=="e"){
newdat$yplus<-exp(newdat$y+seMultiplier*sqrt(pvar1))
newdat$yminus<-exp(newdat$y-seMultiplier*sqrt(pvar1))
newdat$y<-exp(newdat$y)
} else if (logLink=="10"){
newdat$yplus<-10^(newdat$y+seMultiplier*sqrt(pvar1))
newdat$yminus<-10^(newdat$y-seMultiplier*sqrt(pvar1))
newdat$y<-10^(newdat$y)
} else if (logLink=="n") {
newdat$yplus<-(newdat$y+seMultiplier*sqrt(pvar1))
newdat$yminus<-(newdat$y-seMultiplier*sqrt(pvar1))
newdat$y<-(newdat$y)
} else if (logLink=="b"){
newdat$yplus<-1/(1+exp(-(newdat$y+seMultiplier*sqrt(pvar1))))
newdat$yminus<-1/(1+exp(-(newdat$y-seMultiplier*sqrt(pvar1))))
newdat$y<-1/(1+exp(-(newdat$y)))
} else {
stop("Error: the specified log link is not supported")
}
if (!is.null(outDir)){
png(paste(outDir,"/",effects," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mgp=c(2.5,1,0))
par(mar=c(4,4,1,1))
par(las=1)
par(cex.lab=1.5)
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
if (!transformX){
xlims <- c(min(newdat[,effects]),max(newdat[,effects]))
} else {
xlims <- c(exp(min(newdat[,effects])),exp(max(newdat[,effects])))
}
if (!is.null(xlim)){
xlims <- xlim
}
if (!add){
if (is.null(ylim)){
if(!transformX){
plot(-1e5,-1e5,xlim=xlims,
ylim=c(min(newdat$yminus/yDiv),max(newdat$yplus/yDiv)),
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=c(min(newdat$yminus/yDiv),max(newdat$yplus/yDiv)),
xlab=xlab,ylab=ylab,
log="x",main=main)
}
} else {
if (!transformX){
plot(-1e5,-1e5,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main,
log="x")
}
}
}
line.col<-ifelse(is.null(line.cols),"#000000",line.cols[1])
if(plotUncertainty){
if (transformX){
X.Vec <- exp(c(newdat[,effects], max(newdat[,effects]),
rev(newdat[,effects]), min(newdat[,effects])))
} else {
X.Vec <- c(newdat[,effects], max(newdat[,effects]),
rev(newdat[,effects]), min(newdat[,effects]))
}
Y.Vec <- c(newdat$yminus/yDiv, tail(newdat$yplus/yDiv, 1), rev(newdat$yplus/yDiv), (newdat$yminus/yDiv)[1])
polygon(X.Vec, Y.Vec, col = paste(line.col,"33",sep=""), border = NA)
}
if(transformX){
points(exp(newdat[,effects]),
newdat$y/yDiv,type="l",col=line.col)
} else {
points(newdat[,effects],
newdat$y/yDiv,type="l",col=line.col)
}
if(plotRug){
if(transformX){
rug(eval(substitute(exp(data$x),list(x=effects[1]))))
} else {
rug(eval(substitute(data$x,list(x=effects[1]))))
}
}
if (!is.null(outDir)){
dev.off()
}
} else if ((length(effects)==1) && (!is.null(byFactor))){
if((!is.null(line.cols)) & (length(line.cols)!=length(levels(model@frame[,byFactor])))){
stop("Specified colours must be of the same number as the number of factor levels")
}
if((!is.null(line.types)) & (length(line.types)!=length(levels(model@frame[,byFactor])))){
stop("Specified line types must be of the same number as the number of factor levels")
}
eval(substitute(newdat<-data.frame(seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length.out=100)),list(x=effects)))
names(newdat)<-effects
# Add a column to newdat for each of otherFactors
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
newdat[,names(model@frame)[1]] <- 0
y<-data.frame(x=newdat[,1])
yplus<-data.frame(x=newdat[,1])
yminus<-data.frame(x=newdat[,1])
for (l in levels(model@frame[,byFactor])){
newdat[,byFactor]<-factor(l,levels=levels(model@frame[,byFactor]))
for (o in otherContEffects){
newdat[,o] <- median(data[,o][(data[,byFactor]==l)])
}
mm<-model.matrix(terms(model),newdat)
if(ncol(mm)>length(fixef(model))){
mm <- mm[,-which(!(names(mm[1,]) %in% names(fixef(model))))]
}
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
y[,l]<-mm %*% fixef(model)
if (logLink=="e"){
yplus[,l]<-exp(y[,l]+seMultiplier*sqrt(pvar1))
yminus[,l]<-exp(y[,l]-seMultiplier*sqrt(pvar1))
y[,l]<-exp(y[,l])
} else if (logLink=="10"){
yplus[,l]<-10^(y[,l]+seMultiplier*sqrt(pvar1))
yminus[,l]<-10^(y[,l]-seMultiplier*sqrt(pvar1))
y[,l]<-10^(y[,l])
} else if (logLink=="b"){
yplus[,l]<-1/(1+exp(-(y[,l]+seMultiplier*sqrt(pvar1))))
yminus[,l]<-1/(1+exp(-(y[,l]-seMultiplier*sqrt(pvar1))))
y[,l]<-1/(1+exp(-(y[,l])))
} else if (logLink=="n") {
yplus[,l]<-(y[,l]+seMultiplier*sqrt(pvar1))
yminus[,l]<-(y[,l]-seMultiplier*sqrt(pvar1))
} else {
stop("Error: the specified log link is not supported")
}
try(y[y$x>quantile(data[data[,byFactor]==l,][,effects[1]],0.975),][,l]<-NA,silent=TRUE)
try(yplus[yplus$x>quantile(data[data[,byFactor]==l,][,effects[1]],0.975),][,l]<-NA,silent=TRUE)
try(yminus[yminus$x>quantile(data[data[,byFactor]==l,][,effects[1]],0.975),][,l]<-NA,silent=TRUE)
try(y[y$x<quantile(data[data[,byFactor]==l,][,effects[1]],0.025),][,l]<-NA,silent=TRUE)
try(yplus[yplus$x<quantile(data[data[,byFactor]==l,][,effects[1]],0.025),][,l]<-NA,silent=TRUE)
try(yminus[yminus$x<quantile(data[data[,byFactor]==l,][,effects[1]],0.025),][,l]<-NA,silent=TRUE)
}
if (!is.null(outDir)){
png(paste(outDir,"/",effects," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mgp=c(2.5,1,0))
par(mar=c(4,4,1,1))
par(las=1)
par(cex.lab=1.5)
if (nPanels>1) par(mfrow=c(1,nPanels))
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
allLevels<-1:length(levels(model@frame[,byFactor]))
nEffectsPerPlot<-ceiling(length(allLevels)/nPanels)
effectsPerPlot<-split(allLevels,ceiling(allLevels/nEffectsPerPlot))
if(!is.null(line.cols)){
cols<-split(line.cols,ceiling(allLevels/nEffectsPerPlot))
}
if(!is.null(line.types)){
line.types<-split(line.types,ceiling(allLevels/nEffectsPerPlot))
}
for (p in 1:nPanels){
if (!add){
options(scipen=0)
if(!transformX){
xlims<-c(min(y$x),max(y$x))
} else {
xlims<-c(exp(min(y$x)),exp(max(y$x)))
}
if(!is.null(xlim)){
xlims<-xlim
}
if (is.null(ylim)){
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,log="x",main=main)
}
} else {
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,log="x",main=main)
}
}
}
if (is.null(line.cols)){
cols<-list()
cols[1][[1]]<-rep("#000000",length(levels(model@frame[,byFactor])))
} else {
}
o<-c(o,setdiff(1:length(levels(model@frame[,byFactor])),o))
if (plotUncertainty){
i<-1
for (l in levels(model@frame[,byFactor])[effectsPerPlot[p][[1]]]){
y.temp<-y[,c('x',l)]
y.temp<-na.omit(y.temp)
yplus.temp<-yplus[,c('x',l)]
yplus.temp<-na.omit(yplus.temp)
yminus.temp<-yminus[,c('x',l)]
yminus.temp<-na.omit(yminus.temp)
if(transformX){
X.Vec <- exp(c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x)))
} else {
X.Vec <- c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x))
}
Y.Vec <- c((yminus.temp/yDiv)[,l], tail((yplus.temp/yDiv)[,l], 1), rev((yplus.temp/yDiv)[,l]), (yminus.temp/yDiv)[,l][1])
polygon(X.Vec, Y.Vec, col = paste(cols[p][[1]][i],"33",sep=""), border = NA)
i<-i+1
}
}
i<-1
for (l in levels(model@frame[,byFactor])[effectsPerPlot[p][[1]]]){
if(transformX){
points(exp(y$x),(y/yDiv)[,l],type="l",col=cols[p][[1]][i],lty=ifelse(is.null(line.types),1,line.types[p][[1]][i]))
if (plotRug){
rug(exp(data[data[,byFactor]==l,][,effects[1]]),col=cols[p][[1]][i])
}
} else {
points(y$x,(y/yDiv)[,l],type="l",col=cols[p][[1]][i],lty=ifelse(is.null(line.types),1,line.types[p][[1]][i]))
if (plotRug){
rug(data[data[,byFactor]==l,][,effects[1]],col=cols[p][[1]][i])
}
}
i<-i+1
}
}
} else if ((length(effects)==1) && (!is.null(byContEffect))){
if((!is.null(line.cols)) & (length(line.cols)!=3)){
stop("Specified colours must be of the same number as the number of factor levels")
}
if((!is.null(line.types)) & (length(line.types)!=3)){
stop("Specified line types must be of the same number as the number of factor levels")
}
eval(substitute(newdat<-data.frame(seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length.out=100)),
list(x=effects)))
names(newdat)<-effects
# Add a column to newdat for each of otherFactors
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
if (class(model)=="glmmadmb"){
newdat[,names(model$frame)[1]] <- 0
} else {
newdat[,names(model@frame)[1]] <- 0
}
y<-data.frame(x=newdat[,1])
yplus<-data.frame(x=newdat[,1])
yminus<-data.frame(x=newdat[,1])
byContEffectVals <- c(quantile(data[,byContEffect],probs = 0.025,na.rm=TRUE),
median(data[,byContEffect],na.rm=TRUE),
quantile(data[,byContEffect],probs = 0.975,na.rm=TRUE))
binWidth <- 0.1
prbs <- c(binWidth,0.5-(0.5*binWidth),0.5+(0.5*binWidth),1-binWidth)
brks <- c(-Inf,quantile(data[,byContEffect],
probs = prbs,na.rm = TRUE),Inf)
brks <- brks + seq_along(brks) * 10 * .Machine$double.eps
data$ByContEffectTerciles <- cut(data[,byContEffect],
breaks = brks,
labels = c(1,"NA1",2,"NA2",3))
emptyBin <- any(unlist(lapply(as.list(c(1,2,3)),function(x){
return(length(which(data$ByContEffectTerciles==x))==0)
})))
while(emptyBin){
binWidth <- binWidth + 0.01
prbs <- c(binWidth,0.5-(0.5*binWidth),0.5+(0.5*binWidth),1-binWidth)
brks <- c(-Inf,quantile(data[,byContEffect],
probs = prbs,na.rm = TRUE),Inf) +
seq_along(brks) * 10 * .Machine$double.eps
data$ByContEffectTerciles <- cut(data[,byContEffect],
breaks = brks,
labels = c(1,"NA1",2,"NA2",3))
emptyBin <- any(unlist(lapply(as.list(c(1,2,3)),function(x){
return(length(which(data$ByContEffectTerciles==x))==0)
})))
}
for (i in 1:length(byContEffectVals)){
newdat[,byContEffect] <- byContEffectVals[i]
for (o in otherContEffects){
newdat[,o] <- median(data[,o][(data$ByContEffectTerciles==i)],na.rm=TRUE)
}
mm<-model.matrix(terms(model),newdat)
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
y[,paste(i)]<-mm %*% fixef(model)
if (logLink=="e"){
yplus[,paste(i)]<-exp(y[,paste(i)]+seMultiplier*sqrt(pvar1))
yminus[,paste(i)]<-exp(y[,paste(i)]-seMultiplier*sqrt(pvar1))
y[,paste(i)]<-exp(y[,paste(i)])
} else if (logLink=="10"){
yplus[,paste(i)]<-10^(y[,paste(i)]+seMultiplier*sqrt(pvar1))
yminus[,paste(i)]<-10^(y[,paste(i)]-seMultiplier*sqrt(pvar1))
y[,paste(i)]<-10^(y[,paste(i)])
} else if (logLink=="b"){
yplus[,paste(i)]<-1/(1+exp(-(y[,paste(i)]+seMultiplier*sqrt(pvar1))))
yminus[,paste(i)]<-1/(1+exp(-(y[,paste(i)]-seMultiplier*sqrt(pvar1))))
y[,paste(i)]<-1/(1+exp(-(y[,paste(i)])))
} else if (logLink=="n") {
yplus[,paste(i)]<-(y[,paste(i)]+seMultiplier*sqrt(pvar1))
yminus[,paste(i)]<-(y[,paste(i)]-seMultiplier*sqrt(pvar1))
} else {
stop("Error: the specified log link is not supported")
}
try(y[y$x>quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.975),][,paste(i)]<-NA,silent=TRUE)
try(yplus[yplus$x>quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.975),][,paste(i)]<-NA,silent=TRUE)
try(yminus[yminus$x>quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.975),][,paste(i)]<-NA,silent=TRUE)
try(y[y$x<quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.025),][,paste(i)]<-NA,silent=TRUE)
try(yplus[yplus$x<quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.025),][,paste(i)]<-NA,silent=TRUE)
try(yminus[yminus$x<quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.025),][,paste(i)]<-NA,silent=TRUE)
}
if (!is.null(outDir)){
png(paste(outDir,"/",effects," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mgp=c(2.5,1,0))
par(mar=c(4,4,1,1))
par(las=1)
par(cex.lab=1.5)
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
if (!add){
options(scipen=0)
if(!transformX){
xlims<-c(min(y$x),max(y$x))
} else {
xlims<-c(exp(min(y$x)),exp(max(y$x)))
}
if(!is.null(xlim)){
xlims<-xlim
}
if (is.null(ylim)){
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,log="x",main=main)
}
} else {
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,log="x",main=main)
}
}
}
if (is.null(line.cols)){
cols<-rep("#000000",length(byContEffectVals))
} else {
cols <- line.cols
}
if (plotUncertainty){
for (i in 1:length(byContEffectVals)){
y.temp<-y[,c('x',paste(i))]
y.temp<-na.omit(y.temp)
yplus.temp<-yplus[,c('x',paste(i))]
yplus.temp<-na.omit(yplus.temp)
yminus.temp<-yminus[,c('x',paste(i))]
yminus.temp<-na.omit(yminus.temp)
if(transformX){
X.Vec <- exp(c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x)))
} else {
X.Vec <- c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x))
}
Y.Vec <- c((yminus.temp/yDiv)[,paste(i)], tail((yplus.temp/yDiv)[,paste(i)], 1), rev((yplus.temp/yDiv)[,paste(i)]), (yminus.temp/yDiv)[,paste(i)][1])
polygon(X.Vec, Y.Vec, col = paste(cols[i],"33",sep=""), border = NA)
i<-i+1
}
}
i<-1
for (i in 1:length(byContEffectVals)){
if(transformX){
points(exp(y$x),(y/yDiv)[,paste(i)],type="l",col=cols[i],lty=ifelse(is.null(line.types),1,line.types[i]))
# if (plotRug){
# rug(exp(data[data[,byFactor]==l,][,effects[1]]),col=cols[p][[1]][i])
# }
} else {
points(y$x,(y/yDiv)[,paste(i)],type="l",col=cols[i],lty=ifelse(is.null(line.types),1,line.types[i]))
# if (plotRug){
# rug(data[data[,byFactor]==l,][,effects[1]],col=cols[p][[1]][i])
# }
}
}
} else if (length(effects==2)){
if (!is.null(outDir)){
png(paste(outDir,"/",effects[1]," and ",effects[2]," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mar=c(2,2,2,2))
par(mgp=c(2.5,1,0))
par(las=0)
par(cex.lab=1.5)
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
eval(substitute(e1<-seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length=30),list(x=effects[1])))
eval(substitute(e2<-seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length=30),list(x=effects[2])))
z<-array(dim=c(30,30))
zu<-array(dim=c(30,30))
zl<-array(dim=c(30,30))
gx<-rep(e1,30)
gy<-rep(e2,rep(30,30))
eval(substitute(tf<-matrix(exclude.too.far(gx,gy,data$x,data$y,0.2),
byrow=F,ncol=30),list(x=effects[1],y=effects[2])))
for (x in 1:length(e1)){
for (y in 1:length(e2)){
v1<-e1[x]
v2<-e2[y]
newdat<-data.frame(v1,v2)
names(newdat)<-effects
# Add a column to newdat for each other continuous effect
for (o in otherContEffects){
newdat[,o] <- median(data[,o])
}
# Add a column to newdat for each other factor
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
newdat[,names(model@frame)[1]] <- 0
mm<-model.matrix(terms(model),newdat)
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
if(tf[x,y]==FALSE){
pred<-mm %*% fixef(model)
if (logLink=="e"){
z[x,y]<-exp(pred)
zu[x,y]<-exp(pred+seMultiplier*sqrt(pvar1))
zl[x,y]<-exp(pred-seMultiplier*sqrt(pvar1))
} else if (logLink=="10"){
z[x,y]<-10^(pred)
zu[x,y]<-10^(pred+seMultiplier*sqrt(pvar1))
zl[x,y]<-10^(pred-seMultiplier*sqrt(pvar1))
} else if (logLink=="b"){
z[x,y]<-1/(1+exp(-(pred)))
zu[x,y]<-1/(1+exp(-(pred+seMultiplier*sqrt(pvar1))))
zl[x,y]<-1/(1+exp(-(pred-seMultiplier*sqrt(pvar1))))
} else if (logLink=="n") {
z[x,y]<-mm %*% fixef(model)
zu[x,y]<-z[x,y]+seMultiplier*sqrt(pvar1)
zl[x,y]<-z[x,y]-seMultiplier*sqrt(pvar1)
} else {
stop("Error: the specified log link is not supported")
}
} else {
z[x,y]<-NA
zu[x,y]<-NA
zl[x,y]<-NA
}
}
}
z[(zu-zl)>(quantile(z,0.975,na.rm=TRUE)-min(quantile(z,0.025,na.rm=TRUE)))]<-NA
# ry.cols<-colorRampPalette(c("red","orange","yellow","green"))
# colour<-ry.cols(64)
colour<-brewer.pal(n=11,name="RdYlBu")
nrz <- nrow(z)
ncz <- ncol(z)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, 11)
if (plotUncertainty){
plims=c(min(zl,na.rm=TRUE),max(zu,na.rm=TRUE))
} else {
plims=c(min(z,na.rm=TRUE),max(z,na.rm=TRUE))
}
if(!is.null(zlim)){
plims<-zlim
}
if (plotUncertainty){
persp(e1,e2,zl,theta=45,phi=35,xlab=NA,ylab=NA,zlab=NA,zlim=plims,
col=NA,border="#CCCCCC44",axes=FALSE)
par(new=T)
}
persp(e1,e2,z,theta=45,phi=35,xlab=paste("\n",xlab,sep=""),
ylab=paste("\n",ylab,sep=""),zlab=paste("\n\n",zlab,sep=""),
zlim=plims,ticktype="detailed",nticks=6,col=colour[facetcol])
if(plotUncertainty){
par(new=T)
persp(e1,e2,zu,theta=45,phi=35,xlab=NA,ylab=NA,zlab=NA,zlim=plims,
col=NA,border="#CCCCCC44",axes=FALSE)
}
if (!is.null(outDir)){
dev.off()
}
} else {
stop("Error: a maximum of 2 effects can be specified")
}
}
timnewbold/StatisticalModels documentation built on Aug. 25, 2023, 4:58 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions PlotGLMERContinuous
Documented in PlotGLMERContinuous
PlotGLMERContinuous<-function(model,data,effects,byContEffect=NULL,
otherContEffects=character(0),
otherFactors=character(0),xlab,ylab,
byFactor=NULL,zlab=NULL,outDir=NULL,
logLink="n",plotRug=FALSE,seMultiplier=1.96,
params=list(),add=FALSE,xlim=NULL,ylim=NULL,
zlim=NULL,line.cols=NULL,line.types=NULL,
plotUncertainty=TRUE,nPanels=1,main=NULL,
yDiv=1,transformX=FALSE){
if ((length(effects)>1) && (!is.null(byFactor))){
stop("If plotting by a factor is specified, only one effect can be plotted")
}
if ((!is.null(byFactor)) && (!is.null(byContEffect))){
stop("You cannot specify both a factor and a continuous effect to divide the results")
}
if ((length(effects)==1) && (is.null(byFactor)) && (is.null(byContEffect))){
eval(substitute(newdat<-data.frame(seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length.out=100)),list(x=effects)))
names(newdat)<-effects
for (o in otherContEffects){
newdat[,o] <- median(data[,o])
}
# Add a column to newdat for each of otherFactors
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
newdat[,names(model@frame)[1]] <- 0
mm<-model.matrix(terms(model),newdat)
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
newdat$y<-mm %*% fixef(model)
if (logLink=="e"){
newdat$yplus<-exp(newdat$y+seMultiplier*sqrt(pvar1))
newdat$yminus<-exp(newdat$y-seMultiplier*sqrt(pvar1))
newdat$y<-exp(newdat$y)
} else if (logLink=="10"){
newdat$yplus<-10^(newdat$y+seMultiplier*sqrt(pvar1))
newdat$yminus<-10^(newdat$y-seMultiplier*sqrt(pvar1))
newdat$y<-10^(newdat$y)
} else if (logLink=="n") {
newdat$yplus<-(newdat$y+seMultiplier*sqrt(pvar1))
newdat$yminus<-(newdat$y-seMultiplier*sqrt(pvar1))
newdat$y<-(newdat$y)
} else if (logLink=="b"){
newdat$yplus<-1/(1+exp(-(newdat$y+seMultiplier*sqrt(pvar1))))
newdat$yminus<-1/(1+exp(-(newdat$y-seMultiplier*sqrt(pvar1))))
newdat$y<-1/(1+exp(-(newdat$y)))
} else {
stop("Error: the specified log link is not supported")
}
if (!is.null(outDir)){
png(paste(outDir,"/",effects," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mgp=c(2.5,1,0))
par(mar=c(4,4,1,1))
par(las=1)
par(cex.lab=1.5)
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
if (!transformX){
xlims <- c(min(newdat[,effects]),max(newdat[,effects]))
} else {
xlims <- c(exp(min(newdat[,effects])),exp(max(newdat[,effects])))
}
if (!is.null(xlim)){
xlims <- xlim
}
if (!add){
if (is.null(ylim)){
if(!transformX){
plot(-1e5,-1e5,xlim=xlims,
ylim=c(min(newdat$yminus/yDiv),max(newdat$yplus/yDiv)),
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=c(min(newdat$yminus/yDiv),max(newdat$yplus/yDiv)),
xlab=xlab,ylab=ylab,
log="x",main=main)
}
} else {
if (!transformX){
plot(-1e5,-1e5,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main,
log="x")
}
}
}
line.col<-ifelse(is.null(line.cols),"#000000",line.cols[1])
if(plotUncertainty){
if (transformX){
X.Vec <- exp(c(newdat[,effects], max(newdat[,effects]),
rev(newdat[,effects]), min(newdat[,effects])))
} else {
X.Vec <- c(newdat[,effects], max(newdat[,effects]),
rev(newdat[,effects]), min(newdat[,effects]))
}
Y.Vec <- c(newdat$yminus/yDiv, tail(newdat$yplus/yDiv, 1), rev(newdat$yplus/yDiv), (newdat$yminus/yDiv)[1])
polygon(X.Vec, Y.Vec, col = paste(line.col,"33",sep=""), border = NA)
}
if(transformX){
points(exp(newdat[,effects]),
newdat$y/yDiv,type="l",col=line.col)
} else {
points(newdat[,effects],
newdat$y/yDiv,type="l",col=line.col)
}
if(plotRug){
if(transformX){
rug(eval(substitute(exp(data$x),list(x=effects[1]))))
} else {
rug(eval(substitute(data$x,list(x=effects[1]))))
}
}
if (!is.null(outDir)){
dev.off()
}
} else if ((length(effects)==1) && (!is.null(byFactor))){
if((!is.null(line.cols)) & (length(line.cols)!=length(levels(model@frame[,byFactor])))){
stop("Specified colours must be of the same number as the number of factor levels")
}
if((!is.null(line.types)) & (length(line.types)!=length(levels(model@frame[,byFactor])))){
stop("Specified line types must be of the same number as the number of factor levels")
}
eval(substitute(newdat<-data.frame(seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length.out=100)),list(x=effects)))
names(newdat)<-effects
# Add a column to newdat for each of otherFactors
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
newdat[,names(model@frame)[1]] <- 0
y<-data.frame(x=newdat[,1])
yplus<-data.frame(x=newdat[,1])
yminus<-data.frame(x=newdat[,1])
for (l in levels(model@frame[,byFactor])){
newdat[,byFactor]<-factor(l,levels=levels(model@frame[,byFactor]))
for (o in otherContEffects){
newdat[,o] <- median(data[,o][(data[,byFactor]==l)])
}
mm<-model.matrix(terms(model),newdat)
if(ncol(mm)>length(fixef(model))){
mm <- mm[,-which(!(names(mm[1,]) %in% names(fixef(model))))]
}
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
y[,l]<-mm %*% fixef(model)
if (logLink=="e"){
yplus[,l]<-exp(y[,l]+seMultiplier*sqrt(pvar1))
yminus[,l]<-exp(y[,l]-seMultiplier*sqrt(pvar1))
y[,l]<-exp(y[,l])
} else if (logLink=="10"){
yplus[,l]<-10^(y[,l]+seMultiplier*sqrt(pvar1))
yminus[,l]<-10^(y[,l]-seMultiplier*sqrt(pvar1))
y[,l]<-10^(y[,l])
} else if (logLink=="b"){
yplus[,l]<-1/(1+exp(-(y[,l]+seMultiplier*sqrt(pvar1))))
yminus[,l]<-1/(1+exp(-(y[,l]-seMultiplier*sqrt(pvar1))))
y[,l]<-1/(1+exp(-(y[,l])))
} else if (logLink=="n") {
yplus[,l]<-(y[,l]+seMultiplier*sqrt(pvar1))
yminus[,l]<-(y[,l]-seMultiplier*sqrt(pvar1))
} else {
stop("Error: the specified log link is not supported")
}
try(y[y$x>quantile(data[data[,byFactor]==l,][,effects[1]],0.975),][,l]<-NA,silent=TRUE)
try(yplus[yplus$x>quantile(data[data[,byFactor]==l,][,effects[1]],0.975),][,l]<-NA,silent=TRUE)
try(yminus[yminus$x>quantile(data[data[,byFactor]==l,][,effects[1]],0.975),][,l]<-NA,silent=TRUE)
try(y[y$x<quantile(data[data[,byFactor]==l,][,effects[1]],0.025),][,l]<-NA,silent=TRUE)
try(yplus[yplus$x<quantile(data[data[,byFactor]==l,][,effects[1]],0.025),][,l]<-NA,silent=TRUE)
try(yminus[yminus$x<quantile(data[data[,byFactor]==l,][,effects[1]],0.025),][,l]<-NA,silent=TRUE)
}
if (!is.null(outDir)){
png(paste(outDir,"/",effects," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mgp=c(2.5,1,0))
par(mar=c(4,4,1,1))
par(las=1)
par(cex.lab=1.5)
if (nPanels>1) par(mfrow=c(1,nPanels))
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
allLevels<-1:length(levels(model@frame[,byFactor]))
nEffectsPerPlot<-ceiling(length(allLevels)/nPanels)
effectsPerPlot<-split(allLevels,ceiling(allLevels/nEffectsPerPlot))
if(!is.null(line.cols)){
cols<-split(line.cols,ceiling(allLevels/nEffectsPerPlot))
}
if(!is.null(line.types)){
line.types<-split(line.types,ceiling(allLevels/nEffectsPerPlot))
}
for (p in 1:nPanels){
if (!add){
options(scipen=0)
if(!transformX){
xlims<-c(min(y$x),max(y$x))
} else {
xlims<-c(exp(min(y$x)),exp(max(y$x)))
}
if(!is.null(xlim)){
xlims<-xlim
}
if (is.null(ylim)){
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,log="x",main=main)
}
} else {
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,log="x",main=main)
}
}
}
if (is.null(line.cols)){
cols<-list()
cols[1][[1]]<-rep("#000000",length(levels(model@frame[,byFactor])))
} else {
}
o<-c(o,setdiff(1:length(levels(model@frame[,byFactor])),o))
if (plotUncertainty){
i<-1
for (l in levels(model@frame[,byFactor])[effectsPerPlot[p][[1]]]){
y.temp<-y[,c('x',l)]
y.temp<-na.omit(y.temp)
yplus.temp<-yplus[,c('x',l)]
yplus.temp<-na.omit(yplus.temp)
yminus.temp<-yminus[,c('x',l)]
yminus.temp<-na.omit(yminus.temp)
if(transformX){
X.Vec <- exp(c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x)))
} else {
X.Vec <- c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x))
}
Y.Vec <- c((yminus.temp/yDiv)[,l], tail((yplus.temp/yDiv)[,l], 1), rev((yplus.temp/yDiv)[,l]), (yminus.temp/yDiv)[,l][1])
polygon(X.Vec, Y.Vec, col = paste(cols[p][[1]][i],"33",sep=""), border = NA)
i<-i+1
}
}
i<-1
for (l in levels(model@frame[,byFactor])[effectsPerPlot[p][[1]]]){
if(transformX){
points(exp(y$x),(y/yDiv)[,l],type="l",col=cols[p][[1]][i],lty=ifelse(is.null(line.types),1,line.types[p][[1]][i]))
if (plotRug){
rug(exp(data[data[,byFactor]==l,][,effects[1]]),col=cols[p][[1]][i])
}
} else {
points(y$x,(y/yDiv)[,l],type="l",col=cols[p][[1]][i],lty=ifelse(is.null(line.types),1,line.types[p][[1]][i]))
if (plotRug){
rug(data[data[,byFactor]==l,][,effects[1]],col=cols[p][[1]][i])
}
}
i<-i+1
}
}
} else if ((length(effects)==1) && (!is.null(byContEffect))){
if((!is.null(line.cols)) & (length(line.cols)!=3)){
stop("Specified colours must be of the same number as the number of factor levels")
}
if((!is.null(line.types)) & (length(line.types)!=3)){
stop("Specified line types must be of the same number as the number of factor levels")
}
eval(substitute(newdat<-data.frame(seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length.out=100)),
list(x=effects)))
names(newdat)<-effects
# Add a column to newdat for each of otherFactors
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
if (class(model)=="glmmadmb"){
newdat[,names(model$frame)[1]] <- 0
} else {
newdat[,names(model@frame)[1]] <- 0
}
y<-data.frame(x=newdat[,1])
yplus<-data.frame(x=newdat[,1])
yminus<-data.frame(x=newdat[,1])
byContEffectVals <- c(quantile(data[,byContEffect],probs = 0.025,na.rm=TRUE),
median(data[,byContEffect],na.rm=TRUE),
quantile(data[,byContEffect],probs = 0.975,na.rm=TRUE))
binWidth <- 0.1
prbs <- c(binWidth,0.5-(0.5*binWidth),0.5+(0.5*binWidth),1-binWidth)
brks <- c(-Inf,quantile(data[,byContEffect],
probs = prbs,na.rm = TRUE),Inf)
brks <- brks + seq_along(brks) * 10 * .Machine$double.eps
data$ByContEffectTerciles <- cut(data[,byContEffect],
breaks = brks,
labels = c(1,"NA1",2,"NA2",3))
emptyBin <- any(unlist(lapply(as.list(c(1,2,3)),function(x){
return(length(which(data$ByContEffectTerciles==x))==0)
})))
while(emptyBin){
binWidth <- binWidth + 0.01
prbs <- c(binWidth,0.5-(0.5*binWidth),0.5+(0.5*binWidth),1-binWidth)
brks <- c(-Inf,quantile(data[,byContEffect],
probs = prbs,na.rm = TRUE),Inf) +
seq_along(brks) * 10 * .Machine$double.eps
data$ByContEffectTerciles <- cut(data[,byContEffect],
breaks = brks,
labels = c(1,"NA1",2,"NA2",3))
emptyBin <- any(unlist(lapply(as.list(c(1,2,3)),function(x){
return(length(which(data$ByContEffectTerciles==x))==0)
})))
}
for (i in 1:length(byContEffectVals)){
newdat[,byContEffect] <- byContEffectVals[i]
for (o in otherContEffects){
newdat[,o] <- median(data[,o][(data$ByContEffectTerciles==i)],na.rm=TRUE)
}
mm<-model.matrix(terms(model),newdat)
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
y[,paste(i)]<-mm %*% fixef(model)
if (logLink=="e"){
yplus[,paste(i)]<-exp(y[,paste(i)]+seMultiplier*sqrt(pvar1))
yminus[,paste(i)]<-exp(y[,paste(i)]-seMultiplier*sqrt(pvar1))
y[,paste(i)]<-exp(y[,paste(i)])
} else if (logLink=="10"){
yplus[,paste(i)]<-10^(y[,paste(i)]+seMultiplier*sqrt(pvar1))
yminus[,paste(i)]<-10^(y[,paste(i)]-seMultiplier*sqrt(pvar1))
y[,paste(i)]<-10^(y[,paste(i)])
} else if (logLink=="b"){
yplus[,paste(i)]<-1/(1+exp(-(y[,paste(i)]+seMultiplier*sqrt(pvar1))))
yminus[,paste(i)]<-1/(1+exp(-(y[,paste(i)]-seMultiplier*sqrt(pvar1))))
y[,paste(i)]<-1/(1+exp(-(y[,paste(i)])))
} else if (logLink=="n") {
yplus[,paste(i)]<-(y[,paste(i)]+seMultiplier*sqrt(pvar1))
yminus[,paste(i)]<-(y[,paste(i)]-seMultiplier*sqrt(pvar1))
} else {
stop("Error: the specified log link is not supported")
}
try(y[y$x>quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.975),][,paste(i)]<-NA,silent=TRUE)
try(yplus[yplus$x>quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.975),][,paste(i)]<-NA,silent=TRUE)
try(yminus[yminus$x>quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.975),][,paste(i)]<-NA,silent=TRUE)
try(y[y$x<quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.025),][,paste(i)]<-NA,silent=TRUE)
try(yplus[yplus$x<quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.025),][,paste(i)]<-NA,silent=TRUE)
try(yminus[yminus$x<quantile(data[data$ByContEffectTerciles==i,][,effects[1]],0.025),][,paste(i)]<-NA,silent=TRUE)
}
if (!is.null(outDir)){
png(paste(outDir,"/",effects," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mgp=c(2.5,1,0))
par(mar=c(4,4,1,1))
par(las=1)
par(cex.lab=1.5)
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
if (!add){
options(scipen=0)
if(!transformX){
xlims<-c(min(y$x),max(y$x))
} else {
xlims<-c(exp(min(y$x)),exp(max(y$x)))
}
if(!is.null(xlim)){
xlims<-xlim
}
if (is.null(ylim)){
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=c(min((yminus/yDiv)[,-1],na.rm=TRUE),
max((yplus/yDiv)[,-1],na.rm=TRUE)),
xlab=xlab,ylab=ylab,log="x",main=main)
}
} else {
if(!transformX){
plot(-99999,-99999,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,main=main)
} else {
plot(9e99,9e99,xlim=xlims,
ylim=ylim,
xlab=xlab,ylab=ylab,log="x",main=main)
}
}
}
if (is.null(line.cols)){
cols<-rep("#000000",length(byContEffectVals))
} else {
cols <- line.cols
}
if (plotUncertainty){
for (i in 1:length(byContEffectVals)){
y.temp<-y[,c('x',paste(i))]
y.temp<-na.omit(y.temp)
yplus.temp<-yplus[,c('x',paste(i))]
yplus.temp<-na.omit(yplus.temp)
yminus.temp<-yminus[,c('x',paste(i))]
yminus.temp<-na.omit(yminus.temp)
if(transformX){
X.Vec <- exp(c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x)))
} else {
X.Vec <- c(y.temp$x, max(y.temp$x), rev(y.temp$x), min(y.temp$x))
}
Y.Vec <- c((yminus.temp/yDiv)[,paste(i)], tail((yplus.temp/yDiv)[,paste(i)], 1), rev((yplus.temp/yDiv)[,paste(i)]), (yminus.temp/yDiv)[,paste(i)][1])
polygon(X.Vec, Y.Vec, col = paste(cols[i],"33",sep=""), border = NA)
i<-i+1
}
}
i<-1
for (i in 1:length(byContEffectVals)){
if(transformX){
points(exp(y$x),(y/yDiv)[,paste(i)],type="l",col=cols[i],lty=ifelse(is.null(line.types),1,line.types[i]))
# if (plotRug){
# rug(exp(data[data[,byFactor]==l,][,effects[1]]),col=cols[p][[1]][i])
# }
} else {
points(y$x,(y/yDiv)[,paste(i)],type="l",col=cols[i],lty=ifelse(is.null(line.types),1,line.types[i]))
# if (plotRug){
# rug(data[data[,byFactor]==l,][,effects[1]],col=cols[p][[1]][i])
# }
}
}
} else if (length(effects==2)){
if (!is.null(outDir)){
png(paste(outDir,"/",effects[1]," and ",effects[2]," effect_",names(model@frame)[1],".png",
sep=""),width=22.86/2.54,height=12.57/2.54,units="in",res=1200)
}
par(mar=c(2,2,2,2))
par(mgp=c(2.5,1,0))
par(las=0)
par(cex.lab=1.5)
for (p in names(params)){
eval(parse(text=gsub("x",p,"par(x=params$x)")))
}
eval(substitute(e1<-seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length=30),list(x=effects[1])))
eval(substitute(e2<-seq(from=min(data$x,na.rm=TRUE),to=max(data$x,na.rm=TRUE),length=30),list(x=effects[2])))
z<-array(dim=c(30,30))
zu<-array(dim=c(30,30))
zl<-array(dim=c(30,30))
gx<-rep(e1,30)
gy<-rep(e2,rep(30,30))
eval(substitute(tf<-matrix(exclude.too.far(gx,gy,data$x,data$y,0.2),
byrow=F,ncol=30),list(x=effects[1],y=effects[2])))
for (x in 1:length(e1)){
for (y in 1:length(e2)){
v1<-e1[x]
v2<-e2[y]
newdat<-data.frame(v1,v2)
names(newdat)<-effects
# Add a column to newdat for each other continuous effect
for (o in otherContEffects){
newdat[,o] <- median(data[,o])
}
# Add a column to newdat for each other factor
for(i in 1:length(otherFactors)) {
col <- names(otherFactors)[i]
newdat[,col] <- factor(otherFactors[i], levels=levels(data[,col]))
}
newdat[,names(model@frame)[1]] <- 0
mm<-model.matrix(terms(model),newdat)
pvar1 <- diag(mm %*% base::tcrossprod(as.matrix(vcov(model)),mm))
if(tf[x,y]==FALSE){
pred<-mm %*% fixef(model)
if (logLink=="e"){
z[x,y]<-exp(pred)
zu[x,y]<-exp(pred+seMultiplier*sqrt(pvar1))
zl[x,y]<-exp(pred-seMultiplier*sqrt(pvar1))
} else if (logLink=="10"){
z[x,y]<-10^(pred)
zu[x,y]<-10^(pred+seMultiplier*sqrt(pvar1))
zl[x,y]<-10^(pred-seMultiplier*sqrt(pvar1))
} else if (logLink=="b"){
z[x,y]<-1/(1+exp(-(pred)))
zu[x,y]<-1/(1+exp(-(pred+seMultiplier*sqrt(pvar1))))
zl[x,y]<-1/(1+exp(-(pred-seMultiplier*sqrt(pvar1))))
} else if (logLink=="n") {
z[x,y]<-mm %*% fixef(model)
zu[x,y]<-z[x,y]+seMultiplier*sqrt(pvar1)
zl[x,y]<-z[x,y]-seMultiplier*sqrt(pvar1)
} else {
stop("Error: the specified log link is not supported")
}
} else {
z[x,y]<-NA
zu[x,y]<-NA
zl[x,y]<-NA
}
}
}
z[(zu-zl)>(quantile(z,0.975,na.rm=TRUE)-min(quantile(z,0.025,na.rm=TRUE)))]<-NA
# ry.cols<-colorRampPalette(c("red","orange","yellow","green"))
# colour<-ry.cols(64)
colour<-brewer.pal(n=11,name="RdYlBu")
nrz <- nrow(z)
ncz <- ncol(z)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
facetcol <- cut(zfacet, 11)
if (plotUncertainty){
plims=c(min(zl,na.rm=TRUE),max(zu,na.rm=TRUE))
} else {
plims=c(min(z,na.rm=TRUE),max(z,na.rm=TRUE))
}
if(!is.null(zlim)){
plims<-zlim
}
if (plotUncertainty){
persp(e1,e2,zl,theta=45,phi=35,xlab=NA,ylab=NA,zlab=NA,zlim=plims,
col=NA,border="#CCCCCC44",axes=FALSE)
par(new=T)
}
persp(e1,e2,z,theta=45,phi=35,xlab=paste("\n",xlab,sep=""),
ylab=paste("\n",ylab,sep=""),zlab=paste("\n\n",zlab,sep=""),
zlim=plims,ticktype="detailed",nticks=6,col=colour[facetcol])
if(plotUncertainty){
par(new=T)
persp(e1,e2,zu,theta=45,phi=35,xlab=NA,ylab=NA,zlab=NA,zlim=plims,
col=NA,border="#CCCCCC44",axes=FALSE)
}
if (!is.null(outDir)){
dev.off()
}
} else {
stop("Error: a maximum of 2 effects can be specified")
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.