orderedBeta
#163
Replies: 1 comment
-
Yep, probit link at the moment. |
Beta Was this translation helpful? Give feedback.
0 replies
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
-
Yep, probit link at the moment. |
Beta Was this translation helpful? Give feedback.
-
Hello,
I'm trying to match up the results of a glmmTMB and a gllvm for a simulated ordered beta GLM with only one response variable.
I'm simulating some data below.
#' load required packages
library(ggplot2)
library(glmmTMB)
library(DHARMa)
library(ordbetareg)
library(gllvm)
#' The code below carries out several tasks:
#' - Set the random seed.
#' - Define a covariate consisting of 1,000 values.
#' - Standardises this covariate to prevent numerical estimation problems.
#' - Calculate pi_i.
#' - Set the threshold parameters k_1 and k_2.
#' - Set the beta variance parameter theta.
#' - Simulate ordered beta-distributed data using the rordbeta function.
set.seed(123)
X <- seq(0, 4, length = 1000)
X <- (X - mean(X)) / sd(X)
Pi <- exp(-3 + 2 * X) / (1 + exp(-3 + 2 * X))
k1 <- -2
k2 <- 0.25
theta <- 25
Yordb <- rordbeta(n = 1000, mu = Pi,
phi = theta, cutpoints = c(k1, k2))
If you plot this data, then it looks like this:
#' Plot the simulated data.
plot(x = X,
y = Yordb,
xlab = "Covariate",
ylab = "Simulated ordered beta distributed data")
lines(x = X, y = Pi, lwd = 3)
I then apply an ordered beta GLM using glmmTMB:
#* Subsection 2.5: Apply glmmTMB----
#' We will fit an ordered beta GLM to this data using the glmmTMB function.
MyData <- data.frame(Yordb = Yordb, X = X)
M2 <- glmmTMB(Yordb ~ X,
family = ordbeta(link = "logit"),
data = MyData)
This is the output.
Family: ordbeta ( logit )
Formula: Yordb ~ X
Data: MyData
466.5 491.0 -228.2 456.5 995
Dispersion parameter for ordbeta family (): 21.2
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.91508 0.06570 -44.37 <2e-16 ***
X 1.99512 0.05573 35.80 <2e-16 ***
#' The estimated values for the intercept, slope, and theta are
#' relatively close to the ones I used to simulate the data.
#' The estimates for k_1 and k_2 can be located in the last
#' two elements of the M2$fit$parfull output.
ks <- tail(M2$fit$parfull,2)
ks
psi psi
-1.9680935 0.3048462
#' So..that is all good.
#' Now I run the same model in gllvm:
#' Prepare the data
YData <- matrix(MyData[,"Yordb"])
MyDataX <- data.frame(X = MyData[,"X"])
Testb <- gllvm(y = YData,
X = MyDataX,
formula = ~ X,
num.lv = 0,
link = "logit",
#method = "VA",
#seed = 12345,
#control.start = list(n.init = 5, jitter.var = 0.1),
family = "orderedBeta")
#' gllvm output:
summary(Testb)
Call:
gllvm(y = YData, X = MyDataX, formula = ~X, family = "orderedBeta",
num.lv = 0, method = "VA", link = "logit")
Family: orderedBeta
AIC: 479.0084 AICc: 479.0688 BIC: 503.5472 LL: -234.5 df: 5
Informed LVs: 0
Constrained LVs: 0
Unconstrained LVs: 0
Formula: ~X
LV formula: ~1
Coefficients predictors:
Estimate Std. Error z value Pr(>|z|)
X: 1.08903 0.02856 38.13 <2e-16 **
That is a rather different slope. Same for the intercept:
Testb$params$b
Col1
-1.629175
My question is: What am I doing wrong? Is the gllvm actually doing a logit link? The gllvm code has this:
if (family %in% c("beta", "betaH", "orderedBeta")) {
out$link <- link
}
If I use link = "donaldduck" in the gllvm call, it gives me the same results as above. And if I use a probit link in the glmmTMB, then the glmmTMB and gllvm results are rather similar.
Hence, I get the gut feeling that gllvm can only use a probit link for this distribution, at the moment?
And my second question is as follows.
I think that the expected values of an ordered beta GLM are as follows:
#' The expression for the expected values of a variable y_i following
#' an ordered beta distribution, conditional on the regression parameters,
#' is given below
#' E(y_i | Intercept, beta, theta, k_1, k_2) =
#' 0 * alpha + (1 - alpha - gamma) * Pi + gamma
glmmTMB gives Pi when you type fitted(M2)....and I think that gllvm with predict(Testb,type = "response") does the same.
I think it should be like this:
#' This necessitates first extracting all estimated parameters: the intercept
#' and slope (referred to as Betas), k_1 (k1.est), k_2 (k2.est), and theta.
Betas <- fixef(M2)$cond
k1.est <- ks[1]
k2.est <- ks[2]
theta <- summary(M2)$sigma
theta
#' Next, we compute the pi_i values (although fitted(M2) could have been
#' used for this purpose), as well as alpha_i and gamma_i.
MyLogit <- function(x) {exp(x) / (1 + exp(x))}
X <- model.matrix(M2)
Pi <- MyLogit(X %% Betas)
alpha <- 1 - MyLogit(X %% Betas - k1.est)
gamma <- MyLogit(X %*% Betas - k2.est)
#' We now have all three components required to calculate the expected
#' values of the ordered beta.
MyData$Exp.OrdBeta <- (1 - alpha - gamma) * Pi + gamma
head(MyData)
I can email the full code if that helps? Though I am anticipating that the answer is 'probit' link.
Kind regards,
Alain
Beta Was this translation helpful? Give feedback.
All reactions