# R/lnl.ordinal.R In pglm: Panel Generalized Linear Models

#### Defines functions lnl.ordinal

```lnl.ordinal <- function(param, y, X, id, model, link, rn, start.sigma = FALSE){

if (start.sigma) gradient <- hessian <- TRUE

# F, f, e la fonction de distribution de proba et ses derivees
# premiere et seconde
F <- pnorm
f <- dnorm
e <- function(x) - x * f(x)
}
F <- function(x) exp(x)/(1+exp(x))
f <- function(x) exp(x)/(1+exp(x))^2
e <- function(x) (exp(x) - exp(3*x))/(1+exp(x))^4
}
K <- ncol(X)
n <- length(unique(id))
J <- length(unique(y))
beta <- param[1L:K]
mu <- c(-100, 0, param[(K+1L):(K+J-2L)], 100)
bX <- as.numeric(crossprod(t(X), beta))
if (model == "pooling"){
z1 <- mu[y+1] - bX                 ; z2 <- mu[y] - bX
F1 <- F(z1)                        ; F2 <- F(z2)
Li <- F1 - F2
}
if (model == "random"){
sigma <- param[K+J-1L]
Pitr <- lapply(rn\$nodes, function(x)
F(mu[y+1] - bX - sigma * x) - F(mu[y] - bX - sigma * x))
Pir <- lapply(Pitr, function(x) tapply(x, id, prod))
Li <- Reduce("+", mapply("*", Pir, rn\$weights, SIMPLIFY = FALSE))/sqrt(pi)
}
lnL <- sum(log(Li))
W <- matrix(0, nrow(X), J)
# The first two cols are not estimated coefficients (-infty and 0)
# so remove them
W1 <- (col(W) == (y + 1))[, -c(1, 2)]  ; W2 <- (col(W) == y)[, -c(1, 2)]
if (model == "pooling"){
f1 <- f(mu[y+1] - bX)                ; f2 <- f(mu[y] - bX)
gmu <- (W1 * f1 - W2 * f2)
gb <- - (f1 - f2) * X
gradi <- cbind(gb, gmu) / Li
}
if (model == "random"){
gitr <- mapply(
function(w, x, p){
f1 <- f(mu[y+1] - bX - sigma * x) ; f2 <- f(mu[y] - bX - sigma * x)
F1 <- F(mu[y+1] - bX - sigma * x) ; F2 <- F(mu[y] - bX - sigma * x)
w * as.numeric(p[as.character(id)]) *
(f1 * cbind(- 1, W1, - x) -
f2 * cbind(- 1, W2, - x)) / (F1 - F2)
},
rn\$weights, rn\$nodes, Pir, SIMPLIFY = FALSE)
g <- Reduce("+", gitr)
gradi <- cbind(X * g[, 1], g[, -c(1)])/as.numeric(Li[as.character(id)])/sqrt(pi)
}

if (model == "pooling"){
e1 <- e(mu[y+1] - bX)              ; e2 <- e(mu[y] - bX)
M1 <- cbind(- X, W1)
M2 <- cbind(- X, W2)
H <- crossprod(M1 * e1 / Li, M1) - crossprod(M2 * e2 / Li, M2) - crossprod(gradi)
}
if (model == "random"){
H <- mapply(
function(w, x, p){
p <- p / Li
P <- p[as.character(id)]
p <- as.numeric(p)
P <- as.numeric(P)
z1 <- mu[y+1] - bX - sigma * x
z2 <- mu[y] - bX - sigma * x
f1 <- f(z1) ; f2 <- f(z2)
e1 <- e(z1) ; e2 <- e(z2)
F1 <- F(z1) ; F2 <- F(z2)
M1 <- cbind(- X, W1, - x) ; M2 <- cbind(- X, W2, - x)
git <- (f1 * M1 - f2 * M2) / (F1 - F2)
gi <- apply(git, 2, tapply, id, sum)
H1 <- crossprod(p * gi, gi)
H2 <- crossprod(P * (e1 * M1) / (F1 - F2), M1) -
crossprod(P * (e2 * M2) / (F1 - F2), M2) -
crossprod( git * sqrt(P))
(H1 + H2) * w},
rn\$weights, rn\$nodes, Pir, SIMPLIFY = FALSE)
H <- Reduce("+", H) / sqrt(pi) - crossprod(apply(gradi, 2, tapply, id, sum))
}
if (start.sigma){
gradit <- (- f1 + f2 ) / Li
hessit <- (e1 - e2) / Li - gradit^2
mu <- tapply(gradit, id, sum) / tapply(hessit, id, sum)
return(sqrt(2) * sd(mu))
}