#coefficient plot the easy way
varnames <- c("int", "Sex, F", "age")
mod1coefse <- matrix(scan(), ncol = 2, byrow=TRUE)
-2.1 0.3   1.1 0.9  NA NA

mod1coefse

mod2coefse <- matrix(scan(), ncol = 2, byrow=TRUE)
-2.3 1.3   2.1 1.2  1.3 3.3

mod2coefse

mod3coefse <- matrix(scan(), ncol = 2, byrow=TRUE)
-3.0 1.1
1.9 0.9
2.1 0.9

mod3coefse

mod1df <- data.frame(Coefficient = varnames, Value = mod1coefse[,1], SE = mod1coefse[,2], Model = "Baseline")
mod2df <- data.frame(Coefficient = varnames, Value = mod2coefse[,1], SE = mod2coefse[,2], Model = "Plus 1")
mod3df <- data.frame(Coefficient = varnames, Value = mod3coefse[,1], SE = mod3coefse[,2], Model = "Plus 2")
allModelFrame <- data.frame(rbind(mod1df, mod2df, mod3df))

interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
zp1 <- ggplot(allModelFrame, aes(colour = Model))
zp1 <- zp1 + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
zp1 <- zp1 + geom_linerange(aes(x = Coefficient, ymin = Value - SE*interval1,
                                ymax = Value + SE*interval1),
                            lwd = 1, position = position_dodge(width = 1/2))
zp1 <- zp1 + geom_pointrange(aes(x = Coefficient, y = Value, ymin = Value - SE*interval2,
                                 ymax = Value + SE*interval2),
                             lwd = 1/2, position = position_dodge(width = 1/2),
                             shape = 21, fill = "WHITE")
zp1 <- zp1 + coord_flip() + theme_bw()
zp1 <- zp1 + ggtitle("Comparing several models")
print(zp1) # The trick to these is position_dodge().

#alternatively using coefplot
moddf <- data.frame(Coefficient = varnames, Value = mod1coefse[,1],
                    LowInner = mod1coefse[,1] - 1.64*mod1coefse[,2],
                    LowOuter = mod1coefse[,1] - 1.96*mod1coefse[,2],
                    HighInner = mod1coefse[,1] + 1.64*mod1coefse[,2],
                    HighOuter = mod1coefse[,1] + 1.96*mod1coefse[,2],
                    Model = "Baseline")
moddf2 <- data.frame(Coefficient = varnames, Value = mod2coefse[,1],
                    LowInner = mod2coefse[,1] - 1.64*mod2coefse[,2],
                    LowOuter = mod2coefse[,1] - 1.96*mod2coefse[,2],
                    HighInner = mod2coefse[,1] + 1.64*mod2coefse[,2],
                    HighOuter = mod2coefse[,1] + 1.96*mod2coefse[,2],
                    Model = "Model 2")
moddf3<- data.frame(Coefficient = varnames, Value = mod3coefse[,1],
                    LowInner = mod3coefse[,1] - 1.64*mod3coefse[,2],
                    LowOuter = mod3coefse[,1] - 1.96*mod3coefse[,2],
                    HighInner = mod3coefse[,1] + 1.64*mod3coefse[,2],
                    HighOuter = mod3coefse[,1] + 1.96*mod3coefse[,2],
                    Model = "Model 3")
coefplot(moddf)  #it suitably plot age with NA