1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
# 12/11/2017: S. Weisberg. This file contains all the Effect methods that call
# Effect.default. Excluded are Effect.lm, Effect.polr, and Effect.multinom,
# and for now Effect.svyglm.
# 06/08/2018: rewrote method for betareg, removing the 'link' argument from sources
# 11/28/2018: modified Effect.gls to ignore the weights argument by
# deleting it from sources$call.
# 11/30/2018: fixed bug in Effect.merMod() specifying fam$family explicitly.
# 7/5/2019: clm clm2 and clmm were not passing the estimated threshholds to polr
# 3/22/2020: added Effect.glmmPQL (from MASS package)
# 4/27/2020: require 'insight' package for find_formula and get_coefficients
# so formula and coefficients are generally not needed
# 2020-06-13: fix typo (omitted ') in an error message
# 2020-06-23: All the Effect.* methods previously in this file have been removed
# and replaced by effSources.* methods.
effSources <- function(mod){
UseMethod("effSources", mod)
}
effSources.default <- function(mod){NULL}
# lme, nlme package - default works
# gls, nlme package
effSources.gls <- function(mod){
cl <- mod$call
cl$weights <- NULL
list(call = cl)
}
# glmmPQL method 3/22/2020
effSources.glmmPQL <- function(mod) {list(family = mod$family)}
# lme4 -- handled via an Effect method to allow for KR argument
# effSources.merMod <- function(mod){NULL}
# rlmer in robustlmm package, not really needed
effSources.rlmerMod <- function(mod){NULL}
# clm in the ordinal package. clm is not supported by insight package
effSources.clm <- function(mod){
if (requireNamespace("MASS", quietly=TRUE)){
polr <- MASS::polr} else stop("MASS package is required")
polr.methods <- c("logistic", "probit", "loglog",
"cloglog", "cauchit")
method <- mod$link
if(method == "logit") method <- "logistic"
if(!(method %in% polr.methods))
stop("'link' must be a 'method' supported by polr; see help(polr)")
if(mod$threshold != "flexible")
stop("Effects only supports the 'flexible' threshold")
numTheta <- length(mod$Theta)
numBeta <- length(mod$beta)
or <- c( (numTheta+1):(numTheta + numBeta), 1:(numTheta))
list(
type = "polr",
coefficients = mod$beta,
zeta = mod$alpha,
method=method,
vcov = as.matrix(vcov(mod)[or, or]))
}
# clm2, this is supported by insight package
effSources.clm2 <- function(mod){
if (requireNamespace("MASS", quietly=TRUE)){
polr <- MASS::polr}
polr.methods <- c("logistic", "probit", "loglog",
"cloglog", "cauchit")
method <- mod$link
if(!(method %in% polr.methods))
stop("'link' must be a 'method' supported by polr; see help(polr)")
if(is.null(mod$Hessian)){
message("\nRe-fitting to get Hessian\n")
mod <- update(mod, Hess=TRUE)}
if(mod$threshold != "flexible")
stop("Effects only supports the flexible threshold")
numTheta <- length(mod$Theta)
numBeta <- length(mod$beta)
or <- c( (numTheta+1):(numTheta + numBeta), 1:(numTheta))
list(
type = "polr",
formula = mod$call$location,
coefficients = mod$beta,
zeta = mod$Theta,
method=method,
vcov = as.matrix(vcov(mod)[or, or]))
}
#clmm in ordinal package
effSources.clmm <- function(mod){
if (requireNamespace("MASS", quietly=TRUE)){
polr <- MASS::polr}
else stop("The MASS package must be installed")
polr.methods <- c("logistic", "probit", "loglog",
"cloglog", "cauchit")
method <- mod$link
if(method == "logit") method <- "logistic"
if(!(method %in% polr.methods))
stop("'link' must be a 'method' supported by polr; see help(polr)")
if(is.null(mod$Hessian)){
message("\nRe-fitting to get Hessian\n")
mod <- update(mod, Hess=TRUE)}
if(mod$threshold != "flexible")
stop("Only threshold='flexible' is supported by effects")
numTheta <- length(mod$Theta)
numBeta <- length(mod$beta)
or <- c( (numTheta+1):(numTheta + numBeta), 1:(numTheta))
Vcov <- as.matrix(vcov(mod)[or, or])
list(
type = "polr",
formula = insight::find_formula(mod)$conditional,
coefficients = mod$beta,
zeta=mod$alpha,
method=method,
vcov = as.matrix(Vcov))
}
# betareg from the betareg package
effSources.betareg <- function(mod){
coef <- mod$coefficients$mean
vco <- vcov(mod)[1:length(coef), 1:length(coef)]
# betareg uses beta errors with mean link given in mod$link$mean.
# Construct a family based on the binomial() family
fam <- binomial(link=mod$link$mean)
# adjust the variance function to account for beta variance
fam$variance <- function(mu){
f0 <- function(mu, eta) (1-mu)*mu/(1+eta)
do.call("f0", list(mu, mod$coefficient$precision))}
# adjust initialize
fam$initialize <- expression({mustart <- y})
# collect arguments
args <- list(
call = mod$call,
formula = formula(mod),
family=fam,
coefficients = coef,
vcov = vco)
args
}
|
