# Dynamic Panel Data Models (Arellano-Bond, etc.) # for fixed effects in short-T panels # # Helper functions # # Christopher Adolph faculty.washington.edu/cadolph # 8 June 2015 labelBottomK <- function(l,x,k) { labels <- rep("",length(l)) ox <- rank(x) for (i in 1:length(l)) { if (ox[i]<=k) labels[i] <- l[i] } labels } labelTopK <- function(l,x,k) { labels <- rep("",length(l)) ox <- rank(x) for (i in 1:length(l)) { if (ox[i]>(length(l)-k)) labels[i] <- l[i] } labels } labelTopBottomK <- function(l,x,k) { labels <- rep("",length(l)) ox <- rank(x) for (i in 1:length(l)) { if ( (ox[i]>(length(l)-k))|(ox[i]<=k)) labels[i] <- l[i] } labels } # Graph expected values, baseline, and first diffs for all 8 models using tile... cigLineplots <- function(EV, EVbase, FD, RR, limitsEV, limitsFD, limitsRR, initialY, main, file) { # Get 3 nice colors for traces require(RColorBrewer) col <- brewer.pal(4,"Dark2") periods <- 1:3 # Set up lineplot trace of EV packpc EVtrace <- lineplot(x=c(0, periods), y=c(initialY, EV$pe), lower=c(initialY, EV$lower), upper=c(initialY, EV$upper), col=col[1], plot=1) EVptrace <- pointsTile(x=c(0, periods), y=c(initialY, EV$pe), pch=16, col=col[1], plot=1) BASEtrace <- lineplot(x=c(0, periods), y=c(initialY, EVbase$pe), lower=c(initialY, EVbase$lower), upper=c(initialY, EVbase$upper), ci=list(mark="dashed"), col=col[2], plot=1) BASEptrace <- pointsTile(x=c(0, periods), y=c(initialY, EVbase$pe), pch=16, col=col[2], plot=1) FDtrace <- lineplot(x=c(0, periods), y=c(0, FD$pe), lower=c(0, FD$lower), upper=c(0, FD$upper), col=col[3], plot=2) FDptrace <- pointsTile(x=c(0, periods), y=c(0, FD$pe), pch=16, col=col[3], plot=2) RRtrace <- lineplot(x=c(0, periods), y=c(0, 100*(RR$pe-1)), lower=c(0, 100*(RR$lower-1)), upper=c(0, 100*(RR$upper-1)), col=col[4], plot=3) RRptrace <- pointsTile(x=c(0, periods), y=c(0, 100*(RR$pe-1)), pch=16, col=col[4], plot=3) # Set up text labels of lines text1 <- textTile(x=3.9, y=EV$pe[length(EV$pe)], labels="+$.60/pack", col=col[1], cex=0.75, plot=1) text2 <- textTile(x=3.85, y=EVbase$pe[length(EVbase$pe)], labels="baseline", col=col[2], cex=0.75, plot=1) text3 <- textTile(x=3.85, y=FD$pe[length(FD$pe)], labels="change\nin packs", col=col[3], cex=0.75, plot=2) text4 <- textTile(x=3.85, y=100*(RR$pe[length(RR$pe)] -1), labels="percent\nchange", col=col[4], cex=0.75, plot=3) # Set up baseline: for first difference, this is 0 baselineFD <- linesTile(x=c(-1,5), y=c(0,0), plot=2) # Set up baseline: for percent change, this is 0 baselineRR <- linesTile(x=c(-1,5), y=c(0, 0), plot=3) # Plot all traces using tile tile(EVtrace, BASEtrace, FDtrace, RRtrace, EVptrace, BASEptrace, FDptrace, RRptrace, text1, text2, text3, text4, baselineFD, baselineRR, limits=rbind(c(-0.25,4.8,limitsEV), c(-0.25,4.8,limitsFD), c(-0.25,4.8,limitsRR)), xaxis=list(at=c(0,periods), labels=c("t", "t+1", "t+2", "t+3")), yaxis=list(label.loc=-0.5, major=FALSE), xaxistitle=list(labels="Forecast Year"), yaxistitle=list(labels1="Expected Packs per capita", labels2="Change, Packs per capita", labels3="%Change, Packs per capita"), maintitle=list(labels=main), gridlines=list(type="y"), height=list(maintitle=4), width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5), output=list(file=file,width=10) ) } collectUnitSims <- function(x,newsims) { require(abind) if (is.null(newsims$simulates)) stop("newsims must be created using simulates=TRUE") if (is.null(x$units)) x$units <- list() x$units$pe <- rbind(x$units$pe, t(newsims$pe)) x$units$lower <- rbind(x$units$lower, t(newsims$lower)) x$units$upper <- rbind(x$units$upper, t(newsims$upper)) if (is.null(x$simulates)) x$simulates <- list() x$simulates$units <- abind(x$simulates$units, newsims$simulates, along=3) x } aggregateUnitSims <- function(x,fn=mean,ci=0.95,...) { agg <- deparse(substitute(fn)) x$simulates[[agg]] <- apply(x$simulates$units, 1:2, fn, ...) x[[agg]]$pe <- apply(x$simulates[[agg]], 2, mean) x[[agg]]$lower <- apply(x$simulates[[agg]], 2, quantile, probs=(1-ci)/2) x[[agg]]$upper <- apply(x$simulates[[agg]], 2, quantile, probs=(1+ci)/2) x } cigUnitLineplots <- function(FD, RR, compFD=NULL, compRR=NULL, periods=NULL, units=NULL, showUnits=TRUE, showMean=TRUE, avg=NULL, avgLab=NULL, compLab=NULL, unitLabX=NULL, labSet=c("tb", "t", "b"), labK=3, initialZero=FALSE, CI=FALSE, compCI=FALSE, limitsFD=c(NA,NA,NA,NA), limitsRR=c(NA,NA,NA,NA), col=c("black","black"), dots=FALSE, RRasPD=FALSE, xat=NA, yat=NA, plottitles=NULL, xtitles=NULL, ytitles=NULL, file=NULL ) { require(tile) # create list vector and start counter of traces added to list vector collectedTraces <- vector(mode="list") ctr <- 0 # Force compFD and compRR to vectors if (!is.null(compFD)) compFD <- lapply(compFD, function(x) {as.vector(x)}) if (!is.null(compRR)) compRR <- lapply(compRR, function(x) {as.vector(x)}) #transform RR if requested if (RRasPD) RR <- rapply(RR, function(x) {100*(x-1)}, classes=c("numeric","matrix"), how="replace") if (RRasPD&!is.null(compRR)) compRR <- rapply(compRR, function(x) {100*(x-1)}, classes=c("numeric","matrix"), how="replace") nunit <- nrow(FD$units$pe) nper <- ncol(FD$units$pe) if (is.null(periods)) periods <- 1:nper # function to insert zeros to start of all elements of a list, recursively insertZeros <- function(x,z) { x <- rapply(x, function(x) {c(z,x)}, classes=c("numeric"), how="replace") rapply(x, function(x) {cbind(rep(z, nrow(x)), x)}, classes="matrix", how="replace") } # Add initial period zeros if requested if (initialZero) { periods <- c(periods[1]-1, periods) nper <- nper+1 fd0 <- 0 FD <- insertZeros(FD, fd0) if (RRasPD) {rr0 <- 0} else {rr0 <- 1} RR <- insertZeros(RR, rr0) if (!is.null(compFD)) compFD <- insertZeros(compFD, fd0) if (!is.null(compRR)) compRR <- insertZeros(compRR, fd0) } if (showUnits) { FDunitTrace <- polylinesTile(y=t(FD$units$pe), x=matrix(periods, nrow=nper, ncol=nunit), id=matrix(1:nunit, nrow=nper, ncol=nunit, byrow=TRUE), col="gray70", plot=1) ctr <- ctr+1 collectedTraces[[ctr]] <- FDunitTrace RRunitTrace <- polylinesTile(y=t(RR$units$pe), x=matrix(periods, nrow=nper, ncol=nunit), id=matrix(1:nunit, nrow=nper, ncol=nunit, byrow=TRUE), col="gray75", plot=2) ctr <- ctr+1 collectedTraces[[ctr]] <- RRunitTrace if (!is.null(units)) { if (is.null(labSet)) { fuLab <- rep("", nunit) } else { labSet <- labSet[1] if (labSet=="tb") { fuLab <- labelTopBottomK(units, FD$units$pe[,ncol(FD$units$pe)], labK) } if (labSet=="t") { fuLab <- labelTopK(units, FD$units$pe[,ncol(FD$units$pe)], labK) } if (labSet=="b") { fuLab <- labelBottomK(units, FD$units$pe[,ncol(FD$units$pe)], labK) } } } FDunitLabels <- textTile(x=rep(unitLabX, nunit), y=FD$units$pe[,ncol(FD$units$pe)], labels=fuLab, cex=0.8, col="gray75", plot=1) ctr <- ctr+1 collectedTraces[[ctr]] <- FDunitLabels if (!is.null(units)) { if (is.null(labSet)) { fuLab <- rep("", nunit) } else { labSet <- labSet[1] if (labSet=="tb") { fuLab <- labelTopBottomK(units, RR$units$pe[,ncol(RR$units$pe)], labK) } if (labSet=="t") { fuLab <- labelTopK(units, RR$units$pe[,ncol(RR$units$pe)], labK) } if (labSet=="b") { fuLab <- labelBottomK(units, RR$units$pe[,ncol(RR$units$pe)], labK) } } } RRunitLabels <- textTile(x=rep(unitLabX, nunit), y=RR$units$pe[,ncol(RR$units$pe)], labels=fuLab, cex=0.8, col="gray75", plot=2) ctr <- ctr+1 collectedTraces[[ctr]] <- RRunitLabels } if (showMean) { if (CI) { FDmeanTrace <- lineplot(y=FD[[avg]]$pe, x=periods, lower=FD[[avg]]$lower, upper=FD[[avg]]$upper, col=col[1], lex=3, plot=1) } else { FDmeanTrace <- lineplot(y=FD[[avg]]$pe, x=periods, col=col[1], lex=3, plot=1) } ctr <- ctr+1 collectedTraces[[ctr]] <- FDmeanTrace if (CI) { RRmeanTrace <- lineplot(y=RR[[avg]]$pe, x=periods, lower=RR[[avg]]$lower, upper=RR[[avg]]$upper, col=col[1], lex=3, plot=2) } else { RRmeanTrace <- lineplot(y=RR[[avg]]$pe, x=periods, col=col[1], lex=3, plot=2) } ctr <- ctr+1 collectedTraces[[ctr]] <- RRmeanTrace if (dots) { FDmeanpTrace <- pointsTile(y=FD[[avg]]$pe, x=periods, col=col[1], pch=16, size=.6, alpha=1, plot=1) ctr <- ctr+1 collectedTraces[[ctr]] <- FDmeanpTrace RRmeanpTrace <- pointsTile(y=RR[[avg]]$pe, x=periods, col=col[1], pch=16, size=.6, alpha=1, plot=2) ctr <- ctr+1 collectedTraces[[ctr]] <- RRmeanpTrace } FDmeanLabel <- textTile(y=FD[[avg]]$pe[nper], x=unitLabX, labels=avgLab, cex=0.8, col=col[1], plot=1) ctr <- ctr+1 collectedTraces[[ctr]] <- FDmeanLabel RRmeanLabel <- textTile(y=RR[[avg]]$pe[nper], x=unitLabX, labels=avgLab, cex=0.8, col=col[1], plot=2) ctr <- ctr+1 collectedTraces[[ctr]] <- RRmeanLabel } if (!is.null(compFD)) { if (compCI) { FDcompTrace <- lineplot(y=compFD$pe, x=periods, lower=compFD$lower, upper=compFD$upper, ci=list(mark="dashed"), col=col[2], lex=3, plot=1) } else { FDcompTrace <- lineplot(y=compFD$pe, x=periods, col=col[2], lex=3, plot=1) } ctr <- ctr+1 collectedTraces[[ctr]] <- FDcompTrace if (dots) { FDcomppTrace <- pointsTile(y=compFD$pe, x=periods, col=col[2], pch=16, size=.6, alpha=1, plot=1) ctr <- ctr+1 collectedTraces[[ctr]] <- FDcomppTrace } FDcompLabel <- textTile(y=compFD$pe[nper], x=unitLabX, labels=compLab, cex=0.8, col=col[2], plot=1) ctr <- ctr+1 collectedTraces[[ctr]] <- FDcompLabel } if (!is.null(compRR)) { if (compCI) { RRcompTrace <- lineplot(y=compRR$pe, x=periods, lower=compRR$lower, upper=compRR$upper, ci=list(mark="dashed"), col=col[2], lex=3, plot=2) } else { RRcompTrace <- lineplot(y=compRR$pe, x=periods, col=col[2], lex=3, plot=2) } ctr <- ctr+1 collectedTraces[[ctr]] <- RRcompTrace if (dots) { RRcomppTrace <- pointsTile(y=compRR$pe, x=periods, col=col[2], pch=16, size=.6, alpha=1, plot=2) ctr <- ctr+1 collectedTraces[[ctr]] <- RRcomppTrace } RRcompLabel <- textTile(y=compRR$pe[nper], x=unitLabX, labels=compLab, cex=0.8, col=col[2], plot=2) ctr <- ctr+1 collectedTraces[[ctr]] <- RRcompLabel } # Set up baseline: for first difference, this is 0 baselineFD <- linesTile(x=limitsFD[1:2], y=c(0,0), plot=1) # Set up baseline: for percent change, this is 0; # for relative risk, this is 1 if (RRasPD) { baselineRR <- linesTile(x=limitsRR[1:2], y=c(0, 0), plot=2) } else { baselineRR <- linesTile(x=limitsRR[1:2], y=c(1, 1), plot=2) } # Plot all traces using tile tile(baselineFD, baselineRR, collectedTraces, limits=rbind(limitsFD,limitsRR), xaxis=list(at=xat),#, labels=c("t", "t+1", "t+2", "t+3")), yaxis=list(label.loc=-0.5, major=FALSE, yat=yat), xaxistitle=list(labels=xtitles), yaxistitle=list(labels=ytitles), plottitle=list(labels=plottitles), gridlines=list(type="y"), #height=list(maintitle=4), width=list(null=5,yaxistitle=4,yaxis.labelspace=-0.5), output=list(file=file,width=10) ) } #add in aggregate #mark aggregate with dashed line #sizeSig <- 0.57*1.3 #sizeNSig <- 0.57*1.15 cigUnitDotplot <- function(x, units, period=1, showRank=FALSE, showSelected=NULL, reverse=FALSE, pchSig=15, pchNSig=0, sizeSig=0.741, sizeNSig=0.741, col="black", RRasPD=FALSE, limits=NULL, labspace=0.35, at=NA, vertmark=NA, entryheight=0.05, plottitle=NULL, xtitle=NULL, toptitle=NULL, file=NULL ) { # convert to % change from relative risk? if (RRasPD) x <- rapply(x, function(x) {100*(x-1)}, classes=c("numeric","matrix"), how="replace") # reverse? if (reverse) reverse <- -1 else reverse <- 1 # Order on results ord <- order(as.matrix(x$pe)[,period]) if (showRank) { labels <- paste(1:length(ord), ". ", units[ord], sep="") } else { labels <- units[ord] } pe <- reverse*as.matrix(x$pe)[ord, period] lower <- reverse*as.matrix(x$lower)[ord, period] upper <- reverse*as.matrix(x$upper)[ord, period] sig <- (lower*upper)>0 symb <- pchSig*sig + pchNSig*!sig sizethis <- sizeSig*sig + sizeNSig*!sig # Create limits if (is.null(limits)) limits <- c(min(pe, lower, upper, vertmark, na.rm=TRUE), max(pe, lower, upper, vertmark, na.rm=TRUE)) # Show only selected ranks if (!is.null(showSelected)) { pe <- pe[showSelected] lower <- lower[showSelected] upper <- upper[showSelected] labels <- labels[showSelected] } # Make this trace, with options unitTrace <- ropeladder(labels=labels, x=pe, lower=lower, upper=upper, size=sizethis, entryheight=entryheight, pch=symb, col=col, #fill=bg[i], group=1, layer=i, plot=1) # Make shaded row traces labspace <- labspace*(limits[2]-limits[1]) pctr <- 0 collectedPoly <- vector("list", ceiling(length(labels)/2)) for (j in 1:length(labels)) { if (j%%2) { pctr <- pctr+1 yup <- 1-1/length(labels)*(j-1) ylo <- 1-1/length(labels)*j collectedPoly[[pctr]] <- polygonTile(x=c(limits[1]-labspace, limits[1]-labspace, limits[2], limits[2]), y=c(yup, ylo, ylo, yup), col=NA, fill="gray80", clip="off", plot=1) } } # Make annotations vertmark <- linesTile(x=c(0,0), y=c(0,1), plot=1) # Send to tile tile(vertmark, unitTrace, collectedPoly, limits=c(limits,limits), topaxis=list(add=TRUE, at=at, major=FALSE), xaxis=list(major=FALSE, at=at), xaxistitle=list(labels=xtitle, y=0.25), topaxistitle=list(labels=toptitle, y=0.75), gridlines=list(type1="xt"), width=list(plot=1), output=list(outfile=file, width=8) ) } # Summary function for small T PGMM # C. Lanfear - 6-8-2015 # x is a pgmm object produced by pgmm() summary2 <- function(x){ if (x$args$transformation=="d") { l <- 1 } else { l <- 0 } if (length(x$args$namest) < 4-l) { order<-1 } else { order <- 2 } if (length(x$args$namest) > 2-l) { MTest<-mtest(x,order=order,vcovHC(x)) } else { MTest<-"MTest requires T > 3" } Sargan<-sargan(x) std.err <- sqrt(diag(vcovHC(x))) b <- coef(x) z <- b/std.err p <- 2 * pnorm(abs(z), lower.tail = FALSE) Summary <- cbind(round(b,4), round(std.err,4), round(z,3), round(p,3)) colnames(Summary) <- c("b","se","z","p") ChiSq<-plm:::wald(x, "time", vcovHC(x)) list<-list(MTest,Sargan,Summary,ChiSq) names(list)<-c("MTest","Sargan","Summary","ChiSq") list }