###{{{ Objective function
##' @export
tobit_method.lvm <- "nlminb1"
##' @export
tobit_objective.lvm <- function(x,p,data,weights,indiv=FALSE,
algorithm=lava::lava.options()$tobitAlgorithm,
seed=lava::lava.options()$tobitseed,...) {
if (!is.null(seed)) {
if (!exists(".Random.seed")) stats::runif(1)
save.seed <- .Random.seed
set.seed(seed)
}
requireNamespace("mvtnorm")
zz <- lava::manifest(x)
d <- as.matrix(rbind(data)[,zz,drop=FALSE]);
colnames(d) <- zz
yy <- lava::endogenous(x)
yy.w <- intersect(yy,colnames(weights))
yy.idx <- match(yy.w,zz)
Status <- matrix(0,ncol=length(zz),nrow=nrow(d))
Status[,yy.idx]<- as.matrix(weights)[,yy.w,drop=FALSE]
patterns <- unique(Status,MARGIN=1)
cens.type <- apply(Status,1,
function(x) which(apply(patterns,1,function(y) identical(x,y))))
mp <- lava::modelVar(x,p,data=as.data.frame(d))
Sigma <- mp$C ## Model specific covariance matrix
xi <- mp$xi ## Model specific mean-vector
val <- c()
for (i in 1:nrow(patterns)) {
## Usual marginal likelihood for status==1
pat <- patterns[i,]
idx <- which(cens.type==i)
noncens.idx <- which(pat==0)
left.cens.idx <- which(pat==-1)
right.cens.idx <- which(pat==1)
cens.idx <- which(pat!=0)
cens.which.left <- match(left.cens.idx,cens.idx)
cens.which.right <- match(right.cens.idx,cens.idx)
noncens.y <- zz[noncens.idx]
left.cens.y <- zz[left.cens.idx]
right.cens.y <- zz[right.cens.idx] ##setdiff(zz,noncens.y)
y <- d[idx,,drop=FALSE]
val1 <- val0 <- 0;
if (length(noncens.y)>0) {
## p(y), using: int[p(y,y*)]dy* = p(y) int[p(y*|y)]dy*
val1 <- mvtnorm::dmvnorm(d[idx,noncens.y,drop=FALSE],
mean=xi[noncens.idx],
sigma=Sigma[noncens.idx,noncens.idx,drop=FALSE],
log=TRUE)
}
if (length(cens.idx)>0) {
L <- diag(length(cens.idx))
L[cbind(cens.which.left,cens.which.left)] <- (-1)
## L[cens.which.left,cens.which.left] <- (-1)
if (length(noncens.y)==0) {
## suppressMessages(browser())
##low <- L%*%t(d[idx,cens.idx,drop=FALSE])-as.numeric(L%*%xi)
##lower=as.numeric(L%*%d[ii,cens.idx]),
val0 <- sapply(idx,
function(ii)
log(mvtnorm::pmvnorm(lower=as.numeric(L%*%d[ii,cens.idx]),
mean=as.numeric(L%*%as.vector(xi)),
sigma=L%*%Sigma%*%L,algorithm=algorithm))
)
} else {
M <- mom.cens(x,p,data=y,cens.idx,conditional=TRUE,deriv=FALSE)
val0 <- c()
for (j in 1:length(idx)) {
val0 <- c(val0,log(mvtnorm::pmvnorm(lower=as.numeric(L%*%y[j,cens.idx]),mean=as.numeric(L%*%M$mu.censIobs[j,]),sigma=L%*%M$S.censIobs%*%L,algorithm=algorithm)) )
}
}
}
val <- c(val,-val1-val0)
}
if (!is.null(seed))
.Random.seed <<- save.seed
if (!indiv)
return(sum(val))
val
}
###}}} Objective
###{{{ Gradient & Hessian
##' @export
tobit_gradient.lvm <- function(x,p,data,weights,data2=NULL,indiv=FALSE,
algorithm=lava::lava.options()$tobitAlgorithm,
seed=lava::lava.options()$tobitseed,...) {
if (!is.null(seed)) {
if (!exists(".Random.seed")) stats::runif(1)
save.seed <- .Random.seed
set.seed(seed)
}
requireNamespace("mvtnorm")
zz <- lava::manifest(x)
d <- as.matrix(data[,zz,drop=FALSE]); colnames(d) <- zz
yy <- lava::endogenous(x)
yy.w <- intersect(yy,colnames(weights))
yy.idx <- match(yy.w,zz)
Status <- matrix(0,ncol=length(zz),nrow=nrow(d))
Status[,yy.idx]<- as.matrix(weights)[,yy.w,drop=FALSE]
if (is.vector(data2)) data2 <- cbind(data2)
W0 <- data2
if (!is.null(W0)) {
yy.w2 <- intersect(yy,colnames(data2))
yy.idx2 <- match(yy.w2,zz)
W0 <- matrix(1,ncol=length(zz),nrow=nrow(d))
if (length(yy.idx2)==0) {
if (ncol(data2)==ncol(W0)) W0 <- data2
else {
for (i in seq(ncol(W0))) W0[,i] <- data2[,1]
}
} else {
W0[,yy.idx2] <- data2[,yy.w2,drop=FALSE]
}
colnames(W0) <- zz
}
## yy <- endogenous(x)
## yy.idx <- match(yy,zz)
## Status <- matrix(0,ncol=length(zz),nrow=nrow(d))
## Status[,yy.idx] <- as.matrix(weights)[,yy,drop=FALSE]
patterns <- unique(Status,MARGIN=1)
cens.type <- apply(Status,1,
function(x) which(apply(patterns,1,function(y) identical(x,y))))
val <- 0
## score <- c()
score <- matrix(ncol=length(p),nrow=nrow(data))
for (i in 1:nrow(patterns)) {
## Usual marginal likelihood for status==1
pat <- patterns[i,]
idx <- which(cens.type==i)
noncens.idx <- which(pat==0)
left.cens.idx <- which(pat==-1)
right.cens.idx <- which(pat==1)
cens.idx <- which(pat!=0)
cens.which.left <- match(left.cens.idx,cens.idx)
cens.which.right <- match(right.cens.idx,cens.idx)
noncens.y <- zz[noncens.idx]
left.cens.y <- zz[left.cens.idx]
right.cens.y <- zz[right.cens.idx] ##setdiff(zz,noncens.y)
y <- d[idx,,drop=FALSE]
w <- W0[idx,,drop=FALSE]
dummy <- cens.score(x,p,data=y,cens.idx=cens.idx,cens.which.left=cens.which.left, algorithm=algorithm,weights=w)
score[idx,] <- dummy
## y.pat <- unique(y,MARGIN=1)
## system.time(
## y.type <- apply(y,1,
## function(x) which(apply(y.pat,1,function
## (z) identical(x,z))))
## )
## print(system.time(
## dummy <- cens.score(x,p,data=y.pat,cens.idx=cens.idx, cens.which.left=cens.which.left, algorithm=algorithm)
## ))
## score <- rbind(score,dummy[y.type,,drop=FALSE])
## score <- rbind(0)
## dummy <- cens.score(x,p,data=y,cens.idx=cens.idx, cens.which.left=cens.which.left)
## score <- rbind(score,dummy)
}
if (!is.null(seed))
.Random.seed <<- save.seed
if (indiv)
return(-score)
return(-colSums(score))
}
##' @export
tobit_hessian.lvm <- function(x,p,data,weights,...) {
S <- -tobit_gradient.lvm(x,p=p,data=data,weights=weights,indiv=TRUE,...)
J <- t(S)%*%S
attributes(J)$grad <- colSums(S)
return(J)
}
###}}} Gradient & hessian
###{{{ Log-likelihood
##' @export
tobit_logLik.lvm <- function(object,p,data,weights,...) {
res <- -tobit_objective.lvm(x=object,p=p,data=data,weights=weights,...)
args <- list(...)
args$object <- object; args$p <- p; args$data <- data; args$type <- "exo"; args$weights <- NULL
##xl <- lava:::gaussian_logLik.lvm(object,p=p,data=data,type="exo",weights=NULL,...)
res <- res - do.call(lava::gaussian_logLik.lvm,args)
return(res)
}
###}}} Log-likelihood
###{{{ score for fixed censoring pattern
cens.score <- function(x,p,data,cens.idx,cens.which.left,weights,...) {
obs.idx <- setdiff(1:NCOL(data),cens.idx)
n <- NROW(data)
## print(system.time(
M <- mom.cens(x,p,data=data,cens.idx=cens.idx,conditional=TRUE,deriv=TRUE)
##))
## Censored part:
if (length(cens.idx)>0) {
combcens <- 1*(length(cens.which.left)>0) +
-1*(length(cens.which.left)<length(cens.idx))
## 0: left and right, -1: left only, 1: right only
## mu <- M$mu.cens
## S <- M$S.cens
## dmu <- M$dmu.cens
## dS <- M$dS.cens
z <- matrix(data[,cens.idx],nrow=n)
## z <- data[,cens.idx,drop=FALSE]
DCDFs <- c()
S <- M$S.censIobs
dS <- M$dS.censIobs
if (combcens==0) {
L <- -diag(length(cens.idx))
L[cbind(cens.which.left,cens.which.left)] <- 1
z <- z%*%L
S <- L%*%S%*%L
dS <- (L%x%L)%*%dS
}
w <- NULL
for (i in 1:n) {
mu <- M$mu.censIobs[i,]
dmu <- matrix(M$dmu.censIobs[,,i],nrow=length(cens.idx))
if (combcens==0) {
mu <- L%*%mu
dmu <- L%*%dmu
}
zi <- z[i,]
if (!is.null(weights))
w <- diag(weights[i,cens.idx],nrow=length(cens.idx))
if (combcens==-1) {
zi <- -zi
mu <- -mu
dmu <- -dmu
}
yy <- c(1,1); yy[cens.which.left] <- 0
## print(yy)
## print(mu)
## print(system.time(
DCDF <- Dthetapmvnorm(zi,
## mu=L%*%M$mu.censIobs[i,],
mu=mu,
## S=L%*%M$S.censIobs%*%L,
S=S,
## dS=(L%x%L)%*%M$dS.censIobs,
dS=dS,
## dmu=L%*%matrix(M$dmu.censIobs[,,i],nrow=length(cens.idx)),
dmu=dmu,
weights=w,
...)
alpha <- attributes(DCDF)$cdf
DCDFs <- rbind(DCDFs, 1/alpha*DCDF)
}
## DCDF <- Dthetapmvnorm(z,mu=M$mu.cens,S=S,dmu=M$dmu.cens,dS=M$dS.cens)
## S0 <- 1/attributes(DCDF)$cdf*DCDF
S0 <- DCDFs
} else S0 <- 0
## Observed part:
## S1 <- matrix() ...
if (length(obs.idx)>0) {
y1 <- as.matrix(data[,obs.idx,drop=FALSE])
if (n==1) y. <- y1 else y. <- colMeans(y1)
mu <- M$mu.obs
S <- M$S.obs
iS <- lava::Inverse(S)
dS <- M$dS.obs
dmu <- M$dmu.obs
S1 <- c()
part0 <- -1/2*as.vector(iS)%*%dS
for (i in 1:n) {
z <- as.numeric(y1[i,])
u <- z-mu
if (!is.null(weights)) {
W <- diag(weights[i,obs.idx],nrow=length(obs.idx))
S1 <- rbind(S1,
as.numeric(crossprod(W%*%u,iS)%*%dmu -
1/2*as.vector((iS-iS%*%tcrossprod(u)%*%iS)%*%W)%*%dS))
} else {
S1 <- rbind(S1,
as.numeric(part0 + crossprod(as.vector(u),iS)%*%dmu +
1/2*as.vector(iS%*%tcrossprod(as.vector(u))%*%iS)%*%dS))
}
}
} else S1 <- 0
## cat(rep("*",50),sep="")
## print(S1)
## cat(rep("-",50),sep="")
## print(S0)
## cat(rep("#",50),sep="")
## print(dim(data))
return((S1+S0))
}
###}}} score for fixed censoring pattern
###{{{ Derivatives of normal CDF
## Calculates first and second order partial derivatives of normal CDF
Dpmvnorm <- function(Y,S,mu=rep(0,NROW(S)),std=FALSE,
seed=lava::lava.options()$tobitseed,
algorithm=lava::lava.options()$tobitAlgorithm, ...) {
k <- NROW(S)
if (!is.null(seed) & k>1) {
if (!exists(".Random.seed")) stats::runif(1)
save.seed <- .Random.seed
}
requireNamespace("mvtnorm")
if (!std) {
L <- diag(S)^0.5
Li <- diag(1/L,NROW(S))
L <- diag(L,NROW(S))
S0 <- S
S <- Li%*%S%*%Li
## Y0 <- apply(Y,1,function(x)-mu)
Y <- (Li%*%(Y-mu))
}
## Y <- as.numeric(Y)
Y <- as.vector(Y)
if (k==1) {
D <- stats::dnorm(Y,sd=as.vector(S)^0.5)
H <- -Y*D
return(list(grad=D,hessian=H))
}
## Cond. var of Y[-j] given Y[j]
D <- numeric(k)
for (j in 1:k) {
## tcrossprod(S[-j,j,drop=FALSE])
Sj <- S[-j,-j,drop=FALSE] - tcrossprod(S[-j,j,drop=FALSE]) ##/S[j,j] S=correlation
muj <- Y[-j] - S[-j,j,drop=FALSE]*Y[j]
# set.seed(seed)
D[j] <- stats::dnorm(Y[j])*mvtnorm::pmvnorm(upper=as.vector(muj),sigma=Sj,algorithm=algorithm)
}
H <- matrix(0,k,k)
if (k<3) {
H[1,2] <- H[2,1] <- mvtnorm::dmvnorm(unlist(Y),sigma=S)
diag(H) <- -Y*D -H[1,2]*S[1,2]
##as.vector((H*S)%*%rep(1,k))
} else {
phis <- Phis <- H
for (i in 1:(k-1)) {
for (j in (i+1):k) {
Snij <- S[-c(i,j),c(i,j),drop=FALSE]
B <- Snij%*%lava::Inverse(S[c(i,j),c(i,j)])
Sij <- S[-c(i,j),-c(i,j),drop=FALSE] - B%*%t(Snij)
muij <- Y[-c(i,j)] - B%*%Y[c(i,j)]
Phis[i,j] <- Phis[j,i] <- mvtnorm::pmvnorm(upper=as.vector(muij),sigma=Sij,algorithm=algorithm)
phis[i,j] <- phis[j,i] <- mvtnorm::dmvnorm(Y[c(i,j)],sigma=S[c(i,j),c(i,j)])
}
}
H <- Phis*phis
diag(H) <- -Y*D - as.vector((S*H)%*%rep(1,k))
}
if (!std) {
if (!is.null(seed))
.Random.seed <<- save.seed
a <- mvtnorm::pmvnorm(upper=Y,mean=as.numeric(mu),sigma=S,algorithm=algorithm)
return(list(grad=Li%*%D, hessian=Li%*%H%*%Li, R=S, CDF=a, S=S0, mu=mu, L=L, Li=Li))
}
if (!is.null(seed))
.Random.seed <<- save.seed
return(list(grad=D,hessian=H))
}
## Calculates first and second order partial derivatives of normal CDF
## w.r.t. parameter-vector!
Dthetapmvnorm <- function(yy,mu,S,dmu,dS,
seed=lava::lava.options()$tobitseed,
algorithm=lava::lava.options()$tobitAlgorithm, weights,
...) {
if (!is.null(seed)) {
if (!exists(".Random.seed")) stats::runif(1)
save.seed <- .Random.seed
}
##yy <- as.matrix(yy)
## pp <- modelPar(x,p)
## M <- moments(x,p)
## mu <- M$xi
## S <- M$C
iS <- lava::Inverse(S)
## DCDF <- Dpmvnorm(y,S,mu)
L <- diag(S)^0.5
Li <- diag(1/L,NROW(S))
L <- diag(L,NROW(S))
R <- Li%*%S%*%Li
LR <- L%*%R ## = S%*%Li
if (!is.null(weights)) {
K1 <- -0.5*t(dS)%*%cbind(as.vector(iS%*%weights))
K2 <- 0.5*t(dS)%*%(iS%x%iS)
K3 <- t(dmu)%*%(iS%*%weights)
} else {
K1 <- -0.5*t(dS)%*%cbind(as.vector(iS))
K2 <- 0.5*t(dS)%*%(iS%x%iS)
K3 <- t(dmu)%*%(iS)
}
S0 <- function(y) {
z <- Li%*%(y-mu)
# set.seed(seed)
a <- mvtnorm::pmvnorm(upper=y,mean=as.numeric(mu),sigma=S,algorithm=algorithm)
DC <- Dpmvnorm(z,R,std=TRUE,algorithm=algorithm)
MM <- -LR%*%(DC$grad)
VV <- LR%*%(DC$hessian)%*%t(LR) + a*S
## message("DC")
## print(DC$hessian)
part1 <- K1*a
if (!is.null(weights)) {
VV <- VV%*%weights
}
part2 <- K2%*%as.vector(VV)
part3 <- K3%*%as.vector(MM)
res <- part1+part2+part3
return(c(a,res))
}
## DD <- t(apply(yy,1,S0))
DD <- matrix(S0(as.numeric(yy)),nrow=1)
a <- DD[,1]
res <- DD[,-1,drop=FALSE]
attributes(res)$cdf <- a
if (!is.null(seed))
.Random.seed <<- save.seed
return(res)
}
###}}} Derivatives of normal CDF
###{{{ Derivatives of Marginal and Conditional moments of a normal distribution
mom.cens <- function(x,p,cens.idx,data,deriv=TRUE,conditional=TRUE,right=TRUE,...) {
obs.idx <- setdiff(1:NCOL(data),cens.idx)
M <- lava::moments(x,p,data=as.data.frame(data))
if (deriv)
D <- stats::deriv(x,p=p,mom=M,meanpar=TRUE) ##,mu=colMeans(data))
if (length(cens.idx)<1) {
res <- list(S.obs=M$C, mu.obs=M$xi, S.cens=NULL, mu.cens=NULL,
S.censIobs=NULL, mu.censIobs=NULL)
if (deriv)
res <- c(res, list(dS.obs=D$dS, dmu.obs=D$dxi, dS.cens=NULL, dmu.cens=NULL,
dS.censIobs=NULL, dmu.censIobs=NULL
))
return(res)
}
if (length(obs.idx)<1) {
res <- list(S.obs=NULL, mu.obs=NULL, S.cens=M$C, mu.cens=M$xi,
S.censIobs=M$C, mu.censIobs=matrix(M$xi,ncol=length(M$xi),nrow=NROW(data),byrow=TRUE)
)
if (deriv)
res <- c(res, list(dS.obs=NULL, dmu.obs=NULL, dS.cens=D$dS, dmu.cens=D$dxi,
dS.censIobs=D$dS, dmu.censIobs=array(D$dxi,dim=c(dim(D$dxi),NROW(data)))
))
return(res)
}
S.cens <- M$C[cens.idx,cens.idx,drop=FALSE]
S.obs <- M$C[obs.idx,obs.idx,drop=FALSE]
S.censobs <- M$C[cens.idx,obs.idx,drop=FALSE]
S.obscens <- M$C[obs.idx,cens.idx,drop=FALSE]
mu.obs <- M$xi[obs.idx]
mu.cens <- M$xi[cens.idx]
iS.obs <- lava::Inverse(S.obs)
S01iS0 <- S.censobs%*%iS.obs
S.cond <- S.cens - S01iS0%*%S.obscens
y0.obs <- apply(data,1,function(z) z[obs.idx]-mu.obs)
Zy0 <- S01iS0%*%y0.obs
mu.cond <- t(rbind(apply(Zy0,2,function(z) mu.cens+z)))
res <- list(S.obs=S.obs, mu.obs=mu.obs, S.cens=S.cens, mu.cens=mu.cens)
if (conditional)
res <- c(res, list(S.censIobs=S.cond, mu.censIobs=mu.cond))
if (!deriv)
return(res)
M.idx <- matrix(1:nrow(M$C)^2,nrow(M$C))
Cens <- as.vector(M.idx[cens.idx,cens.idx])
Obs <- as.vector(M.idx[obs.idx,obs.idx])
CensObs <- as.vector(M.idx[cens.idx,obs.idx])
ObsCens <- as.vector(M.idx[obs.idx,cens.idx])
dS.cens <- D$dS[Cens,,drop=FALSE]
dS.obs <- D$dS[Obs,,drop=FALSE]
dS.censobs <- D$dS[CensObs,,drop=FALSE]
dS.obscens <- D$dS[ObsCens,,drop=FALSE]
res <- c(res, list(##dT.obs=D$dT[Obs,],
dS.obs=dS.obs, dmu.obs=D$dxi[obs.idx,,drop=FALSE],
dS.cens=dS.cens, dmu.cens=D$dxi[cens.idx,,drop=FALSE])
)
if(!conditional)
return(res)
K0 <- S.censobs%*%iS.obs
K1 <- K0%x%diag(length(cens.idx))
K2 <- K0%x%K0
### K3 <- diag(length(obs.idx))%x%K0
K3 <- diag(length(cens.idx))%x%K0
dSc <-
dS.cens -
K1%*%dS.censobs +
K2%*%dS.obs -
K3%*%dS.obscens
dxic.1 <- D$dxi[cens.idx,] - K0%*%D$dxi[obs.idx,]
dxic <- array(dxic.1, dim=c(nrow(dxic.1),ncol(dxic.1),nrow(data)))
for (k in 1:length(p)) {
dxic[,k,] <- dxic[,k,] +
(
matrix(dS.censobs[,k],nrow=length(cens.idx),ncol=length(obs.idx)) -
K0%*%matrix(dS.obs[,k],ncol=length(obs.idx))
)%*%iS.obs%*%y0.obs
}
## y <- as.matrix(d[1,obs.idx])
## y0 <- y-M$xi[obs.idx]
## Ky <- (y0)%*%solve(S.obs)
## K1 <- Ky%x%diag(length(cens.idx))
## K2 <- Ky%x%K0
## dxic.2 <- K1%*%dS.censobs - K2%*%dS.obs
## dxic <- dxic.1+dxic.2
res <- c(res, list(#D=invisible(D),
dS.censIobs=dSc, dmu.censIobs=dxic))
return(res)
}
###}}} Derivatives of Marginal and Conditional moments of a normal distribution
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.