Nothing
list.coef<-function(model,rounded=3,alpha=.05){
if(is(model,"glmerMod") | is(model,"lmerMod")){
out<-matrix(0,nrow=length(nlme::fixef(model)),ncol=10)
}else{out<-matrix(0,nrow=length(coef(model)),ncol=10)}
if(is(model,"multinom")){out[,1]<-t(coef(model))
}else if(is(model,"lme")){
out[,1]<- nlme::fixef(model)
}else if(is(model,"glmerMod") | is(model,"lmerMod")){
out[,1]<- nlme::fixef(model)}else{out[,1]<-coef(model)}
if(is(model,"lme")){vcov1 <- vcov(model)}else if(is(model,"glmerMod") | is(model,"lmerMod")){vcov1 <- vcov(model)}else{vcov1<-vcov(model)[1:length(coef(model)),1:length(coef(model))]}
out[,2]<-sqrt(diag(vcov1))
out[,3]<-out[,1]/out[,2]
out[,4] <- out[,1]-qnorm(1-(alpha/2),lower.tail=TRUE)*out[,2]
out[,5] <- out[,1]+qnorm(1-(alpha/2),lower.tail=TRUE)*out[,2]
out[,6]<-dnorm(out[,3])
if(is(model,"multinom")){ out[,7]<-exp(t(coef(model)))}else if(is(model,"lme")){out[,7]<-exp(nlme::fixef(model))}else if(is(model,"glmerMod") | is(model,"lmerMod")){out[,7]<-exp(nlme::fixef(model))}else{out[,7]<-exp(coef(model))}
out[,8]<-exp(out[,4])
out[,9]<-exp(out[,5])
out[,10]<-100*(exp(out[,1])-1)
out<-round(out,rounded)
colnames(out)<-c("b","SE","z","ll","ul","p.val","exp.b","ll.exp.b","ul.exp.b","percent")
out <- data.frame(out,CI=paste(100*(1-alpha),"%"))
if(is(model,"multinom")){
cn<-colnames(coef(model))
rn<-rownames(coef(model))
names <- paste(rn[1],cn)
for(i in 2:length(rn)){
namesi <- paste(rn[i],cn)
names <- c(names,namesi)}
out <- data.frame(variables=names,out)
}else if(is(model,"glmerMod") | is(model,"lmerMod")){out<-data.frame(variables=names(nlme::fixef(model)),out)}else{
out <- data.frame(variables=names(coef(model)),out)}
return(out)}
margins.des<-function (mod, ivs, excl = "nonE",data) {
if(is(mod,"lm") | is(mod,"glm") | is(mod,"polr") | is(mod,"multinom") | is(mod,"vglm") |
is(mod,"negbin") | is(mod,"zeroinfl") | is(mod,"hurdle") | is(mod,"glmerMod") |
is(mod,"clmm") | is(mod,"lme") | is(mod,"lmerModLmerTest") | is(mod,"lmerMod")){
if(is(mod,"nnet")){
c1<-as.character(mod$call)
m.polr<-MASS::polr(c1[2],data=data)
mod<-m.polr}
if(is(mod,"lme")){
X5 <- mod$data # pull out all of the data
var.names <- names(mod$fixDF$terms) # get the names of the variables that were used
var.names <- var.names[-1] # remove the intercept
X2 <- X5[,match(var.names,colnames(X5))] # from the full data, keep only the variables in the model
X.mod <- X2 # rename
if (excl == "nonE") {}else {
X.mod <- X2[, -match(excl, colnames(X2))]
var.names <- var.names[match(colnames(X.mod), var.names)]} # If you want to remove any factor-coded variables
X.mod <- X.mod[, -match(names(ivs), colnames(X.mod))] # Exclude "ivs" because those are not set to their means
var.names <- var.names[-match(names(ivs), var.names)] # Exclude "ivs" because those are not set to their means
if (is(X.mod,"numeric")) {
controls <- mean(X.mod)
names(controls) <- var.names
}else {
controls <- apply(X.mod, 2, FUN = "mean")} # Set covariates to their means
}else if(is(mod,"glmerMod")){
dv <- as.character(formula(mod))[[2]]
var.names <- names(nlme::fixef(mod)) #names of variables in the model
var.names <- var.names[-1] #remove intercept
# Seems like I need to add a variable to the function that defines the DV (so I can adjust for missing data)
data <- data[, match(c(dv,var.names), colnames(data))]
X.mod <- na.omit(data) # remove missing responses
X.mod <- X.mod[,-1] # Now, remove the dv from the predictors
if (excl == "nonE") {
}else {
X.mod <- X.mod[, -match(excl, colnames(X.mod))]
var.names <- var.names[match(colnames(X.mod), var.names)]
}
X.mod <- X.mod[, -match(names(ivs), colnames(X.mod))] # remove variables in the ivs statement
var.names <- var.names[-match(names(ivs), var.names)]
if (is(X.mod,"numeric")) {
controls <- mean(X.mod)
names(controls) <- var.names
}else {
controls <- apply(X.mod, 2, FUN = "mean")
}
}else if(is(mod,"lmerModLmerTest") | is(mod,"lmerMod")){
dv <- as.character(formula(mod))[[2]]
var.names <- names(nlme::fixef(mod)) #names of variables in the model
var.names <- var.names[-1] #remove intercept
# Seems like I need to add a variable to the function that defines the DV (so I can adjust for missing data)
data <- data[, match(c(dv,var.names), colnames(data))]
X.mod <- na.omit(data) # remove missing responses
X.mod <- X.mod[,-1] # Now, remove the dv from the predictors
if (excl == "nonE") {
}else {
X.mod <- X.mod[, -match(excl, colnames(X.mod))]
var.names <- var.names[match(colnames(X.mod), var.names)]
}
X.mod <- X.mod[, -match(names(ivs), colnames(X.mod))] # remove variables in the ivs statement
var.names <- var.names[-match(names(ivs), var.names)]
if (is(X.mod,"numeric")) {
controls <- mean(X.mod)
names(controls) <- var.names
}else {
controls <- apply(X.mod, 2, FUN = "mean")
}
}else if(is(mod,"clmm")){
X <- mod$model
X <- X[,-c(1,ncol(X))]
controls <- apply(X, 2, FUN = "mean")
controls <- controls[-match(names(ivs), colnames(X))]
}
else{
X.mod <- mod$model[, -1]
var.names <- colnames(X.mod)
if (excl == "nonE") {
}
else {
X.mod <- X.mod[, -match(excl, colnames(X.mod))]
var.names <- var.names[match(colnames(X.mod), var.names)]
}
X.mod <- X.mod[, -match(names(ivs), colnames(X.mod))]
var.names <- var.names[-match(names(ivs), var.names)]
if (is(X.mod,"numeric")) {
controls <- mean(X.mod)
names(controls) <- var.names
}
else {
controls <- apply(X.mod, 2, FUN = "mean")
}}
design <- ivs
design <- cbind(design, t(controls))
return(design)} else {message("Model type not supported.")}
}
margins.dat <- function (mod, des, alpha = 0.05, rounded = 3, cumulate = "no",
pscl.data = data, num.sample = 1000, prop.sample = 0.9, seed = 1234) {
pr <- NULL; prob <- NULL; SE <- NULL; fitted <- NULL ; se <- NULL
if(is(mod,"lm") | is(mod,"glm") | is(mod,"polr") | is(mod,"multinom") |
is(mod,"negbin") | is(mod,"zeroinfl") | is(mod,"hurdle") | is(mod,"glmerMod") |
is(mod,"clmm") | is(mod,"lme")| is(mod,"lmerModLmerTest") | is(mod,"lmerMod")){
if (cumulate == "no") {
if (sum(class(mod) == "lm") > 0) { # if I have to change this line, I'll need to differentiate various forms of glms
probs <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",
weights = "proportional", at = as.list(des[1,
]))))
if (nrow(des > 1)) {
for (i in 2:nrow(des)) {
probsi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",
weights = "proportional", at = as.list(des[i,
]))))
probs <- rbind(probs, probsi)
}
}
probs <- probs[, -c(1, 4)]
probs[, 3] <- probs[, 1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 2]
probs[, 4] <- probs[, 1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 2]
colnames(probs) <- c("fitted", "SE", "ll", "ul")
marginsdat <- data.frame(des, probs)
marginsdat <- round(marginsdat, rounded)
}
if (is(mod,"polr")) {
probs <- suppressMessages(data.frame(emmeans::emmeans(mod, specs=colnames(mod$model)[1], at = as.list(des[1,
]), weights = "proportional", mode = "prob")))
if (nrow(des) > 1) {
for (i in 2:nrow(des)) {
probsi <- suppressMessages(data.frame(emmeans::emmeans(mod, specs=colnames(mod$model)[1], at = as.list(des[i,
]), weights = "proportional", mode = "prob")))
probs <- rbind(probs, probsi)
}
}
probs[, 5] <- probs[, 2] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 3]
probs[, 6] <- probs[, 2] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 3]
colnames(probs)[5:6] <- c("ll", "ul")
probs <- probs[, -4]
probs[, 2] <- round(probs[, 2], rounded)
probs[, 3] <- round(probs[, 3], rounded)
probs[, 4] <- round(probs[, 4], rounded)
probs[, 5] <- round(probs[, 5], rounded)
des <- des[rep(1:nrow(des), each = length(table(mod$model[,1]))),
]
des <- round(des, rounded)
marginsdat <- data.frame(des, probs)
}
if (is(mod,"multinom")) {
cl <- as.character(mod$call)[2]
dv<-strsplit(cl," ")[[1]][1]
probs <- suppressMessages(data.frame(emmeans::emmeans(mod, specs=dv, at = as.list(des[1,
]), weights = "proportional", mode = "prob")))
if (nrow(des) > 1) {
for (i in 2:nrow(des)) {
probsi <- suppressMessages(data.frame(emmeans::emmeans(mod, specs=dv, at = as.list(des[i,
]), weights = "proportional", mode = "prob")))
probs <- rbind(probs, probsi)
}
}
probs[, 5] <- probs[, 2] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 3]
probs[, 6] <- probs[, 2] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 3]
colnames(probs)[5:6] <- c("ll", "ul")
probs <- probs[, -4]
probs[, 2] <- round(probs[, 2], rounded)
probs[, 3] <- round(probs[, 3], rounded)
probs[, 4] <- round(probs[, 4], rounded)
probs[, 5] <- round(probs[, 5], rounded)
des <- des[rep(1:nrow(des), each = length(mod$lev)),
]
des <- round(des, rounded)
marginsdat <- data.frame(des, probs)
}
if (is(mod,"zeroinfl")) {
p1 <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, at = as.list(des[1,
]), mode = "count", data = pscl.data)))
if (nrow(des) > 1) {
for (i in 2:nrow(des)) {
pi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, at = as.list(des[i,
]), mode = "count", data = pscl.data)))
p1 <- rbind(p1, pi)
}
}
p1 <- p1[, c(2, 3)]
colnames(p1)[1] <- "fitted"
marginsdat <- cbind(des, p1)
marginsdat <- dplyr::mutate(marginsdat, ll = fitted - qnorm(1 -
(alpha/2), lower.tail = TRUE) * SE, ul = fitted +
qnorm(1 - (alpha/2), lower.tail = TRUE) * SE)
marginsdat <- round(marginsdat, rounded)
}
if (is(mod,"hurdle")) {
p1 <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, at = as.list(des[1,
]), mode = "response"))[2:3])
if (nrow(des > 1)) {
for (i in 2:nrow(des)) {
pi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, at = as.list(des[i,
]), mode = "response"))[2:3])
p1 <- rbind(p1, pi)
}
}
ll <- p1$emmean - qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
ul <- p1$emmean + qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
marginsdat <- data.frame(round(des, rounded), fitted = round(p1$emmean,
rounded), se = round(p1$SE, rounded), ll = round(ll,
rounded), ul = round(ul, rounded))
}
if (is(mod,"lme")) {
p1 <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",
weights = "proportional", at = as.list(des[1,])))[2:3])
if (nrow(des > 1)) {
for (i in 2:nrow(des)) {
pi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",
weights = "proportional", at = as.list(des[i,])))[2:3])
p1 <- rbind(p1, pi)
}
}
ll <- p1$emmean - qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
ul <- p1$emmean + qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
marginsdat <- data.frame(round(des, rounded), fitted = round(p1$emmean,
rounded), se = round(p1$SE, rounded), ll = round(ll,
rounded), ul = round(ul, rounded))
}
if (is(mod,"glmerMod")) {
p1 <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",weights = "proportional", at = as.list(des[1,]))))
if (nrow(des > 1)) {
for (i in 2:nrow(des)) {
pi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",weights = "proportional", at = as.list(des[i,]))))
p1 <- rbind(p1, pi)
}
}
ll <- p1[,2] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
ul <- p1[,2] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
marginsdat <- data.frame(round(des, rounded), fitted = round(p1[,2],rounded), se = round(p1[,3], rounded), ll = round(ll,rounded), ul = round(ul, rounded))
}
if (is(mod,"lmerModLmerTest") | is(mod,"lmerMod")){
p1 <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",weights = "proportional", at = as.list(des[1,]))))
if (nrow(des > 1)) {
for (i in 2:nrow(des)) {
pi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~1, type = "response",weights = "proportional", at = as.list(des[i,]))))
p1 <- rbind(p1, pi)
}
}
ll <- p1[,2] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
ul <- p1[,2] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
p1$SE
marginsdat <- data.frame(round(des, rounded), fitted = round(p1[,2],rounded), se = round(p1[,3], rounded), ll = round(ll,rounded), ul = round(ul, rounded))
}
if (is(mod,"clmm")) {
out <- suppressMessages(data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(des[1,]))))
out <- out[,1:3]
colnames(out) <- c("dv.level","fitted","se")
out[,2:3]<- round(out[,2:3],rounded)
out <- suppressWarnings(suppressMessages(data.frame(des[1,],out)))
if(nrow(des) >1){for(i in 2:nrow(des)){
out2 <- suppressMessages(data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(des[i,]))))
out2 <- out2[,1:3]
colnames(out2) <- c("dv.level","fitted","se")
out2[,2:3]<- round(out2[,2:3],rounded)
out2 <- suppressWarnings(suppressMessages(data.frame(des[i,],out2)))
out <- rbind(out,out2)}}
marginsdat <- data.frame(dplyr::mutate(out,
ll=(fitted - qnorm(1 - (alpha/2), lower.tail = TRUE) *se),
ul=(fitted + qnorm(1 - (alpha/2), lower.tail = TRUE) *se)))
}
}
else {
if (is(mod,"polr")) {
probs <- suppressMessages(data.frame(emmeans::emmeans(mod, ~cut, at = as.list(des[1,
]), weights = "proportional", mode = "cum.prob")))
if (nrow(des) > 1) {
for (i in 2:nrow(des)) {
probsi <- suppressMessages(data.frame(emmeans::emmeans(mod, ~cut, at = as.list(des[i,
]), weights = "proportional", mode = "cum.prob")))
probs <- rbind(probs, probsi)
}
}
probs[, 5] <- probs[, 2] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 3]
probs[, 6] <- probs[, 2] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
probs[, 3]
colnames(probs)[5:6] <- c("ll", "ul")
probs <- probs[, -4]
probs[, 2] <- round(probs[, 2], rounded)
probs[, 3] <- round(probs[, 3], rounded)
probs[, 4] <- round(probs[, 4], rounded)
probs[, 5] <- round(probs[, 5], rounded)
des <- des[rep(1:nrow(des), each = length(table(mod$model[,
1])) - 1), ]
des <- round(des, rounded)
marginsdat <- data.frame(des, probs)
}
if (is(mod,"clmm")) {
out <- suppressMessages(data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(des[1,]))))
if(nrow(des) >1){for(i in 2:nrow(des)){
out2 <- suppressMessages(data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(des[i,]))))
out <- rbind(out,out2)}}
ll <- out[,2] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out$SE
ul <- out[,2] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out$SE
marginsdat <- data.frame(round(des, rounded),cut=out[,1], fitted = round(out[,2],rounded), se = round(out[,3], rounded), ll = round(ll,rounded), ul = round(ul, rounded))
}
}
colnames(marginsdat)[match("SE", colnames(marginsdat))] <- "se"
return(marginsdat)
} else {message("Model type not supported.")}}
# Take a Poisson model object of the count process.
count.fit<-function(m1,y.range,rounded=3,use.color="yes"){
se <- NULL; a <- NULL; b <- NULL
if(is(m1,"glm") & m1$family[1]$family=="poisson" ){
outcome<-m1$model[,1]
form <- as.character(m1$formula)
form.2<-as.formula(paste(form[2],form[1],form[3],"|",form[3]))
m1<-glm(m1$formula ,family="poisson",data=m1$data )
m2<- suppressMessages(MASS::glm.nb(m1$formula ,data=m1$data ))
m3 <- suppressMessages(pscl::zeroinfl(form.2,data=m1$data))
m4 <- pscl::zeroinfl(form.2,data=m1$data, dist="negbin")
bics<-c(BIC(m1),BIC(m2),BIC(m3),BIC(m4))
aics<-c(AIC(m1),AIC(m2),AIC(m3),AIC(m4))
out<-matrix(0,nrow=2*(ncol(m1$model)),ncol=16)
out[,1]<-c(coef(m1),rep(0,length(coef(m1))))
out[,5]<-c(coef(m2),rep(0,length(coef(m2))))
out[,9]<-c(coef(m3))
out[,13]<-c(coef(m4))
out[,2]<-c(sqrt(diag(vcov((m1)))),rep(0,length(coef(m1))))
out[,6]<-c(sqrt(diag(vcov((m2)))),rep(0,length(coef(m2))))
out[,10]<-c(sqrt(diag(vcov((m3)))))
out[,14]<-c(sqrt(diag(vcov((m4)))))
out[,3]<-out[,1]/out[,2]
out[,7]<-out[,5]/out[,6]
out[,11]<-out[,9]/out[,10]
out[,15]<-out[,13]/out[,14]
out[,4]<-dnorm(out[,3])
out[,8]<-dnorm(out[,7])
out[,12]<-dnorm(out[,11])
out[,16]<-dnorm(out[,15])
out<-round(out,rounded)
rownames(out)<-names(coef(m3))
colnames(out)<-c("Pcoef","Pse","Pz","Ppval","NBcoef","NBse","NBz","NBpval","ZIPcoef","ZIPse","ZIPz","ZIPpval","ZNBcoef","ZNBse","ZNBz","ZNBpval")
# Observed probabilities
out2<-outcome[which(outcome<(max(y.range)+1))]
obs<-0:max(y.range)
for(i in 1:length(y.range)){
obs[i]<-sum((out2==y.range[i]))/length(out2)}
poi<-round(dpois(y.range,lambda=mean(predict(m1,newdata=,type="response"))),rounded)
nb<-round(dnbinom(y.range,m2$theta,mu=mean(predict(m2,newdata=,type="response"))),rounded)
zpoi<-round(apply(pscl::predprob(m3),2,mean)[1:length(y.range)],rounded)
znb<-round(apply(pscl::predprob(m4),2,mean)[1:length(y.range)],rounded)
poi<-obs-poi
nb<-obs-nb
zpoi<-obs-zpoi
znb<-obs-znb
type<-c(rep("Pois",length(y.range)),rep("NB",length(y.range)),rep("ZIP",length(y.range)),rep("ZNB",length(y.range)))
dat<-data.frame(a=c(poi,nb,zpoi,znb),b=rep(y.range,4),type=type)
if(use.color=="yes"){
suppressMessages(suppressWarnings(
out.plot<- ggplot2::ggplot(data = dat, mapping = ggplot2::aes(x=b, y=a),group=type,linetype=type) +
ggplot2::geom_hline(yintercept=0,color="gray80") +
ggplot2::theme_bw() +
ggplot2::theme(legend.position="bottom") +
ggplot2::labs(x="Outcome",y="Observed - Predicted",color="") +
ggplot2::scale_x_continuous(breaks=seq(min(y.range),max(y.range),round(length(y.range)/4))) +
ggplot2::geom_line(ggplot2::aes(linetype=type,color=type)) + ggplot2::labs(linetype="",color="") +
ggplot2::scale_shape_manual(name="",values=c(22,21,23,24)) +
ggplot2::theme(panel.background = ggplot2::element_blank(), axis.line = ggplot2::element_line(colour = "black"))))}else{
suppressMessages(suppressWarnings(out.plot<- ggplot2::ggplot(data = dat, mapping = ggplot2::aes(x=b, y=a),group=type,linetype=type) +
ggplot2::geom_hline(yintercept=0,color="gray80") +
ggplot2::theme_bw() +
ggplot2::theme(legend.position="bottom") +
ggplot2::xlab("Outcome") +
ggplot2::ylab("Observed - Predicted") +
ggplot2::scale_x_continuous(breaks=seq(min(y.range),max(y.range),round(length(y.range)/4))) +
ggplot2::geom_line(ggplot2::aes(linetype=type)) + ggplot2::labs(linetype="") +
ggplot2::scale_shape_manual(name="",values=c(22,21,23,24)) +
ggplot2::theme(panel.background = ggplot2::element_blank(), axis.line = ggplot2::element_line(colour = "black"))))
}
models <- out[,c(1,2,5,6,9,10,13,14)]
pdat <- data.frame(var=c("Intercept",names(m1$coef)[2:length(names(m1$coefficients))]),
coef=as.numeric(c(models[,1],models[,3],models[,5],models[,7])),
se=as.numeric(c(models[,2],models[,4],models[,6],models[,8])),
model=rep(c("Poisson","Neg Binom","ZIP","ZNB"),each=nrow(models)))
pdat <- dplyr::mutate(pdat,max=coef+1.96*se,min=coef-1.96*se)
pdat$var <- paste(rep(rep(c("c","z"),each=length(m1$coefficients)),4),pdat$var)
pdat1<-pdat[1:(nrow(pdat)/4),]
pdat1$var <- factor(pdat1$var,levels=paste(pdat1$var))
pdat1[((nrow(pdat1)/2)+1):nrow(pdat1),2:3]<-NA
p1<- suppressMessages(suppressWarnings(ggplot2::ggplot(data = pdat1,mapping = ggplot2::aes(y=var,x=coef,xmin=min,xmax=max)) +
ggplot2::geom_vline(xintercept=0,color="gray80") +
ggplot2::geom_pointrange() +
ggplot2::theme_classic() + ggplot2::labs(x="",y="",title="Poisson") + ggplot2::scale_x_continuous(limits=c(min(pdat$min),max(pdat$max)))))
pdat2<-pdat[((nrow(pdat)/4)+1):((nrow(pdat)/2)),]
pdat2$var <- factor(pdat2$var,levels=paste(pdat2$var))
pdat2[((nrow(pdat2)/2)+1):nrow(pdat2),2:3]<-NA
p2<- suppressMessages(suppressWarnings(ggplot2::ggplot(data = pdat2,mapping = ggplot2::aes(y=var,x=coef,xmin=min,xmax=max)) +
ggplot2::geom_vline(xintercept=0,color="gray80") +
ggplot2::geom_pointrange() +
ggplot2::theme_classic() + ggplot2::labs(x="",y="",title="Neg Binom") + ggplot2::scale_x_continuous(limits=c(min(pdat$min),max(pdat$max)))))
pdat3<-pdat[((nrow(pdat)/2)+1):(((nrow(pdat)/4))*3),]
pdat3$var <- factor(pdat3$var,levels=paste(pdat3$var))
p3<- suppressMessages(suppressWarnings(ggplot2::ggplot(data = pdat3,mapping = ggplot2::aes(y=var,x=coef,xmin=min,xmax=max)) +
ggplot2::geom_vline(xintercept=0,color="gray80") +
ggplot2::geom_pointrange() +
ggplot2::theme_classic() + ggplot2::labs(x="",y="",title="ZIP") + ggplot2::scale_x_continuous(limits=c(min(pdat$min),max(pdat$max)))))
pdat4<-pdat[(1+(((nrow(pdat)/4))*3)):nrow(pdat),]
pdat4$var <- factor(pdat4$var,levels=paste(pdat4$var))
p4<- suppressMessages(suppressWarnings(ggplot2::ggplot(data = pdat4,mapping = ggplot2::aes(y=var,x=coef,xmin=min,xmax=max)) +
ggplot2::geom_vline(xintercept=0,color="gray80") +
ggplot2::geom_pointrange() +
ggplot2::theme_classic() + ggplot2::labs(x="",y="",title="ZINB") + ggplot2::scale_x_continuous(limits=c(min(pdat$min),max(pdat$max)))))
coef.pic<- suppressMessages(suppressWarnings(ggpubr::ggarrange(p1,p2,p3,p4,nrow=2,ncol=2)))
information.criterion<-rbind(bics,aics)
rownames(information.criterion)<-c("BIC","AIC")
colnames(information.criterion)<-c("Poisson","Neg Binom","ZIP","ZNB")
outp<-list(ic=information.criterion,pic=out.plot,models=out,models.pic=coef.pic)} else{
outp <- message("Model type should be Poisson regression as estimated by the glm function.")
}
return(outp)}
first.diff.fitted <- function (mod, design.matrix, compare, alpha = 0.05, rounded = 3,
bootstrap = "no", num.sample = 1000, prop.sample = 0.9, data,
seed = 1234,cum.probs="no") {
m1 <- NULL; design <- NULL
if(is(mod,"lm") | is(mod,"glm") | is(mod,"multinom") | is(mod,"vglm") |
is(mod,"negbin") | is(mod,"zeroinfl") | is(mod,"hurdle") |
is(mod,"glmerMod") | is(mod,"clmm") | is(mod,"polr")){
if (bootstrap == "no" & is(mod,"lm") | bootstrap == "no" & is(mod,"glm") |
bootstrap == "no" & is(mod,"negbin") | bootstrap == "no" & is(mod,"lme")){
des1 <- design.matrix[compare[1:2], ]
f2 <- margins.dat(mod, design.matrix[compare[1:2], ])
if (is(mod,"lm")) {
f1 <- function(x) predict(x, type = "response", newdata = des1[1,
]) - predict(x, type = "response", newdata = des1[2,
])
}
else {
f1 <- function(x) predict(x, type = "response", newdata = des1[1,
]) - predict(x, type = "response", newdata = des1[2,
])
}
out <- marginaleffects::hypotheses(mod,f1)
out <- c(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[5] <- out[1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
out[6] <- out[1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
names(out) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- as.matrix(t(out))
out <- round(out, rounded)
out <- as.data.frame(out)
if (sum(class(mod) == "nnet") == 1) {
out <- cbind(out, levels(f2[nrow(design.matrix)/2,
ncol(design.matrix) + 1]))
colnames(out)[ncol(out)] <- "dv"
}
if (length(compare) > 2) {
lc <- length(compare)/2
coms <- seq(1, length(compare), 2)
coms <- coms[-1]
for (j in 1:(lc - 1)) {
compare2 <- compare[coms[j]:(coms[j] + 1)]
f2 <- margins.dat(mod, design.matrix[compare2,
])
if (sum(class(mod) == "lm") > 0) {
f1 <- function(x) predict(x, type = "response",
newdata = design.matrix[compare2[1], ]) -
predict(x, type = "response", newdata = design.matrix[compare2[2],
])
}
else {
f1 <- function(x) predict(x, type = "response",
newdata = design.matrix[compare2[1], ]) -
predict(x, type = "response", newdata = design.matrix[compare2[2],
])
}
out2 <- marginaleffects::hypotheses(mod,f1)
out2 <- c(out2$estimate,out2$std.error,out2$statistic,out2$p.value,out2$conf.low,out2$conf.high)
out2[5] <- out2[1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[2]
out2[6] <- out2[1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[2]
names(out2) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out2 <- as.matrix(t(out2))
out2 <- round(out2, rounded)
out2 <- as.data.frame(out2)
if (sum(class(mod) == "nnet") == 1) {
out2 <- cbind(out2, levels(f2[nrow(design.matrix)/2,
ncol(design.matrix) + 1]))
colnames(out2)[ncol(out2)] <- "dv"
}
out <- rbind(out, out2)
}
}
}else if (bootstrap == "no" & is(mod,"multinom")){
des1 <- design.matrix[compare[1:2], ]
f1 <- function(x) predict(x, type = "probs", newdata = des1[1,]) - predict(x, type = "probs", newdata = des1[2,])
out <- marginaleffects::hypotheses(mod,f1)
out <- data.frame(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[,5] <- out[,1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
out[,6] <- out[,1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
names(out) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- round(out, rounded)
if (length(compare) > 2) {
lc <- length(compare)/2
coms <- seq(1, length(compare), 2)
coms <- coms[-1]
for (j in 1:(lc - 1)) {
compare2 <- compare[coms[j]:(coms[j] + 1)]
f1 <- function(x) predict(x, type = "probs",
newdata = design.matrix[compare2[1], ]) -
predict(x, type = "probs", newdata = design.matrix[compare2[2],
])
out2 <- marginaleffects::hypotheses(mod,f1)
out2 <- data.frame(out2$estimate,out2$std.error,out2$statistic,out2$p.value,out2$conf.low,out2$conf.high)
out2[,5] <- out2[,1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[,2]
out2[,6] <- out2[,1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[,2]
names(out2) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out2 <- round(out2, rounded)
out <- rbind(out, out2)
}
}
out <- data.frame(dv=names(predict(mod,type = "probs", newdata = des1[1,])),out)
names(out) <- c("Outcome Level","First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
} else if (bootstrap == "no" & is(mod,"vglm") | bootstrap == "no" & is(mod,"zeroinfl") | bootstrap == "no" & is(mod,"hurdle")){
out<-("Model not suppported with Delta method. Use bootstrapping.")
} else if (bootstrap == "no" & is(mod,"glmerMod")){
f1 <- function(x) predict(x, type = "response", newdata = design.matrix[compare[1],],re.form=NA) - predict(x, type = "response", newdata = design.matrix[compare[2],],re.form=NA)
out <- marginaleffects::hypotheses(mod,f1)
out <- c(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[5] <- out[1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
out[6] <- out[1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
names(out) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- as.matrix(t(out))
out <- round(out, rounded)
out <- as.data.frame(out)
if (length(compare) > 2) {
lc <- length(compare)/2
coms <- seq(1, length(compare), 2)
coms <- coms[-1]
for (j in 1:(lc - 1)) {
compare2 <- compare[coms[j]:(coms[j] + 1)]
f1 <- function(x) predict(x, type = "response", newdata = design.matrix[compare2[1],],re.form=NA) - predict(x, type = "response", newdata = design.matrix[compare2[2],],re.form=NA)
out2 <- marginaleffects::hypotheses(mod,f1)
out2 <- c(out2$estimate,out2$std.error,out2$statistic,out2$p.value,out2$conf.low,out2$conf.high)
out2[5] <- out2[1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[2]
out2[6] <- out2[1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[2]
names(out2) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out2 <- as.matrix(t(out2))
out2 <- round(out2, rounded)
out2 <- as.data.frame(out2)
out <- rbind(out, out2)}
}
}else if (bootstrap == "no" & is(mod,"polr") & cum.probs=="no"){
des1 <- design.matrix[compare[1:2], ]
f1 <- function(x) predict(x, type = "probs", newdata = des1[1,
]) - predict(x, type = "probs", newdata = des1[2,
])
out <- marginaleffects::hypotheses(mod,f1)
out <- data.frame(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[,5] <- out[,1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
out[,6] <- out[,1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
names(out) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- round(out, rounded)
if (length(compare) > 2) {
lc <- length(compare)/2
coms <- seq(1, length(compare), 2)
coms <- coms[-1]
for (j in 1:(lc - 1)) {
compare2 <- compare[coms[j]:(coms[j] + 1)]
f1 <- function(x) predict(x, type = "probs", newdata = design.matrix[compare2[1],]) - predict(x, type = "probs", newdata = design.matrix[compare2[2],])
out2 <- marginaleffects::hypotheses(mod,f1)
out2 <- data.frame(out2$estimate,out2$std.error,out2$statistic,out2$p.value,out2$conf.low,out2$conf.high)
out2[,5] <- out2[,1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[,2]
out2[,6] <- out2[,1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out2[,2]
names(out2) <- c("First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out2 <- round(out2, rounded)
out <- rbind(out, out2)
}
}
out <- data.frame(outcome.level=mod$lev,out)
names(out) <- c("Outcome Level","First Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
}else if (bootstrap == "no" & is(mod,"polr") & cum.probs=="yes"){
out <- message("Cumulative probabilities are not supported with parametric inference/delta method. Try bootstrapping or non-cumulative probabilities.")
}
else if (bootstrap == "no" & is(mod,"clmm")){out <- "Delta Method is not supported for mixed effects ordinal logistic regression. Use bootstrapping :("}
else if (bootstrap == "yes" & is(mod,"glmerMod")){out<-"Bootstrapping is not supported. Use the delta method." }
else if (bootstrap=="yes"){
# Start Bootstrap code
# Linear Regression
if(is(mod,"lm")){
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- lm(formula(mod),data=fd.model2)
obs.diffi<- predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
fd.dist[,1] <- sort(fd.dist[,1])
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(sd(fd.dist),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"glm")){
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- glm(formula(mod),data=fd.model2,family=family(mod))
obs.diffi<- predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
fd.dist[,1] <- sort(fd.dist[,1])
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(sd(fd.dist),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"polr")){
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- MASS::polr(formula(mod),data=fd.model2, Hess=TRUE)
obs.diffi<- predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"multinom")){
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "probs", newdata = des1[1,]) -
predict(mod, type = "probs", newdata = des1[2,])
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- data[sample(1:nrow(data),round(prop.sample*nrow(data),0),replace=TRUE),]
fd.modi <- nnet::multinom(formula(mod),data=fd.model2)
obs.diffi<- predict(fd.modi, type = "probs", newdata = des1[1,]) -
predict(fd.modi, type = "probs", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
names(out) <- c("First Difference","SD of the distribution","ll.boot","ul.boot")
} else if(is(mod,"vglm")){
out <- ("Partial proportional odds model not suppported in a general way. An example is provided in the following script: [first.diff.fitted.partial.prop.odds.e.g.r]")
} else if(is(mod,"negbin")){
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- MASS::glm.nb(formula(mod),data=fd.model2)
obs.diffi<- predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"zeroinfl")) {
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- pscl::zeroinfl(formula(mod),dist=mod$dist,data=fd.model2)
obs.diffi<- predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
names(out) <- c("First Difference","SD of the distribution","ll.boot","ul.boot")
} else if(is(mod,"hurdle")) {
des1 <- design.matrix[compare[1:2], ]
obs.diff<- predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- pscl::hurdle(formula(mod),dist=mod$dist$count,data=fd.model2)
obs.diffi<- predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
names(out) <- c("First Difference","SD of the distribution","ll.boot","ul.boot")
} else if(is(mod,"clmm") & cum.probs=="no"){
obs.diff<- suppressMessages(data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(design.matrix[compare[1],])))$prob -
data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(design.matrix[compare[2],])))$prob)
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- ordinal::clmm(formula(mod), data=fd.model2,na.action=na.omit)
obs.diffi<- suppressMessages(data.frame(emmeans::emmeans(fd.modi,~dv,mode="prob",at=as.list(design.matrix[compare[1],])))$prob -
data.frame(emmeans::emmeans(fd.modi,~dv,mode="prob",at=as.list(design.matrix[compare[2],])))$prob)
fd.dist[i,] <- obs.diffi}
fd.dist[,1] <- sort(fd.dist[,1])
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,FUN="sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"clmm") & cum.probs=="yes"){
obs.diff<- suppressMessages(data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(design.matrix[compare[1],])))$cumprob -
data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(design.matrix[compare[2],])))$cumprob)
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- ordinal::clmm(formula(mod), data=fd.model2,na.action=na.omit)
obs.diffi<- suppressMessages(data.frame(emmeans::emmeans(fd.modi,~cut,mode="cum.prob",at=as.list(design.matrix[compare[1],])))$cumprob -
data.frame(emmeans::emmeans(fd.modi,~cut,mode="cum.prob",at=as.list(design.matrix[compare[2],])))$cumprob)
fd.dist[i,] <- obs.diffi}
fd.dist <- apply(fd.dist,2,FUN="sort")
out <- data.frame(first.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,FUN="sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
}
}
return(out)
}else{message("Model type is not supported.")}}
second.diff.fitted <- function (mod, design.matrix, compare, alpha = 0.05, rounded = 3,
bootstrap = "no", num.sample = 1000, prop.sample = 0.9, data,
seed = 1234,cum.probs="no") {
if(is(mod,"lm") | is(mod,"glm") | is(mod,"polr") | is(mod,"multinom") | is(mod,"vglm") |
is(mod,"negbin") | is(mod,"zeroinfl") |
is(mod,"hurdle") | is(mod,"glmerMod") | is(mod,"clmm")){
if (bootstrap == "no" & is(mod,"lm") | bootstrap == "no" & is(mod,"glm") |
bootstrap == "no" & is(mod,"negbin")){
if (is(mod,"lm")) {
f1 <- function(x) (predict(x, type = "response",
newdata = design.matrix[compare[1], ]) - predict(x,
type = "response", newdata = design.matrix[compare[2],
])) - (predict(x, type = "response", newdata = design.matrix[compare[3],
]) - predict(x, type = "response", newdata = design.matrix[compare[4],
]))
}
else {
f1 <- function(x) (predict(x, type = "response", newdata = design.matrix[compare[1],
]) - predict(x, type = "response", newdata = design.matrix[compare[2],
])) - (predict(x, type = "response", newdata = design.matrix[compare[3],
]) - predict(x, type = "response", newdata = design.matrix[compare[4],
]))
}
out <- marginaleffects::hypotheses(mod,f1)
out <- c(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[5] <- out[1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
out[6] <- out[1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
names(out) <- c("Second Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- as.matrix(t(out))
out <- round(out, rounded)
out <- as.data.frame(out)
} else if (bootstrap == "no" & is(mod,"multinom")){
f1 <- function(x) (predict(x, type = "probs", newdata = design.matrix[compare[1],]) -
predict(x, type = "probs", newdata = design.matrix[compare[2],])) -
(predict(x, type = "probs", newdata = design.matrix[compare[3],]) -
predict(x, type = "probs", newdata = design.matrix[compare[4],]))
out <- marginaleffects::hypotheses(mod,f1)
out <- data.frame(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[,5] <- out[,1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
out[,6] <- out[,1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
names(out) <- c("Second Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- round(out, rounded)
out <- data.frame(dv=names(predict(mod,type = "probs", newdata = design.matrix[1,])),out)
names(out) <- c("Outcome Level","Second Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
} else if (bootstrap == "no" & is(mod,"glmerMod")){
f1 <- function(x) (predict(mod, type = "response",newdata = design.matrix[compare[1], ],re.form=NA) -
predict(mod,type = "response", newdata = design.matrix[compare[2],],re.form=NA)) -
(predict(mod, type = "response", newdata = design.matrix[compare[3],],re.form=NA) -
predict(x, type = "response", newdata = design.matrix[compare[4], ],re.form=NA))
out <- marginaleffects::hypotheses(mod,f1)
out <- c(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[5] <- out[1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
out[6] <- out[1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[2]
names(out) <- c("Second Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
out <- as.matrix(t(out))
out <- round(out, rounded)
out <- as.data.frame(out)
}else if (bootstrap == "no" & is(mod,"polr") & cum.probs=="no"){
des1 <- design.matrix[compare[1:4], ]
f1 <- function(x) (predict(x, type = "probs", newdata = des1[1,]) - predict(x, type = "probs", newdata = des1[2,])) -
(predict(x, type = "probs", newdata = des1[3,]) - predict(x, type = "probs", newdata = des1[4,]))
out <- marginaleffects::hypotheses(mod,f1)
out <- data.frame(out$estimate,out$std.error,out$statistic,out$p.value,out$conf.low,out$conf.high)
out[,5] <- out[,1] - qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
out[,6] <- out[,1] + qnorm(1 - (alpha/2), lower.tail = TRUE) *
out[,2]
out <- round(out, rounded)
out <- data.frame(outcome.level=mod$lev,out)
names(out) <- c("Outcome Level","Second Difference","Standard Error","Statistic","p-value" ,"ll", "ul")
}else if (bootstrap == "no" & is(mod,"polr") & cum.probs=="yes"){
out <- message("Cumulative probabilities are not supported with parametric inference/delta method. Try bootstrapping or non-cumulative probabilities.")
}else if (bootstrap == "no" & is(mod,"vglm") | bootstrap == "no" & is(mod,"zeroinfl") | bootstrap == "no" & is(mod,"hurdle") | bootstrap == "no" & is(mod,"clmm")){
out<-message("Model not suppported with Delta method. Use bootstrapping.")
} else if (bootstrap=="yes"){
# Start Bootstrap code
# Linear Regression
if(is(mod,"lm")){
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])) - (predict(mod, type = "response", newdata = des1[3,]) -
predict(mod, type = "response", newdata = des1[4,]))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- lm(formula(mod),data=fd.model2)
obs.diffi<- (predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])) - (predict(fd.modi, type = "response", newdata = des1[3,]) -
predict(fd.modi, type = "response", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
fd.dist[,1] <- sort(fd.dist[,1])
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(sd(fd.dist),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"glm")){
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])) - (predict(mod, type = "response", newdata = des1[3,]) -
predict(mod, type = "response", newdata = des1[4,]))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- glm(formula(mod),data=fd.model2,family=family(mod))
obs.diffi<- (predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])) - (predict(fd.modi, type = "response", newdata = des1[3,]) -
predict(fd.modi, type = "response", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
fd.dist[,1] <- sort(fd.dist[,1])
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(sd(fd.dist),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"polr")){
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])) - (predict(mod, type = "response", newdata = des1[3,]) -
predict(mod, type = "response", newdata = des1[4,]))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- MASS::polr(formula(mod),data=fd.model2, Hess=TRUE)
obs.diffi<- (predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])) - (predict(fd.modi, type = "response", newdata = des1[3,]) -
predict(fd.modi, type = "response", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"multinom")){
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "probs", newdata = des1[1,]) -
predict(mod, type = "probs", newdata = des1[2,])) - (predict(mod, type = "probs", newdata = des1[3,]) -
predict(mod, type = "probs", newdata = des1[4,]))
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- data[sample(1:nrow(data),round(prop.sample*nrow(data),0),replace=TRUE),]
fd.modi <- nnet::multinom(formula(mod),data=fd.model2)
obs.diffi<- (predict(fd.modi, type = "probs", newdata = des1[1,]) -
predict(fd.modi, type = "probs", newdata = des1[2,])) - (predict(fd.modi, type = "probs", newdata = des1[3,]) -
predict(fd.modi, type = "probs", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
names(out) <- c("Second Difference","SD of the distribution","ll.boot","ul.boot")
} else if(is(mod,"vglm")){
out <- ("Partial proportional odds model not suppported in a general way. An example is provided in the following script: [first.diff.fitted.partial.prop.odds.e.g.r]")
} else if(is(mod,"negbin")){
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])) - (predict(mod, type = "response", newdata = des1[3,]) -
predict(mod, type = "response", newdata = des1[4,]))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- MASS::glm.nb(formula(mod),data=fd.model2)
obs.diffi<- (predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])) - (predict(fd.modi, type = "response", newdata = des1[3,]) -
predict(fd.modi, type = "response", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"zeroinfl")) {
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])) - (predict(mod, type = "response", newdata = des1[3,]) -
predict(mod, type = "response", newdata = des1[4,]))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- pscl::zeroinfl(formula(mod),dist=mod$dist,data=fd.model2)
obs.diffi<- (predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])) - (predict(fd.modi, type = "response", newdata = des1[3,]) -
predict(fd.modi, type = "response", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
names(out) <- c("Second Difference","SD of the distribution","ll.boot","ul.boot")
} else if(is(mod,"hurdle")) {
des1 <- design.matrix[compare[1:4], ]
obs.diff<- (predict(mod, type = "response", newdata = des1[1,]) -
predict(mod, type = "response", newdata = des1[2,])) - (predict(mod, type = "response", newdata = des1[3,]) -
predict(mod, type = "response", newdata = des1[4,]))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- pscl::hurdle(formula(mod),dist=mod$dist$count,data=fd.model2)
obs.diffi<- (predict(fd.modi, type = "response", newdata = des1[1,]) -
predict(fd.modi, type = "response", newdata = des1[2,])) - (predict(fd.modi, type = "response", newdata = des1[3,]) -
predict(fd.modi, type = "response", newdata = des1[4,]))
fd.dist[i,] <- obs.diffi}
for(i in 1:ncol(fd.dist)){fd.dist[,i] <- sort(fd.dist[,i])}
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(apply(fd.dist,2,"sd"),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
names(out) <- c("Second Difference","SD of the distribution","ll.boot","ul.boot")
}else if(is(mod,"clmm") & cum.probs=="no"){
obs.diff<- suppressMessages((data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(design.matrix[compare[1],])))$prob -
data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(design.matrix[compare[2],])))$prob) -
(data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(design.matrix[compare[3],])))$prob -
data.frame(emmeans::emmeans(mod,~dv,mode="prob",at=as.list(design.matrix[compare[4],])))$prob))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- ordinal::clmm(formula(mod), data=fd.model2,na.action=na.omit)
obs.diffi<- suppressMessages((data.frame(emmeans::emmeans(fd.modi,~dv,mode="prob",at=as.list(design.matrix[compare[1],])))$prob -
data.frame(emmeans::emmeans(fd.modi,~dv,mode="prob",at=as.list(design.matrix[compare[2],])))$prob) -
(data.frame(emmeans::emmeans(fd.modi,~dv,mode="prob",at=as.list(design.matrix[compare[3],])))$prob -
data.frame(emmeans::emmeans(fd.modi,~dv,mode="prob",at=as.list(design.matrix[compare[4],])))$prob))
fd.dist[i,] <- obs.diffi}
fd.dist <- apply(fd.dist,2,FUN="sort")
out <- data.frame(second.diff=round(obs.diff,rounded),sd.boot.dist=round(sd(fd.dist),rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
} else if(is(mod,"clmm") & cum.probs=="yes"){
obs.diff<- suppressMessages((data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(design.matrix[compare[1],])))$cumprob -
data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(design.matrix[compare[2],])))$cumprob) -
(data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(design.matrix[compare[3],])))$cumprob -
data.frame(emmeans::emmeans(mod,~cut,mode="cum.prob",at=as.list(design.matrix[compare[4],])))$cumprob))
fd.model<-mod$model
fd.dist <-matrix(NA,nrow=num.sample,ncol=length(obs.diff))
for(i in 1:num.sample){
set.seed(seed + i); fd.model2 <- fd.model[sample(1:nrow(fd.model),round(prop.sample*nrow(fd.model),0),replace=TRUE),]
fd.modi <- ordinal::clmm(formula(mod), data=fd.model2,na.action=na.omit)
obs.diffi<- suppressMessages((data.frame(emmeans::emmeans(fd.modi,~cut,mode="cum.prob",at=as.list(design.matrix[compare[1],])))$cumprob -
data.frame(emmeans::emmeans(fd.modi,~cut,mode="cum.prob",at=as.list(design.matrix[compare[2],])))$cumprob) -
(data.frame(emmeans::emmeans(fd.modi,~cut,mode="cum.prob",at=as.list(design.matrix[compare[3],])))$cumprob -
data.frame(emmeans::emmeans(fd.modi,~cut,mode="cum.prob",at=as.list(design.matrix[compare[4],])))$cumprob))
fd.dist[i,] <- obs.diffi}
fd.dist <- apply(fd.dist,2,FUN="sort")
out <- data.frame(second.diff=round(obs.diff,rounded),ll.boot=round(fd.dist[nrow(fd.dist)*(alpha/2),],rounded),ul.boot=round(fd.dist[nrow(fd.dist)*(1-(alpha/2)),],rounded))
}else if(is(mod,"glmerMod")) {
out <- "Model type not supported with bootstrapping. Use Delta method."
}
}
return(out)
}else{message("Model type is not supported.")}
}
compare.margins <- function(margins,margins.ses,seed=1234,rounded=3,nsim=10000){
difference <- margins[1] - margins[2]
if(difference>0){
set.seed(seed); p.value<-sum((rnorm(nsim,mean=margins[1],sd=margins.ses[1]) - rnorm(nsim,mean=margins[2],sd=margins.ses[2])) < 0) /nsim
}else{
set.seed(seed); p.value<-sum((rnorm(nsim,mean=margins[1],sd=margins.ses[1]) - rnorm(nsim,mean=margins[2],sd=margins.ses[2])) > 0) /nsim
}
out <- c(difference,p.value)
names(out) <- c("Difference","p-value")
out <- as.data.frame(t(out))
return(out)
}
rubins.rule <- function(std.errors){
r.r.std.error<-sqrt(mean(std.errors^2) + var(std.errors) + var(std.errors)/length(std.errors))
return(r.r.std.error)}
diagn <- function(model){
if(is(model,"glm")){
pearsonres<-residuals.glm(model,type="pearson")
h<-data.frame(influence.measures(model)$infmat)$hat
stdpres<-pearsonres/(sqrt(1-h))
deltabeta<-(pearsonres^2*h)/(1-h)^2
deltabeta <- deltabeta/length(coef(model))
obs = 1:length(pearsonres)
devres<-residuals.glm(model,type="deviance")
out=data.frame(pearsonres,h,stdpres,deltabeta,obs,devres)
}
else if(is(model,"negbin")){
pearsonres<-residuals.glm(model,type="pearson")
h<-data.frame(influence.measures(model)$infmat)$hat
stdpres<-pearsonres/(sqrt(1-h))
deltabeta<-(pearsonres^2*h)/(1-h)^2
deltabeta <- deltabeta/length(coef(model))
obs = 1:length(pearsonres)
devres<-residuals.glm(model,type="deviance")
out=data.frame(pearsonres,h,stdpres,deltabeta,obs,devres)
}
else if(is(model,"polr")){
out=c("Ordinal model: you should binarize it for diagnostics.")
}
else if(is(model,"multinom")){
out=c("Multinomial model: you should binarize it for diagnostics.")
}
else if(is(model,"vglm")){
out=c("Partial Proportional Odds model: you should binarize it for diagnostics.")
}
else if(is(model,"zeroinfl")){
pearsonres<-residuals(model,type="pearson")
obs = 1:length(pearsonres)
out=data.frame(pearsonres,obs)
}
else if(is(model,"hurdle")){
pearsonres<-residuals(model,type="pearson")
obs = 1:length(pearsonres)
out=data.frame(pearsonres,obs)
}
else{out="Model type not supported."}
return(out)}
margins.dat.clogit <-function (mod, design.matrix, run.boot = "no", num.sample = 1000,
prop.sample = 0.9, alpha = 0.05, seed = 1234, rounded = 3) {
'(strata)' <- NULL
if(is(mod,"clogistic")){
coefs <- as.numeric(na.omit(coef(mod)))
des <- dplyr::mutate(design.matrix, lp = exp(as.matrix(design.matrix) %*%
coefs), probs = lp/sum(lp))
sims <- matrix(NA, ncol = nrow(design.matrix), nrow = num.sample)
vcovc <- vcov(mod)
sl <- which(is.na(coef(mod)))
if (length(sl) > 0) {
vcovc <- vcovc[-sl, -sl]
}
for (i in 1:num.sample) {
set.seed(seed + i)
lp <- exp(as.matrix(design.matrix) %*% MASS::mvrnorm(mu = coefs,
Sigma = vcovc))
sims[i, ] <- lp/sum(lp)
}
sims <- apply(sims, 2, FUN = "sort")
des <- dplyr::mutate(des, ll = sims[alpha/2 * nrow(sims), ], ul = sims[(1 -
alpha/2) * nrow(sims), ], se = apply(sims, 2, FUN = "sd"))
out <- round(des, rounded)
if (run.boot == "yes") {
boot.dist <- matrix(NA, nrow = num.sample, ncol = nrow(des))
for (i in 1:num.sample) {
set.seed(seed + i)
mod2 <- mod$model[sample(1:nrow(mod$model), round(prop.sample *
nrow(mod$model), 0), replace = TRUE), ]
m.1 <- Epi::clogistic(mod$formula, strata = `(strata)`,
data = mod2)
coefs2 <- as.numeric(na.omit(coef(m.1)))
design20 <- dplyr::mutate(des, lp = exp(as.matrix(design.matrix) %*%
coefs2), probs = lp/sum(lp))
boot.dist[i, ] <- design20[, ncol(design20)]
}
boot.dist[, 1] <- sort(boot.dist[, 1])
if (ncol(boot.dist) > 1) {
for (i in 2:ncol(boot.dist)) {
boot.dist[, i] <- sort(boot.dist[, i])
}
}
lower.limit <- boot.dist[nrow(boot.dist) * (alpha/2),
]
upper.limit <- boot.dist[nrow(boot.dist) * (1 - (alpha/2)),
]
des <- round(dplyr::mutate(des, ll.boot = lower.limit, ul.boot = upper.limit),
rounded)
out <- list(des = des, boot.dist = boot.dist)
}
return(out)} else if(is(mod,"mlogit")){
des <- dplyr::mutate(design.matrix, probs= predict(mod,newdata=design.matrix))
m2 <- mod
sims <- matrix(NA, ncol = nrow(design.matrix), nrow = num.sample)
vcovc <- vcov(mod)
for (i in 1:num.sample) {
set.seed(seed + i)
m2$coefficients <- MASS::mvrnorm(mu = coef(mod), Sigma = vcov(mod))
sims[i, ] <- predict(m2,newdata=design.matrix)}
sims <- apply(sims, 2, FUN = "sort")
des <- dplyr::mutate(des, ll = sims[alpha/2 * nrow(sims), ], ul = sims[(1 -
alpha/2) * nrow(sims), ], se = apply(sims, 2, FUN = "sd"))
out <- round(des, rounded)
return(out)
} else{message("Model type is not supported.")}}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.