Nothing
start.values.gllvm.TMB <- function(y, X = NULL, lv.X = NULL, TR=NULL, family,
offset= NULL, trial.size = 1, num.lv = 0, num.lv.c = 0, num.RR = 0, start.lvs = NULL,
seed = NULL,Power=NULL,starting.val="res",formula=NULL, lv.formula = NULL,
jitter.var=0,yXT=NULL, row.eff=FALSE, TMB=TRUE,
link = "probit", randomX = NULL, beta0com = FALSE, zeta.struc="species", maxit=4000,max.iter=4000, disp.group = NULL, randomB = FALSE, method="VA", Ntrials = 1) {
if(!is.null(seed)) set.seed(seed)
N <- n <- nrow(y); p = ncol(y); y = as.matrix(y)
num.T = 0; if(!is.null(TR)) num.T = dim(TR)[2]
num.X = 0; Xdesign = NULL
if(!is.null(X)){
if(!is.null(formula) & is.null(TR)){
Xdesign = model.matrix(formula, as.data.frame(X))
Xdesign = Xdesign[, !(colnames(Xdesign) %in% "(Intercept)"), drop=FALSE]
} else {
Xdesign = X
}
num.X = dim(Xdesign)[2]
}
Br <- sigmaB <- sigmaij <- NULL
mu = NULL
if(method=="LA")
method = "VA"
out = list()
sigma = 1
row.params <- rep(0, n);
if(starting.val %in% c("res","random") && !isFALSE(row.eff) || row.eff == "random"){
rmeany <- rowMeans(y/matrix(Ntrials, ncol = ncol(y), nrow=nrow(y), byrow = T), na.rm = TRUE)
if(family=="binomial"){
rmeany=(1e-3+0.99*rmeany)
if(row.eff %in% c("fixed",TRUE)) {
row.params = binomial(link = link)$linkfun(rmeany) - binomial(link = link)$linkfun(rmeany[1])
} else{
row.params = binomial(link = link)$linkfun(rmeany) - binomial(link = link)$linkfun(mean(rmeany))
}
} else if(family=="gaussian"){
rmeany = 1e-3+0.99*rmeany
if(row.eff %in% c("fixed",TRUE)) {
row.params = rmeany - rmeany[1]
} else{
row.params = rmeany - mean(rmeany)
}
} else {
if(row.eff %in% c("fixed",TRUE)) {
row.params = row.params <- log(rmeany)-log(rmeany[1])
} else{
row.params <- row.params <- log(rmeany)-log(mean(y, na.rm = TRUE))
}
}
if(any(abs(row.params)>1.5)) row.params[abs(row.params)>1.5] = 1.5 * sign(row.params[abs(row.params)>1.5])
sigma = sd(row.params)
}
if(!is.numeric(y))
stop("y must a numeric. If ordinal data, please convert to numeric with lowest level equal to 1. Thanks")
# if(family=="ZIP") family = "poisson"
# if(family=="ZINB") family="negative.binomial"
# if(family=="betaH") family = "beta"
# if(family=="orderedBeta") family = "beta"
# if(family=="betaH") family="beta"
if(!(family %in% c("poisson","negative.binomial","binomial","ordinal","tweedie", "gaussian", "gamma", "exponential", "beta", "betaH", "orderedBeta","ZIP","ZINB")))
stop("inputed family not allowed...sorry =(")
if((num.lv+num.lv.c) > 0) {
unique.ind = which(!duplicated(y))
if(is.null(start.lvs)) {
index = MASS::mvrnorm(N, rep(0, num.lv+num.lv.c),diag(num.lv+num.lv.c));
if(num.lv>0&num.lv.c==0)colnames(index) = paste("LV",1:num.lv, sep = "")
if(num.lv==0&num.lv.c>0)colnames(index) = paste("CLV",1:num.lv.c, sep = "")
if(num.lv>0&num.lv.c>0)colnames(index) = c(paste("CLV",1:num.lv.c, sep = ""),paste("LV",1:num.lv, sep = ""))
unique.index = as.matrix(index[unique.ind,])
}
if(!is.null(start.lvs)) {
if(num.lv.c>0){
index.lm = lm(start.lvs ~ 0+lv.X)
b.lv = coef(index.lm)
start.lvs = as.matrix(residuals.lm(index.lm))
}
index = as.matrix(start.lvs)
unique.index = as.matrix(index[unique.ind,])
}
}
if((num.lv+num.lv.c) == 0) { index <- NULL }
y = as.matrix(y)
if(family == "ordinal" && zeta.struc == "species") {
max.levels = apply(y,2,function(x) length(min(x):max(x)));
if(any(max.levels == 1) || all(max.levels == 2)) stop("Ordinal data requires all columns to have at least has two levels. If al columns only have two levels, please use family == binomial instead. Thanks")
}else if(family=="ordinal" && zeta.struc == "common"){
max.levels = length(min(y):max(y))
}
if(family == "orderedBeta") {
zetaOB = matrix(rep(0,p),rep(1,p),p,2)
}
if(is.null(rownames(y))) rownames(y) <- paste("row",1:N,sep="")
if(is.null(colnames(y))) colnames(y) <- paste("col",1:p,sep="")
options(warn = -1)
if(family!="ordinal") { ## Using logistic instead of prbit regession here for binomial, but whatever...
if(starting.val=="res" && is.null(start.lvs) ){# && num.lv>0
if(is.null(TR)){
if(family!="gaussian") {
if(!is.null(X)) fit.mva <- gllvm.TMB(y=y, X=X, formula = formula, lv.X = lv.X, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link =link, Power = Power, disp.group = disp.group, method=method, Ntrials = Ntrials)#mvabund::manyglm(y ~ X, family = family, K = trial.size)
if(is.null(X)) fit.mva <- gllvm.TMB(y=y, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link =link, Power = Power, disp.group = disp.group, method=method, Ntrials = Ntrials)#mvabund::manyglm(y ~ 1, family = family, K = trial.size)
coef = cbind(fit.mva$params$beta0,fit.mva$params$Xcoef)
fit.mva$phi = fit.mva$params$phi
if(family == "orderedBeta") {zetaOB = fit.mva$zeta = fit.mva$params$zeta}
if(family=="ZINB")fit.mva$ZINB.phi = fit.mva$params$ZINB.phi
if(family=="tweedie")Power = fit.mva$Power
resi = NULL
mu = cbind(rep(1,n),fit.mva$X.design)%*%t(cbind(fit.mva$params$beta0, fit.mva$params$Xcoef))
} else {
if(!is.null(X)) fit.mva <- mlm(y, X = Xdesign)
if(is.null(X)) fit.mva <- mlm(y)
mu = cbind(rep(1,nrow(y)),Xdesign) %*% fit.mva$coefficients
# resi = fit.mva$residuals; resi[is.infinite(resi)] <- 0; resi[is.nan(resi)] <- 0
coef = t(fit.mva$coef)
fit.mva$phi = sapply(1:length(unique(disp.group)),function(x)sd(fit.mva$residuals[,which(disp.group==x)]))[disp.group]
}
gamma=NULL
if((num.lv+num.lv.c+num.RR)>0){
lastart <- FAstart(eta=mu, family=family, y=y, num.lv = num.lv, num.lv.c = num.lv.c, num.RR = num.RR, phis=fit.mva$phi, lv.X = lv.X, zeta = fit.mva$zeta, link = link, maxit=maxit,max.iter=max.iter, Power = Power, disp.group = disp.group, randomB = randomB, method = method, Ntrials = Ntrials)
gamma<-lastart$gamma
index<-lastart$index
if(num.lv.c>0){
b.lv<-lastart$b.lv
if(num.RR>0){
b.lv <- cbind(lastart$b.lv, lastart$RRcoef)
}
}else if(num.RR>0){
b.lv <- lastart$b.lv
}
}
} else {
n1 <- colnames(X)
n2 <- colnames(TR)
form1 <-NULL
if(is.null(formula)){
form1 <- "~ 1"
for(m in 1:length(n2)){
for(l in 1:length(n1)){
ni <- paste(n1[l],n2[m],sep = "*")
form1 <- paste(form1,ni,sep = "+")
}
}
formula=form1
}
if(!TMB) fit.mva <- gllvm.VA(y, X = X, TR = TR, formula = formula(formula), family = family, num.lv = 0, Lambda.struc = "diagonal", trace = FALSE, plot = FALSE, sd.errors = FALSE, maxit = 1000, max.iter=200, seed = seed, n.init = 1, starting.val="zero", yXT = yXT)
if(TMB) {
fit.mva <- try(trait.TMB(y, X = X, TR = TR, formula = formula(formula), family = family, num.lv = 0, Lambda.struc = "diagonal", trace = FALSE, sd.errors = FALSE, maxit = 1000, max.iter=200, seed=seed,n.init=1,starting.val="zero",yXT = yXT, row.eff = row.eff, diag.iter = 0, optimizer = "nlminb", beta0com = beta0com, link = link, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials), silent = TRUE);
fit.mva$method = "VA"
if(!is.null(randomX) && !inherits(fit.mva, "try-error") && is.finite(fit.mva$logL)) {
fit.mva <- trait.TMB(y, X = X, TR = TR, formula=formula(formula), family = family, num.lv = 0, Lambda.struc = "diagonal", trace = FALSE, sd.errors = FALSE, maxit = 1000, max.iter=200, seed=seed,n.init=1,starting.val="zero",yXT = yXT, row.eff = row.eff, diag.iter = 0, optimizer = "nlminb", randomX = randomX, beta0com = beta0com, start.params = fit.mva, link = link, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials);
} else {fit.mva <- trait.TMB(y, X = X, TR = TR, formula=formula(formula), family = family, num.lv = 0, Lambda.struc = "diagonal", trace = FALSE, sd.errors = FALSE, maxit = 1000, max.iter=200, seed=seed,n.init=1,starting.val="zero",yXT = yXT, row.eff = row.eff, diag.iter = 0, optimizer = "nlminb", randomX = randomX, beta0com = beta0com, link = link, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials);}
fit.mva$coef=fit.mva$params
if(row.eff=="random") { # !!!!
sigma=c(max(fit.mva$params$sigma[1],sigma),fit.mva$params$sigma[-1])
fit.mva$params$row.params <- fit.mva$params$row.params/sd(fit.mva$params$row.params)*sigma[1]
}
if(family=="tweedie")Power = fit.mva$Power
out$fitstart <- list(A=fit.mva$A, Ab=fit.mva$Ab, TMBfnpar=fit.mva$TMBfn$par) #params = fit.mva$params,
}
if(!is.null(form1)){
if(!is.null(fit.mva$coef$row.params)) row.params=fit.mva$coef$row.params
env <- (fit.mva$coef$B)[n1]
trait <- (fit.mva$coef$B)[n2]
inter <- (fit.mva$coef$B)[!names(fit.mva$coef$B) %in% c(n1,n2)]
B <- c(env,trait,inter)
} else {
B<-fit.mva$coef$B
if(!is.null(fit.mva$coef$row.params)) row.params=fit.mva$coef$row.params
}
if(family == "orderedBeta") {zetaOB = fit.mva$zeta = fit.mva$params$zeta}
if(family=="ZINB") fit.mva$ZINB.phi = fit.mva$params$ZINB.phi
fit.mva$phi <- phi <- fit.mva$coef$phi
ds.res <- matrix(NA, n, p)
rownames(ds.res) <- rownames(y)
colnames(ds.res) <- colnames(y)
mu <- (matrix(fit.mva$X.design%*%fit.mva$coef$B,n,p)+ matrix(fit.mva$coef$beta0,n,p,byrow = TRUE))
if(row.eff %in% c(TRUE, "random", "fixed")) {mu <- mu + row.params }
if(!is.null(randomX)) {
Br <- fit.mva$params$Br
sigmaB <- fit.mva$params$sigmaB
if(ncol(fit.mva$Xrandom)>1) sigmaij <- fit.mva$TMBfn$par[names(fit.mva$TMBfn$par)=="sigmaij"]#fit.mva$params$sigmaB[lower.tri(fit.mva$params$sigmaB)]
mu <- mu + fit.mva$Xrandom%*%Br
}
gamma=NULL
if((num.lv+num.lv.c+num.RR)>0){
lastart <- FAstart(eta=mu, family=family, y=y, num.lv = num.lv, num.lv.c = num.lv.c, phis=fit.mva$phi, lv.X = lv.X, zeta = fit.mva$zeta, link = link, maxit=maxit,max.iter=max.iter, disp.group = disp.group, randomB = randomB, method = method, Ntrials = Ntrials)
gamma<-lastart$gamma
index<-lastart$index
if(num.lv.c>0)
{
b.lv<-lastart$b.lv
if(num.RR>0){
b.lv <- cbind(lastart$b.lv, lastart$RRcoef)
}
}else if(num.RR>0){
b.lv <- lastart$b.lv
}
}
}
if(is.null(TR)){params = cbind(coef,gamma)
} else { params = cbind((fit.mva$coef$beta0),gamma)}
# if(family == "orderedBeta") {
# zetaOB = fit.mva$params$zeta
# }
} else if(starting.val != "zero") {
if(family!="gaussian") {
if(is.null(TR)){
if(!is.null(X) & (num.lv+num.lv.c) > 0) {
if(is.null(formula)) {
if(is.null(colnames(index))) colnames(index) <- paste(index, 1:ncol(index), sep = "")
formula = formula(paste(paste(formula, collapse =""), paste(colnames(index), collapse = " + "), sep = "+"))
}
fit.mva <- gllvm.TMB(y=y, X=cbind(X,index), formula = formula, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, maxit = maxit, max.iter = max.iter, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials)
if(family=="tweedie")Power = fit.mva$Power
}
if(is.null(X) & (num.lv+num.lv.c) > 0) {
fit.mva <- gllvm.TMB(y=y, X=index, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, maxit = maxit, max.iter = max.iter, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials)
if(family=="tweedie")Power = fit.mva$Power
}
if(!is.null(X) & (num.lv+num.lv.c) == 0) {
fit.mva <- gllvm.TMB(y=y, X=X, formula = formula, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, maxit= maxit, max.iter= max.iter, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials)
if(family=="tweedie")Power = fit.mva$Power
}
if(is.null(X) & (num.lv+num.lv.c) == 0) {
fit.mva <- gllvm.TMB(y=y, X=NULL, family = family, num.lv = 0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, maxit = maxit, max.iter = max.iter, Power = Power, disp.group = disp.group, method = method, Ntrials = Ntrials)
if(family=="tweedie")Power = fit.mva$Power
}
} else {
if((num.lv+num.lv.c) > 0) fit.mva <- gllvm.TMB(y=y, X=index, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, maxit = maxit, max.iter = max.iter, disp.group = disp.group, method = method, Ntrials = Ntrials)
if((num.lv+num.lv.c) == 0) fit.mva <- gllvm.TMB(y=y, X=NULL, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, maxit = maxit, max.iter = max.iter, disp.group = disp.group, method = method, Ntrials = Ntrials)
env <- rep(0,num.X)
trait <- rep(0,num.T)
inter <- rep(0, num.T * num.X)
B <- c(env,trait,inter)
}
params <- cbind(fit.mva$params$beta0, fit.mva$params$Xcoef, fit.mva$params$theta)
fit.mva$phi <- fit.mva$params$phi
} else {
if(is.null(TR)){
if(!is.null(X) & (num.lv+num.lv.c) > 0) fit.mva <- mlm(y, X = Xdesign, index = index)
if(is.null(X) & (num.lv+num.lv.c) > 0) fit.mva <- mlm(y, index = index)
if(!is.null(X) & (num.lv+num.lv.c) == 0) fit.mva <- mlm(y, X = Xdesign)
if(is.null(X) & (num.lv+num.lv.c) == 0) fit.mva <- mlm(y)
} else {
if((num.lv+num.lv.c) > 0) fit.mva <- mlm(y, index = index)
if((num.lv+num.lv.c) == 0) fit.mva <- mlm(y)
env <- rep(0,num.X)
trait <- rep(0,num.T)
inter <- rep(0, num.T * num.X)
B <- c(env,trait,inter)
}
params <- t(fit.mva$coefficients)
fit.mva$phi <- sapply(1:length(unique(disp.group)),function(x)sd(fit.mva$residuals[,which(disp.group==x)]))[disp.group]
}
#params <- t(fit.mva$coef) #!!
} else {
fit.mva<-list()
fit.mva$phi <- rep(1, p)
if (family %in% c("betaH", "orderedBeta")) {
fit.mva$phi <- rep(5,p)
}
}
}
if(family == "negative.binomial") {
phi <- fit.mva$phi + 1e-5
} else if(family %in% c("gaussian", "gamma", "beta", "betaH", "orderedBeta","tweedie","ZIP","ZINB")) {
phi <- fit.mva$phi
} else { phi <- NULL }
if(family == "ZINB"){
ZINB.phi <- fit.mva$ZINB.phi + 1e-5
}
#
# if(family == "tweedie") {
# if(is.null(TR)){
# indeX <- cbind(index, X)
# }else{indeX <- cbind(index)}
# if((num.lv+num.lv.c) == 0) indeX <- X
#
# phi0<-NULL
# coefs0<-NULL
# for(j in 1:p){
# fitj <- try(glm( y[,j] ~ indeX, family=statmod::tweedie(var.power=power, link.power=0) ),silent = TRUE)
# if(is.null(X) && (num.lv+num.lv.c) == 0) fitj<-try(glm( y[,j] ~ 1, family=statmod::tweedie(var.power=power, link.power=0) ),silent = TRUE)
# if(!inherits(fitj, "try-error")){
# coefs0<-rbind(coefs0,fitj$coefficients)
# phi0<-c(phi0,summary(fitj)$dispersion)
# } else {coefs0<-rbind(coefs0,rnorm((num.lv+num.lv.c)+dim(X)[2]+1)); phi0<-c(phi0,runif(1,0.2,3));}
# }
#
# if(!is.null(TR)) B=rep(0,num.X+num.T+num.X*num.T)
# params <- as.matrix(coefs0)
# phi <- phi0
#
# if(starting.val%in%c("res") && (num.lv+num.lv.c+num.RR)>0){
# eta.mat <- matrix(params[,1],n,p,byrow=TRUE)
# mu <- eta.mat
# if((num.lv+num.lv.c+num.RR)>0){
# lastart <- FAstart(eta.mat, family=family, y=y, num.lv = num.lv, num.lv.c = num.lv.c, num.RR= num.RR, zeta = zeta, zeta.struc = zeta.struc, lv.X = lv.X, link = link, maxit=maxit,max.iter=max.iter, phis=phi, Power = power)
# gamma<-lastart$gamma
# index<-lastart$indexunt
# params <- cbind(params,gamma)
# if(num.lv.c>0)
# {
# b.lv<-lastart$b.lv
# if(num.RR>0){
# b.lv <- cbind(lastart$b.lv, lastart$RRcoef)
# }
# }else if(num.RR>0){
# b.lv <- lastart$b.lv
# }
# }
# }
#
#
# }
if(family == "ordinal") {
max.levels <- length(unique(c(y)))
params <- matrix(0,p,ncol(cbind(1,Xdesign))+(num.lv+num.lv.c+num.RR))
env <- rep(0,num.X)
trait <- rep(0,num.T)
inter <- rep(0, num.T * num.X)
B=c(env,trait,inter)
if(zeta.struc == "species"){
zeta <- matrix(0,p,max.levels - 1)
zeta[,1] <- 0 ## polr parameterizes as no intercepts and all cutoffs vary freely. Change this to free intercept and first cutoff to zero
}else{
cw.fit <- MASS::polr(factor(y) ~ 1, method = "probit")
zeta <- cw.fit$zeta
zeta[1] <- 0
}
if(starting.val=="res"){
if(!is.null(X)) fit.mva <- gllvm.TMB(y=y, X=X, formula = formula, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, zeta.struc = zeta.struc, maxit = maxit, max.iter = max.iter, method = method, Ntrials = Ntrials)
if(is.null(X)) fit.mva <- gllvm.TMB(y=y, family = family, num.lv=0, starting.val = "zero", sd.errors = FALSE, optimizer = "nlminb", link = link, zeta.struc = zeta.struc, maxit = maxit, max.iter = max.iter, method = method, Ntrials = Ntrials)
params[,1:ncol(cbind(fit.mva$params$beta0,fit.mva$params$Xcoef))] <- cbind(fit.mva$params$beta0,fit.mva$params$Xcoef)
zeta <- fit.mva$params$zeta
resi <- NULL
eta.mat <- cbind(rep(1,n),fit.mva$X.design)%*%t(cbind(fit.mva$params$beta0, fit.mva$params$Xcoef))
if(!is.finite(fit.mva$logL)|is.na(fit.mva$logL)){
stop("Could not calculate starting values for ordinal model. Change starting values or zeta.struc, further simplify your model, or center and scale your predictors.")
}
} else {
for(j in 1:p) {
y.fac <- factor(y[,j])
if(length(levels(y.fac)) > 2) {
if(starting.val%in%c("zero") || (num.lv+num.lv.c)==0){
if(is.null(X) ) try(cw.fit <- MASS::polr(y.fac ~ 1, method = "probit"),silent = TRUE)
if(!is.null(X) ) try(cw.fit <- MASS::polr(y.fac ~ Xdesign, method = "probit"),silent = TRUE)
} else {
if(is.null(X)) try(cw.fit <- MASS::polr(y.fac ~ index, method = "probit"),silent = TRUE)
if(!is.null(X)) try(cw.fit <- MASS::polr(y.fac ~ Xdesign+index, method = "probit"),silent = TRUE)
}
if(starting.val=="random"){
params[j,] <- c(cw.fit$zeta[1],-cw.fit$coefficients)
if(zeta.struc == "species"){
zeta[j,2:length(cw.fit$zeta)] <- cw.fit$zeta[-1]-cw.fit$zeta[1]
}
}else{
params[j,1:length( c(cw.fit$zeta[1],-cw.fit$coefficients))]<- c(cw.fit$zeta[1],-cw.fit$coefficients)
if(zeta.struc == "species"){
zeta[j,2:length(cw.fit$zeta)] <- cw.fit$zeta[-1]-cw.fit$zeta[1]
}
}
}
if(length(levels(y.fac)) == 2) {
if(starting.val%in%c("zero") || (num.lv+num.lv.c)==0){
if(is.null(X)) try(cw.fit <- glm(y.fac ~ 1, family = binomial(link = "probit")),silent = TRUE)
if(!is.null(X)) try(cw.fit <- glm(y.fac ~ Xdesign, family = binomial(link = "probit")),silent = TRUE)
} else { # || (!is.null(TR) && NCOL(TR)>0) & is.null(TR)
if(is.null(X)) try(cw.fit <- glm(y.fac ~ index, family = binomial(link = "probit")),silent = TRUE)
if(!is.null(X)) try(cw.fit <- glm(y.fac ~ Xdesign+index, family = binomial(link = "probit")),silent = TRUE)
}
params[j,1:length(cw.fit$coef)] <- cw.fit$coef
}
}
}
if(starting.val%in%c("res") && (num.lv+num.lv.c+num.RR)>0){
eta.mat <- matrix(params[,1],n,p,byrow=TRUE)
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + (Xdesign %*% matrix(params[,2:(1+num.X)],num.X,p))
mu <- eta.mat
if((num.lv+num.lv.c+num.RR)>0){
lastart <- FAstart(eta.mat, family=family, y=y, num.lv = num.lv, num.lv.c = num.lv.c, num.RR= num.RR, zeta = zeta, zeta.struc = zeta.struc, lv.X = lv.X, link = link, maxit=maxit,max.iter=max.iter, disp.group = disp.group, randomB = randomB, method = method, Ntrials = Ntrials)
gamma<-lastart$gamma
index<-lastart$index
params[,(ncol(cbind(1,X))+1):ncol(params)]=gamma
if(num.lv.c>0)
{
b.lv<-lastart$b.lv
if(num.RR>0){
b.lv <- cbind(lastart$b.lv, lastart$RRcoef)
}
}else if(num.RR>0){
b.lv <- lastart$b.lv
}
}
}
}
if(starting.val=="random"&(num.lv.c+num.RR)>0){
if(num.lv.c>0){
index.lm <- lm(index[,1:num.lv.c,drop=F]~0+lv.X)
b.lv<-coef(index.lm)
index[,1:num.lv.c] <- residuals.lm(index.lm)
}
if(num.RR>0)b.lv2 <- matrix(rnorm(ncol(lv.X)*num.RR),nrow=ncol(lv.X),ncol=num.RR)
if(num.lv.c>0&num.RR>0){b.lv <- cbind(b.lv, b.lv2)}else if(num.lv.c==0&num.RR>0)b.lv<-b.lv2
if(num.lv==0){
params <- cbind(params,matrix(rnorm(num.RR*p),ncol=num.RR,nrow=p))
}else{
params <- cbind(params[,1:(ncol(params)-num.lv)],matrix(rnorm(num.RR*p),ncol=num.RR,nrow=p),params[,(ncol(params)-num.lv+1):ncol(params)])
}
}
if((num.lv.c+num.lv)>0){
if((family!="ordinal" || (family=="ordinal" & starting.val=="res")) & starting.val!="zero"){
if(num.lv==0 &&p>2 & (num.lv.c+num.RR)>1 | (num.lv.c+num.RR) == 0 && p>2 & num.lv>1){
gamma<-as.matrix(params[,(ncol(params) - num.lv - num.lv.c - num.RR + 1):ncol(params)])
qr.gamma <- qr(t(gamma))
if(num.lv.c>0)b.lv <- b.lv%*%qr.Q(qr.gamma)
params[,(ncol(params) - num.lv - num.lv.c - num.RR + 1):ncol(params)]<-t(qr.R(qr.gamma))
if(num.lv.c>1)index<-(index%*%qr.Q(qr.gamma)[1:num.lv.c,1:num.lv.c])
}else if((num.lv.c+num.RR)>1 && num.lv>1 && p>2){
gamma<-as.matrix(params[,(ncol(params) - num.lv.c - num.lv - num.RR + 1):ncol(params)])
qr.gamma1 <- qr(t(gamma[,1:(num.lv.c+num.RR)]))
qr.gamma2 <- qr(t(gamma[,(num.lv.c+num.RR+1):ncol(gamma)]))
params[,(ncol(params) - num.lv.c - num.lv - num.RR + 1):(ncol(params)- num.lv)]<-t(qr.R(qr.gamma1))
params[,(ncol(params) - num.lv + 1):ncol(params)]<-t(qr.R(qr.gamma2))
b.lv <- b.lv%*%qr.Q(qr.gamma1)
if(num.lv.c>0)index[,1:num.lv.c]<-(index[,1:num.lv.c,drop=F]%*%qr.Q(qr.gamma1)[1:num.lv.c,1:num.lv.c])
index[,(num.lv.c+1):ncol(index)]<-(index[,(num.lv.c+1):ncol(index),drop=F]%*%qr.Q(qr.gamma2))
}
}
}
if(starting.val=="zero"){
params=matrix(0,p,1+num.X+(num.lv+num.lv.c+num.RR))
params[,1:(ncol(params) - (num.lv+num.lv.c+num.RR))] <- 0
env <- rep(0,num.X)
trait <- rep(0,num.T)
inter <- rep(0, num.T * num.X)
B=c(env,trait,inter)
if(num.lv > 0 && (num.lv.c+num.RR) == 0) {
gamma <- matrix(1,p,num.lv)
if(num.lv>1)gamma[upper.tri(gamma)]=0
params[,(ncol(params) - num.lv + 1):ncol(params)] <- gamma
index <- matrix(0,n,num.lv)
}else if(num.lv == 0 & (num.lv.c+num.RR) > 0){
gamma <- matrix(1,p,num.lv.c+num.RR)
if((num.lv.c+num.RR)>1)gamma[upper.tri(gamma)]=0
params[,(ncol(params) - num.lv.c - num.RR + 1):ncol(params)] <- gamma
index <- matrix(0,n,num.lv.c)
}else if(num.lv > 0 & num.lv.c>0){
gamma <- matrix(1,p,num.lv+num.lv.c+num.RR)
if((num.lv.c+num.RR)>1)gamma[,1:(num.lv.c+num.RR)][upper.tri(gamma[,1:(num.lv.c+num.RR)])]=0
if(num.lv>1)gamma[,((num.lv.c+num.RR)+1):ncol(gamma)][upper.tri(gamma[,((num.lv.c+num.RR)+1):ncol(gamma)])] <- 0
params[,(ncol(params) - (num.lv.c+num.RR) - num.lv + 1):ncol(params)] <- gamma
index <- matrix(0,n,num.lv.c+num.lv)
}
phi <- rep(1,p)
if((num.lv.c+num.RR)>0){
b.lv <- matrix(0.1,nrow=ncol(lv.X),ncol=(num.lv.c+num.RR))
}
sigma.lv <- rep(1,num.lv+num.lv.c)
}
if((num.lv+num.lv.c+num.RR) > 0 & starting.val!="zero") {
if((num.lv+num.lv.c)>0)index <- index+ MASS::mvrnorm(n, rep(0, num.lv+num.lv.c),diag(num.lv+num.lv.c)*jitter.var);
if(num.lv>0&(num.lv.c+num.RR)==0){
try({
gamma.new <- params[,(ncol(params) - num.lv + 1):(ncol(params)),drop=F];
sigma.lv <- diag(gamma.new)
sign <- sign(diag(gamma.new))
params[,(ncol(params) - num.lv + 1):(ncol(params))] <- t(t(gamma.new)/sigma.lv)
index <- t(t(index)*sig)}, silent = TRUE)
}else if(num.lv==0 & (num.lv.c+num.RR) >0){
try({
if(num.lv.c>0){
gamma.new <- params[,(ncol(params) - (num.lv.c+num.RR) + 1):(ncol(params)),drop=F];
sig <- sign(diag(gamma.new[,1:num.lv.c,drop=F]));
sigma.lv <- diag(gamma.new)
params[,(ncol(params) - (num.lv.c+num.RR) + 1):(ncol(params))] <- t(t(gamma.new)/sigma.lv)
index <- t(t(index)*sig)}}, silent = TRUE)
}else if(num.lv >0 & (num.lv.c+num.RR) >0){
try({
gamma.new2 <- params[,(ncol(params) - num.lv + 1):(ncol(params)),drop=F];
sigma.lv2 <- diag(gamma.new2)
sig2 <- sign(diag(gamma.new2));
params[,(ncol(params) - num.lv + 1):(ncol(params))] <- t(t(gamma.new2)/sigma.lv2)
sigma.lv <- NULL
if(num.lv.c>0){
gamma.new <- params[,(ncol(params) - num.lv.c - num.lv - num.RR + 1):(ncol(params)-num.lv),drop=F];
sig <- sign(diag(gamma.new));
sigma.lv <- diag(gamma.new)
params[,(ncol(params) - num.lv.c - num.lv -num.RR + 1):(ncol(params)-num.lv)] <- t(t(gamma.new)/sigma.lv)
index[,1:num.lv.c,drop=F] <- t(t(index[,1:num.lv.c,drop=F])*sig)
}
index[,(num.lv.c+1):ncol(index),drop=F] <- t(t(index[,(num.lv.c+1):ncol(index),drop=F])*sig2)}, silent = TRUE)
sigma.lv <- c(sigma.lv, sigma.lv2)
}
}
if((num.lv+num.lv.c) == 0){
qqs <- quantile(params)
uppout <- qqs[4] + (qqs[4]- qqs[2])*2
lowout <- qqs[2] - (qqs[4]- qqs[2])*2
params[params < lowout] <- lowout
params[params > uppout] <- uppout
}
out$params = params
if((num.lv.c+num.RR)>0){
if(num.lv.c>0){
if(!is.matrix(b.lv)) b.lv = as.matrix(b.lv)
b.lv[,1:num.lv.c] = t(t(b.lv[,1:num.lv.c])*abs(sigma.lv[1:num.lv.c]))
}
out$b.lv <- b.lv
if(randomB!=FALSE){
if(starting.val!="zero"&randomB=="LV"){
out$b.lv <- scale(b.lv,center = FALSE)
out$sigmab_lv <- log(attr(scale(b.lv,center=F),"scaled:scale"))
}else if(starting.val=="zero"&randomB=="LV"){
out$sigmab_lv <- rep(0.01,num.lv.c+num.RR)
}else if(randomB=="P"&starting.val!="zero"&(num.lv.c+num.RR)>1){
out$sigmab_lv <- log(apply(b.lv,1,sd))
out$b.lv <- t(scale(t(b.lv),center = FALSE))
}else if(randomB=="P"&starting.val=="zero"|randomB=="P"&(num.lv.c+num.RR)==1){
out$sigmab_lv <- rep(0.01,ncol(lv.X))
}else if(starting.val!="zero"&randomB=="single"){
out$sigmab_lv <- log(sd(b.lv))
out$b.lv <- b.lv/sd(b.lv)
}else if(starting.val=="zero"&randomB=="single"){
out$sigmab_lv <- 0.01
}
# else if(randomB=="all"){
# out$sigmab_lv <- rep(0.01,ncol(lv.X)*(num.RR+num.lv.c))
# }
}
}
if(family=="tweedie")out$Power = Power
if(family=="ZINB")out$ZINB.phi <- ZINB.phi
out$phi <- phi
out$mu <- mu
if(!is.null(TR)) { out$B <- B}
if((num.lv+num.lv.c) > 0) {
out$sigma.lv <- abs(sigma.lv)[1:(num.lv+num.lv.c)]
out$index <- index
out$mu<-out$mu+((out$index)%*%t(out$params[,(ncol(out$params) - (num.lv+num.lv.c) + 1):(ncol(out$params))]%*%diag(out$sigma.lv, length(out$sigma.lv), length(out$sigma.lv))))
if(family == "beta" && starting.val=="res"){
mumax<-apply(abs(out$mu),2,max)
if(any(mumax>5))
out$params[mumax>5,(ncol(out$params) - (num.lv+num.lv.c) + 1):(ncol(out$params))]<-out$params[mumax>5,(ncol(out$params) - (num.lv+num.lv.c) + 1):(ncol(out$params))]/mumax[mumax>5]*5
}
}
if(family == "ordinal") out$zeta <- zeta
if(family=="orderedBeta"){
out$zeta <- zetaOB
}
options(warn = 0)
if(row.eff!=FALSE) {
out$row.params=row.params
if(row.eff=="random") out$sigma=sigma
}
if(!is.null(randomX)){
out$Br <- Br
out$sigmaB <- sigmaB
out$sigmaij <- sigmaij
}
return(out)
}
FAstart <- function(eta, family, y, num.lv = 0, num.lv.c = 0, num.RR = 0, zeta = NULL, zeta.struc = "species", phis = NULL,
jitter.var = 0, resi = NULL, row.eff = FALSE, lv.X, link = NULL, maxit=NULL,max.iter=NULL, Power = NULL, disp.group = NULL, randomB = FALSE, method = "VA", Ntrials = 1){
n<-NROW(y); p <- NCOL(y)
row.params <- NULL # !!!!
b.lv <- NULL
RRcoef <- NULL
RRgamma <- NULL
gamma <- NULL
index.c <- NULL
gamma.c <- NULL
if(family=="orderedBeta"){
y <- y*0.99+0.005
family="beta"
}
#generate starting values for constrained LVs
if(num.lv.c>0&num.lv==0){
if(p>2 && n>2){
# start.fit <- gllvm.TMB(y,lv.X=lv.X,num.lv=0,num.lv.c=0,num.RR = num.lv.c, family=family,starting.val="zero",row.eff=row.eff,sd.errors=F,zeta.struc=zeta.struc,max.iter = 50,maxit=50)
# b.lv <- start.fit$params$LvXcoef
# gamma <- start.fit$params$theta
#eta <- eta + lv.X%*%start.fit$params$LvXcoef%*%t(start.fit$params$theta)
if(family %in% c("poisson", "negative.binomial", "gamma", "exponential","tweedie","ZIP","ZINB")) {
mu <- exp(eta)
}else if(family %in% c("binomial","beta","betaH","orderedBeta")) {
mu <- binomial(link = link)$linkinv(eta)
}else {
mu<-eta
}
if(is.null(resi)){
ds.res <- matrix(NA, n, p)
rownames(ds.res) <- rownames(y)
colnames(ds.res) <- colnames(y)
phis <- phis + 1e-05
for (i in 1:n) {
for (j in 1:p) {
if(!is.na(y[i,j])){
if (family == "ZIP") {
b <- pzip(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzip(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ZINB") {
b <- pzinb(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzinb(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "tweedie") {
b <- fishMod::pTweedie(as.vector(unlist(y[i, j])), mu = mu[i, j], phi = phis[j], p = Power)
a <- min(b,fishMod::pTweedie(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], phi = phis[j], p = Power));
if((as.vector(unlist(y[i, j])) - 1)<0)
a<-0
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "poisson") {
b <- ppois(as.vector(unlist(y[i, j])), mu[i,j])
a <- min(b,ppois(as.vector(unlist(y[i, j])) - 1, mu[i,j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "negative.binomial") {
b <- pnbinom(as.vector(unlist(y[i, j])), mu = mu[i, j], size = 1/phis[j])
a <- min(b,pnbinom(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], size = 1/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "binomial") {
b <- pbinom(as.vector(unlist(y[i, j])), Ntrials, mu[i, j])
a <- min(b,pbinom(as.vector(unlist(y[i, j])) - 1, Ntrials, mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gaussian") {
b <- pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j])
a <- min(b,pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
# ds.res[i, j] <- (y[i, j] - mu[i, j])/phis[j]
}
if (family == "gamma") {
b <- pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j])
a <- min(b,pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "exponential") {
b <- pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j])
a <- min(b,pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "beta") {
b <- pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
a <- min(b,pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j])))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "betaH") {
if(j>(p/2)){
if(y[i, j]==0){
b = 1 - mu[i,j]
a = 0
} else {
b <- 1
a <- 1 - (mu[i,j])
}
} else {
if(y[i, j]==0){
b = 1 - (mu[i,(p/2)+j])
a = 0
} else {
if(y[i,j]==1) y[i,j]=1-1e-16
b <- a <- 1 - (mu[i,(p/2)+j]) + (mu[i,(p/2)+j])*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
}
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "orderedBeta") {
if(y[i, j]==1){
b = 1
a = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2])
} else if(y[i, j]==0){
b = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1])
a = 0
} else {
b <- a <- 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) + (binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2]))*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
u <- try({runif(n = 1, min = a, max = b)})
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ordinal") {
if(zeta.struc == "species"){
probK <- NULL
probK[1] <- pnorm(zeta[j,1]-mu[i,j],log.p = FALSE)
probK[max(y[,j])] <- 1 - pnorm(zeta[j,max(y[,j]) - 1] - mu[i,j])
if(max(y[,j]) > 2) {
j.levels <- 2:(max(y[,j])-1)
for(k in j.levels) { probK[k] <- pnorm(zeta[j,k] - mu[i,j]) - pnorm(zeta[j,k - 1] - mu[i,j]) }
}
probK <- c(0,probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}else{
probK <- NULL
probK[1] <- pnorm(zeta[1] - mu[i, j], log.p = FALSE)
probK[max(y)] <- 1 - pnorm(zeta[max(y) - 1] - mu[i,j])
levels <- 2:(max(y) - min(y))#
for (k in levels) {
probK[k] <- pnorm(zeta[k] - mu[i, j]) - pnorm(zeta[k - 1] - mu[i, j])
}
probK <- c(0, probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}
}
}
}
}
} else {
ds.res <- resi
}
resi <- as.matrix(ds.res); resi[is.na(resi)] <- 0; resi[is.infinite(resi)] <- 0; resi[is.nan(resi)] <- 0
if(n>p){
fa <- try(factanal(resi,factors=num.lv.c,scores = "regression"),silent=T)
if(family=="gaussian"&inherits(fa,"try-error")){
fa <- princomp(resi)
fa$scores <- fa$scores[,1:num.lv.c,drop=F]
fa$loadings <- fa$loadings[,1:num.lv.c,drop=F]
}
if(inherits(fa,"try-error")) stop("Calculating starting values failed. Try centering and scaling your predictors, a smaller 'num.lv.c' value, or change 'starting.val' to 'zero' or 'random'.")
index <- fa$scores
} else if(n<p) {
fa <- try(factanal(t(resi),factors=num.lv.c,scores = "regression"),silent=T)
if(family=="gaussian"&inherits(fa,"try-error")){
fa <- princomp(t(resi))
fa$loadings <- fa$loadings[,1:num.lv.c, drop=F]
fa$scores <- fa$scores[,1:num.lv.c, drop=F]
}
if(inherits(fa,"try-error")) stop("Calculating starting values failed. Try centering and scaling your predictors, a smaller 'num.lv.c' value, or change 'starting.val' to 'zero' or 'random'.")
} else {
tryfit <- TRUE; tryi <- 1
while(tryfit && tryi<5) {
fa <- try(factanal(rbind(resi,rnorm(p,0,0.01)),factors=num.lv.c,scores = "regression"), silent = TRUE)
tryfit <- inherits(fa,"try-error"); tryi <- tryi + 1;
}
if(inherits(fa,"try-error")) {
warning(attr(fa,"condition")$message, "\n Factor analysis for Calculating starting values failed. Try centering and scaling your predictors, a smaller 'num.lv.c' value, or change 'starting.val' to 'zero' or 'random'. Using solution from Principal Component Analysis instead. /n")
fa <- princomp(resi)
}
}
if(n>p){
index<-as.matrix(fa$scores)
gamma <- as.matrix(fa$loadings)[1:p,,drop=F]
}else{
index<-as.matrix(fa$loadings)
gamma <- as.matrix(fa$scores[1:p,,drop=F])
}
index.lm <- lm(index~0+lv.X)
b.lv <- coef(index.lm)
#Ensures independence of LVs with predictors
index <- matrix(residuals.lm(index.lm),ncol=num.lv.c)
colnames(gamma) <- paste("CLV",1:num.lv.c,sep='')
if(num.lv.c>1 && p>2){
qr.gamma <- qr(t(gamma))
gamma.new<-t(qr.R(qr.gamma))
sig <- sign(diag(gamma.new))
gamma <- t(t(gamma.new)*sig)
index<-(index%*%qr.Q(qr.gamma))
b.lv<-b.lv%*%qr.Q(qr.gamma)
b.lv<-t(t(b.lv)*sig)
index <- t(t(index)*sig)
} else if(p>n & num.lv.c>1) {
sdi <- sqrt(diag(cov(index)))
sdt <- sqrt(diag(cov(gamma)))
indexscale <- diag(x = 0.8/sdi, nrow = length(sdi))
index <- index%*%indexscale
gammascale <- diag(x = 1/sdt, nrow = length(sdi))
gamma <- gamma%*%gammascale
}
} else {
b.lv <- matrix(1,ncol=num.lv.c,nrow=ncol(lv.X))
gamma <- matrix(1,p,num.lv.c)
gamma[upper.tri(gamma)]=0
index <- matrix(0,n,num.lv.c)
}
eta <- eta+(index+lv.X%*%b.lv)%*%t(gamma)
}else if(num.lv.c>0&num.lv>0){
if(family!="ordinal"){
zeta.struc<-"species"
}
if(num.lv.c>1)start.fit <- suppressWarnings(gllvm.TMB(y, lv.X = lv.X, num.lv = 0, num.lv.c = num.lv.c, family = family, starting.val = "zero", row.eff = row.eff, sd.errors=F, zeta.struc = zeta.struc, offset = eta, disp.group = disp.group, optimizer = "alabama", method = method, Ntrials = Ntrials))
if(num.lv.c<=1)start.fit <- suppressWarnings(gllvm.TMB(y, lv.X = lv.X, num.lv = 0, num.lv.c = num.lv.c, family = family, starting.val = "zero", row.eff = row.eff, sd.errors=F, zeta.struc = zeta.struc, offset = eta, disp.group = disp.group, optimizer = "nlminb", method = method, Ntrials = Ntrials))
gamma <- start.fit$params$theta
index <- start.fit$lvs
b.lv <- start.fit$params$LvXcoef
# To ensure we do not start off at a point that fully satisfies the constraints
# Especially optimizer="alabama" seems to not like that
if(num.lv.c>1){
mat <- matrix(0,ncol=num.lv.c,nrow=num.lv.c)
mat[upper.tri(mat)]<- 0.01
diag(mat) <- 1
b.lv <- b.lv%*%mat
}
eta <- eta+(index+lv.X%*%b.lv)%*%t(gamma)
if(num.lv>0){
gamma.c <- gamma
index.c <- index
}
}
if(num.RR>0){
#recalculate residual if we have added something to the linear predictor (i.e. num.lv.c)
if(family %in% c("poisson", "negative.binomial", "gamma", "exponential","tweedie","ZIP","ZINB")) {
mu <- exp(eta)
}else if(family %in% c("binomial","beta","betaH","orderedBeta")) {
mu <- binomial(link = link)$linkinv(eta)
}else{
mu <- eta
}
if(num.lv.c>0|is.null(resi)){
ds.res <- matrix(NA, n, p)
rownames(ds.res) <- rownames(y)
colnames(ds.res) <- colnames(y)
phis <- phis + 1e-05
for (i in 1:n) {
for (j in 1:p) {
if(!is.na(y[i,j])){
if (family == "ZIP") {
b <- pzip(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzip(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ZINB") {
b <- pzinb(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzinb(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "tweedie") {
b <- fishMod::pTweedie(as.vector(unlist(y[i, j])), mu = mu[i, j], phi = phis[j], p = Power)
a <- min(b,fishMod::pTweedie(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], phi = phis[j], p = Power));
if((as.vector(unlist(y[i, j])) - 1)<0)
a<-0
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "poisson") {
b <- ppois(as.vector(unlist(y[i, j])), mu[i,j])
a <- min(b,ppois(as.vector(unlist(y[i, j])) - 1, mu[i,j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "negative.binomial") {
b <- pnbinom(as.vector(unlist(y[i, j])), mu = mu[i, j], size = 1/phis[j])
a <- min(b,pnbinom(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], size = 1/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "binomial") {
b <- pbinom(as.vector(unlist(y[i, j])), Ntrials, mu[i, j])
a <- min(b,pbinom(as.vector(unlist(y[i, j])) - 1, Ntrials, mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gaussian") {
b <- pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j])
a <- min(b,pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gamma") {
b <- pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j])
a <- min(b,pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "exponential") {
b <- pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j])
a <- min(b,pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "beta") {
b <- pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
a <- min(b,pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j])))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "betaH") {
if(j>(p/2)){
if(y[i, j]==0){
b = 1 - mu[i,j]
a = 0
} else {
b <- 1
a <- 1 - (mu[i,j])
}
} else {
if(y[i, j]==0){
b = 1 - (mu[i,(p/2)+j])
a = 0
} else {
if(y[i,j]==1) y[i,j]=1-1e-16
b <- a <- 1 - (mu[i,(p/2)+j]) + (mu[i,(p/2)+j])*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
}
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "orderedBeta") {
if(y[i, j]==1){
b = 1
a = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2])
} else if(y[i, j]==0){
b = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1])
a = 0
} else {
b <- a <- 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) + (binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2]))*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
u <- try({runif(n = 1, min = a, max = b)})
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ordinal") {
if(zeta.struc == "species"){
probK <- NULL
probK[1] <- pnorm(zeta[j,1]-mu[i,j],log.p = FALSE)
probK[max(y[,j])] <- 1 - pnorm(zeta[j,max(y[,j]) - 1] - mu[i,j])
if(max(y[,j]) > 2) {
j.levels <- 2:(max(y[,j])-1)
for(k in j.levels) { probK[k] <- pnorm(zeta[j,k] - mu[i,j]) - pnorm(zeta[j,k - 1] - mu[i,j]) }
}
probK <- c(0,probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}else{
probK <- NULL
probK[1] <- pnorm(zeta[1] - mu[i, j], log.p = FALSE)
probK[max(y)] <- 1 - pnorm(zeta[max(y) - 1] - mu[i,j])
levels <- 2:(max(y) - min(y))#
for (k in levels) {
probK[k] <- pnorm(zeta[k] - mu[i, j]) - pnorm(zeta[k - 1] - mu[i, j])
}
probK <- c(0, probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}
}
}
}
}
} else {
ds.res <- resi
}
resi <- as.matrix(ds.res); resi[is.na(resi)] <- 0; resi[is.infinite(resi)] <- 0; resi[is.nan(resi)] <- 0
}
if(num.RR>0){
#Alternative for RRcoef is factor analysis of the predictors.
RRmod <- lm(resi~0+lv.X)
beta <- t(coef(RRmod))
# if(ncol(beta)>1&(randomB==FALSE)){
# betaSD=apply(beta,1,sd)
# beta=beta/betaSD
# }
qr.beta=qr(t(beta))
R=t(qr.R(qr.beta))[,1:num.RR,drop=F]
# R=R[,1:num.RR,drop=F]/sqrt(num.RR*nrow(RRmod$coefficients))
Q=t(qr.Q(qr.beta))[1:num.RR,,drop=F]
# To ensure we do not start off at a point that fully satisfies the constraints
# Especially optimizer="alabama" seems to not like that
if(num.RR>1){
mat <- matrix(0,ncol=num.RR,nrow=num.RR)
mat[upper.tri(mat)]<- 0.01
diag(mat) <- 1
Q <- (mat%*%Q)
}
# RRcoef = Q*sqrt(num.RR*nrow(RRmod$coefficients))
RRcoef <- Q
sgn=t(sign(diag(R)))
if(num.RR>1){
RRgamma = R%*%diag(c(sgn))
RRcoef = t(diag(c(sgn))%*%RRcoef)
}else{
RRgamma = R*c(sgn)
RRcoef = t(c(sgn)*RRcoef)
}
#transfer RRgamma diagonals to RRcoef
RRcoef <- t(t(RRcoef)*diag(RRgamma))
RRgamma<-t(t(RRgamma)/diag(RRgamma))
# qr.gam <- qr(coef(RRmod))
# RRcoef <- qr.Q(qr.gam)[,1:num.RR]
# RRcoef <- t(t(RRcoef)*diag(qr.R(qr.gam))[1:num.RR])
# RRgamma <- t(qr.R(qr.gam)/diag(qr.R(qr.gam)))[,1:num.RR]
eta <- eta + lv.X%*%RRcoef%*%t(RRgamma)
}
if((num.lv.c+num.RR)>0&num.lv>0){
# recalculate if we have added something to the linear predictor (i.e. num.lv.c. or num.RR)
if(family %in% c("poisson", "negative.binomial", "gamma", "exponential","tweedie","ZIP","ZINB")) {
mu <- exp(eta)
}else if(family %in% c("binomial","beta","betaH","orderedBeta")) {
mu <- binomial(link = link)$linkinv(eta)
}else{
mu <- eta
}
if(is.null(resi)){
ds.res <- matrix(NA, n, p)
rownames(ds.res) <- rownames(y)
colnames(ds.res) <- colnames(y)
phis <- phis + 1e-05
for (i in 1:n) {
for (j in 1:p) {
if(!is.na(y[i,j])){
if (family == "ZIP") {
b <- pzip(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzip(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ZINB") {
b <- pzinb(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzinb(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "tweedie") {
b <- fishMod::pTweedie(as.vector(unlist(y[i, j])), mu = mu[i, j], phi = phis[j], p = Power)
a <- min(b,fishMod::pTweedie(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], phi = phis[j], p = Power));
if((as.vector(unlist(y[i, j])) - 1)<0)
a<-0
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "poisson") {
b <- ppois(as.vector(unlist(y[i, j])), mu[i,j])
a <- min(b,ppois(as.vector(unlist(y[i, j])) - 1, mu[i,j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "negative.binomial") {
b <- pnbinom(as.vector(unlist(y[i, j])), mu = mu[i, j], size = 1/phis[j])
a <- min(b,pnbinom(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], size = 1/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "binomial") {
b <- pbinom(as.vector(unlist(y[i, j])), Ntrials, mu[i, j])
a <- min(b,pbinom(as.vector(unlist(y[i, j])) - 1, Ntrials, mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gaussian") {
b <- pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j])
a <- min(b,pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gamma") {
b <- pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j])
a <- min(b,pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "exponential") {
b <- pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j])
a <- min(b,pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "beta") {
b <- pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
a <- min(b,pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j])))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "betaH") {
if(j>(p/2)){
if(y[i, j]==0){
b = 1 - mu[i,j]
a = 0
} else {
b <- 1
a <- 1 - (mu[i,j])
}
} else {
if(y[i, j]==0){
b = 1 - (mu[i,(p/2)+j])
a = 0
} else {
if(y[i,j]==1) y[i,j]=1-1e-16
b <- a <- 1 - (mu[i,(p/2)+j]) + (mu[i,(p/2)+j])*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
}
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "orderedBeta") {
if(y[i, j]==1){
b = 1
a = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2])
} else if(y[i, j]==0){
b = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1])
a = 0
} else {
b <- a <- 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) + (binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2]))*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
u <- try({runif(n = 1, min = a, max = b)})
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ordinal") {
if(zeta.struc == "species"){
probK <- NULL
probK[1] <- pnorm(zeta[j,1]-mu[i,j],log.p = FALSE)
probK[max(y[,j])] <- 1 - pnorm(zeta[j,max(y[,j]) - 1] - mu[i,j])
if(max(y[,j]) > 2) {
j.levels <- 2:(max(y[,j])-1)
for(k in j.levels) { probK[k] <- pnorm(zeta[j,k] - mu[i,j]) - pnorm(zeta[j,k - 1] - mu[i,j]) }
}
probK <- c(0,probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}else{
probK <- NULL
probK[1] <- pnorm(zeta[1] - mu[i, j], log.p = FALSE)
probK[max(y)] <- 1 - pnorm(zeta[max(y) - 1] - mu[i,j])
levels <- 2:(max(y) - min(y))#
for (k in levels) {
probK[k] <- pnorm(zeta[k] - mu[i, j]) - pnorm(zeta[k - 1] - mu[i, j])
}
probK <- c(0, probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}
}
}
}
}
} else {
ds.res <- resi
}
resi <- as.matrix(ds.res); resi[is.na(resi)] <- 0; resi[is.infinite(resi)] <- 0; resi[is.nan(resi)] <- 0
if(p>2 && n>2){
if(any(is.nan(resi))){stop("Method 'res' for starting values can not be used, when glms fit too poorly to the data. Try other starting value methods 'zero' or 'random' or change the model.")}
if(n>p){
fa <- try(factanal(resi,factors=num.lv,scores = "regression"))
if(family=="gaussian"&inherits(fa,"try-error")){
fa <- princomp(resi)
fa$scores <- fa$scores[,1:num.lv,drop=F]
fa$loadings <- fa$loadings[,1:num.lv,drop=F]
}
if(inherits(fa,"try-error")) stop("Calculating starting values failed. Maybe too many latent variables. Try smaller 'num.lv' value or change 'starting.val' to 'zero' or 'random'.")
gamma<-matrix(fa$loadings,p,num.lv)
index <- fa$scores
} else if(n<p) {
fa <- try(factanal(t(resi),factors=num.lv,scores = "regression"))
if(family=="gaussian"&inherits(fa,"try-error")){
fa <- princomp(t(resi))
fa$loadings <- fa$loadings[,1:num.lv,drop=F]
fa$scores <- fa$scores[,1:num.lv,drop=F]
}
if(inherits(fa,"try-error")) stop("Calculating starting values failed. Maybe too many latent variables. Try smaller 'num.lv' value or change 'starting.val' to 'zero' or 'random'.")
gamma<-fa$scores
index <- matrix(fa$loadings,n,num.lv)
} else {
tryfit <- TRUE; tryi <- 1
while(tryfit && tryi<5) {
fa <- try(factanal(rbind(resi,rnorm(p,0,0.01)),factors=num.lv,scores = "regression"), silent = TRUE)
tryfit <- inherits(fa,"try-error"); tryi <- tryi + 1;
}
if(tryfit) {
warning(attr(fa,"condition")$message, "\n Factor analysis for calculating starting values failed. Maybe too many latent variables. Try smaller 'num.lv' value or change 'starting.val' to 'zero' or 'random'. Using solution from Principal Component Analysis instead./n")
pr <- princomp(resi)
gamma<-matrix(pr$loadings[,1:num.lv],p,num.lv)
index<-matrix(pr$scores[,1:num.lv],n,num.lv)
}else{
gamma<-matrix(fa$loadings,p,num.lv)
index <- fa$scores[1:num.lv,]
}
}
} else {
gamma <- matrix(1,p,num.lv)
gamma[upper.tri(gamma)]=0
index <- matrix(0,n,num.lv)
}
index <- cbind(index.c,index)
gamma<-cbind(gamma.c,gamma)
}
if(num.lv>0&(num.lv.c+num.RR)==0){
if(family %in% c("poisson", "negative.binomial", "gamma", "exponential","tweedie","ZIP","ZINB")) {
mu <- exp(eta)
}else if(family %in% c("binomial","beta", "betaH", "orderedBeta")) {
mu <- binomial(link = link)$linkinv(eta)
}else{
mu <- eta
}
if(is.null(resi)){
ds.res <- matrix(NA, n, p)
rownames(ds.res) <- rownames(y)
colnames(ds.res) <- colnames(y)
phis <- phis + 1e-05
for (i in 1:n) {
for (j in 1:p) {
if(!is.na(y[i,j])){
if (family == "ZIP") {
b <- pzip(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzip(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ZINB") {
b <- pzinb(as.vector(unlist(y[i, j])), mu = mu[i, j], sigma = phis[j])
a <- min(b,pzinb(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], sigma = phis[j]))
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "tweedie") {
b <- fishMod::pTweedie(as.vector(unlist(y[i, j])), mu = mu[i, j], phi = phis[j], p = Power)
a <- min(b,fishMod::pTweedie(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], phi = phis[j], p = Power));
if((as.vector(unlist(y[i, j])) - 1)<0)
a<-0
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "poisson") {
b <- ppois(as.vector(unlist(y[i, j])), mu[i,j])
a <- pmin(b,ppois(as.vector(unlist(y[i, j])) - 1, mu[i,j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "negative.binomial") {
b <- pnbinom(as.vector(unlist(y[i, j])), mu = mu[i, j], size = 1/phis[j])
a <- min(b,pnbinom(as.vector(unlist(y[i, j])) - 1, mu = mu[i, j], size = 1/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "binomial") {
b <- pbinom(as.vector(unlist(y[i, j])), Ntrials, mu[i, j])
a <- min(b,pbinom(as.vector(unlist(y[i, j])) - 1, Ntrials, mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gaussian") {
b <- pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j])
a <- min(b,pnorm(as.vector(unlist(y[i, j])), mu[i, j], sd = phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "gamma") {
b <- pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j])
a <- min(b,pgamma(as.vector(unlist(y[i, j])), shape = phis[j], scale = mu[i, j]/phis[j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "exponential") {
b <- pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j])
a <- min(b,pexp(as.vector(unlist(y[i, j])), rate = 1/mu[i, j]))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "beta") {
b <- pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
a <- min(b,pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j])))
u <- runif(n = 1, min = a, max = b)
ds.res[i, j] <- qnorm(u)
}
if (family == "betaH") {
if(j>(p/2)){
if(y[i, j]==0){
b = 1 - mu[i,j]
a = 0
} else {
b <- 1
a <- 1 - (mu[i,j])
}
} else {
if(y[i, j]==0){
b = 1 - (mu[i,(p/2)+j])
a = 0
} else {
if(y[i,j]==1) y[i,j]=1-1e-16
b <- a <- 1 - (mu[i,(p/2)+j]) + (mu[i,(p/2)+j])*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
}
u <- runif(n = 1, min = a, max = b)
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "orderedBeta") {
if(y[i, j]==1){
b = 1
a = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2])
} else if(y[i, j]==0){
b = 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1])
a = 0
} else {
b <- a <- 1 - binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) + (binomial(link = link)$linkinv(eta[i,j] - zeta[j,1]) - binomial(link = link)$linkinv(eta[i,j] - zeta[j,2]))*pbeta(as.vector(unlist(y[i, j])), shape1 = phis[j]*mu[i, j], shape2 = phis[j]*(1-mu[i, j]))
}
u <- try({runif(n = 1, min = a, max = b)})
if(u==1) u=1-1e-16
if(u==0) u=1e-16
ds.res[i, j] <- qnorm(u)
}
if (family == "ordinal") {
if(zeta.struc == "species"){
probK <- NULL
probK[1] <- pnorm(zeta[j,1]-mu[i,j],log.p = FALSE)
probK[max(y[,j])] <- 1 - pnorm(zeta[j,max(y[,j]) - 1] - mu[i,j])
if(max(y[,j]) > 2) {
j.levels <- 2:(max(y[,j])-1)
for(k in j.levels) { probK[k] <- pnorm(zeta[j,k] - mu[i,j]) - pnorm(zeta[j,k - 1] - mu[i,j]) }
}
probK <- c(0,probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y[,j])==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}else{
probK <- NULL
probK[1] <- pnorm(zeta[1] - mu[i, j], log.p = FALSE)
probK[max(y)] <- 1 - pnorm(zeta[max(y) - 1] - mu[i,j])
levels <- 2:(max(y) - min(y))#
for (k in levels) {
probK[k] <- pnorm(zeta[k] - mu[i, j]) - pnorm(zeta[k - 1] - mu[i, j])
}
probK <- c(0, probK)
cumsum.b <- sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0) + 1)])
cumsum.a <- min(cumsum.b, sum(probK[1:(y[i,j]+ifelse(min(y)==0,1,0))]))
u <- runif(n = 1, min = cumsum.a, max = cumsum.b)
if (abs(u - 1) < 1e-05)
u <- 1
if (abs(u - 0) < 1e-05)
u <- 0
ds.res[i, j] <- qnorm(u)
}
}
}
}
}
} else {
ds.res <- resi
}
resi <- as.matrix(ds.res); resi[is.na(resi)] <- 0; resi[is.infinite(resi)] <- 0; resi[is.nan(resi)] <- 0
if(p>2 && n>2){
if(any(is.nan(resi))){stop("Method 'res' for starting values can not be used, when glms fit too poorly to the data. Try other starting value methods 'zero' or 'random' or change the model.")}
if(n>p){
fa <- try(factanal(resi,factors=num.lv,scores = "regression"))
if(inherits(fa,"try-error")) stop("Factor analysis for calculating starting values failed. Maybe too many latent variables. Try smaller 'num.lv' value or change 'starting.val' to 'zero' or 'random'.")
gamma<-matrix(fa$loadings,p,num.lv)
index <- fa$scores
} else if(n<p) {
fa <- try(factanal(t(resi),factors=num.lv,scores = "regression"))
if(inherits(fa,"try-error")) stop("Factor analysis for calculating starting values failed. Maybe too many latent variables. Try smaller 'num.lv' value or change 'starting.val' to 'zero' or 'random'.")
gamma<-fa$scores
index <- matrix(fa$loadings,n,num.lv)
} else {
tryfit <- TRUE; tryi <- 1
while(tryfit && tryi<5) {
fa <- try(factanal(rbind(resi,rnorm(p,0,0.01)),factors=num.lv,scores = "regression"), silent = TRUE)
tryfit <- inherits(fa,"try-error"); tryi <- tryi + 1;
}
if(tryfit) {
warning(attr(fa,"condition")$message, "\n Factor analysis for calculating starting values failed. Maybe too many latent variables. Try smaller 'num.lv' value or change 'starting.val' to 'zero' or 'random'. Using solution from Principal Component Analysis instead./n")
fa <- princomp(resi)
}else{
gamma<-matrix(fa$loadings[,1:num.lv],p,num.lv)
index<-matrix(fa$scores[,1:num.lv],n,num.lv)
}
}
} else {
gamma <- matrix(1,p,num.lv)
gamma[upper.tri(gamma)]=0
index <- matrix(0,n,num.lv)
}
b.lv <- NULL
}
if(row.eff == "random") { # !!!!
row.params <- index[,(num.lv+num.lv.c)]* (sd(matrix(index[,(num.lv+num.lv.c)])%*%matrix(gamma[, (num.lv+num.lv.c)], nrow=1)))/sd(index[,(num.lv+num.lv.c)])
index <- as.matrix(index[,-(num.lv+num.lv.c)])
gamma <- as.matrix(gamma[,-(num.lv+num.lv.c)])
num.lv <- num.lv-1
}
if((num.lv+num.lv.c)>0){
if((num.lv+num.lv.c)>1 && p>2){
if(num.lv.c==0&num.lv>0|num.lv==0&num.lv.c>0){
qr.gamma <- qr(t(gamma))
gamma.new<-t(qr.R(qr.gamma))
sig <- sign(diag(gamma.new))
if(num.lv.c>0){
b.lv<-b.lv%*%qr.Q(qr.gamma)
b.lv <- t(t(b.lv)*sig)
}
gamma <- t(t(gamma.new)*sig)
index<-(index%*%qr.Q(qr.gamma))
if(num.lv.c>0)b.lv<-b.lv%*%qr.Q(qr.gamma)
index <- t(t(index)*sig)
}else{
qr.gamma1 <- qr(t(gamma[,1:num.lv.c,drop=F]))
qr.gamma2 <- qr(t(gamma[,(num.lv.c+1):ncol(gamma),drop=F]))
gamma.new1<-t(qr.R(qr.gamma1))
gamma.new2<-t(qr.R(qr.gamma2))
sig1 <- sign(diag(gamma.new1))
b.lv<-b.lv%*%qr.Q(qr.gamma1)
b.lv <- t(t(b.lv)*sig1)
sig2 <- sign(diag(gamma.new2))
gamma.new1 <- t(t(gamma.new1)*sig1)
gamma.new2 <- t(t(gamma.new2)*sig2)
gamma <- cbind(gamma.new1,gamma.new2)
index[,1:num.lv.c]<-(index[,1:num.lv.c,drop=F]%*%qr.Q(qr.gamma1)[1:num.lv.c,1:num.lv.c])
index[,(num.lv.c+1):ncol(index)]<-(index[,(num.lv.c+1):ncol(index)]%*%qr.Q(qr.gamma2))
index[,1:num.lv.c] <- t(t(index[,1:num.lv.c])*sig1[1:num.lv.c])
index[,(num.lv.c+1):ncol(index)] <- t(t(index[,(num.lv.c+1):ncol(index)])*sig2)
}
} else {
if((num.lv.c>0&num.lv==0)|(num.lv>0&num.lv.c==0)){
sig <- sign(diag(gamma))
if(num.lv.c>0)b.lv <- t(t(b.lv)*sig)
gamma <- t(t(gamma)*sig)
index <- t(t(index)*sig)
}else if(num.lv.c>0&num.lv>0){
sig1 <- sign(diag(gamma[,1:num.lv.c,drop=F]));
sig2 <- sign(diag(gamma[,(num.lv.c+1):ncol(gamma),drop=F]))
b.lv <- t(t(b.lv)*sig1)
gamma[,1:num.lv.c] <- t(t(gamma[,1:num.lv.c,drop=F])*sig1)
gamma[,(num.lv.c+1):ncol(gamma)] <- t(t(gamma[,(num.lv.c+1):ncol(gamma),drop=F])*sig2)
index[,1:num.lv.c] <- t(t(index[,1:num.lv.c,drop=F])*sig1)
index[,(num.lv.c+1):ncol(index)] <- t(t(index[,(num.lv.c+1):ncol(index),drop=F])*sig2)
}
}
if(p>n) {
if(num.lv==0&num.lv.c>0|num.lv.c==0&num.lv>0){
sdi <- sqrt(diag(cov(index)))
sdt <- sqrt(diag(cov(gamma)))
indexscale <- diag(x = 0.8/sdi, nrow = length(sdi))
index <- index%*%indexscale
gammascale <- diag(x = 1/sdt, nrow = length(sdi))
gamma <- gamma%*%gammascale
}else if(num.lv>0&num.lv.c>0){
sdi <- sqrt(diag(cov(index[,1:num.lv.c,drop=F])))
sdt <- sqrt(diag(cov(gamma[,1:num.lv.c,drop=F])))
indexscale <- diag(x = 0.8/sdi, nrow = length(sdi))
index[,1:num.lv.c] <- index[,1:num.lv.c,drop=F]%*%indexscale
gammascale <- diag(x = 1/sdt, nrow = length(sdi))
gamma[,1:num.lv.c] <- gamma[,1:num.lv.c,drop=F]%*%gammascale
sdi <- sqrt(diag(cov(index[,(num.lv.c+1):ncol(index),drop=F])))
sdt <- sqrt(diag(cov(gamma[,(num.lv.c+1):ncol(index),drop=F])))
indexscale <- diag(x = 0.8/sdi, nrow = length(sdi))
index[,(num.lv.c+1):ncol(index)] <- index[,(num.lv.c+1):ncol(index),drop=F]%*%indexscale
gammascale <- diag(x = 1/sdt, nrow = length(sdi))
gamma[,(num.lv.c+1):ncol(index)] <- gamma[,(num.lv.c+1):ncol(index),drop=F]%*%gammascale
}
} else {
if(num.lv==0|num.lv.c==0){
sdi <- sqrt(diag(cov(index)))
sdt <- sqrt(diag(cov(gamma)))
if( any(is.na(c(sdi,sdt)))|any(c(sdi,sdt)==0) ) { sdi<-rep(1,length(sdi)); sdt<-rep(1,length(sdt)) }
indexscale <- diag(x = 1/sdi, nrow = length(sdi))
index <- index%*%indexscale
gammascale <- diag(x = 0.8/sdt, nrow = length(sdi))
gamma <- gamma%*%gammascale
}else if(num.lv>0&num.lv.c>0){
sdi <- sqrt(diag(cov(index[,1:num.lv.c,drop=F])))
sdt <- sqrt(diag(cov(gamma[,1:num.lv.c,drop=F])))
indexscale <- diag(x = 1/sdi, nrow = length(sdi))
index[,1:num.lv.c] <- index[,1:num.lv.c,drop=F]%*%indexscale
gammascale <- diag(x = 0.8/sdt, nrow = length(sdi))
gamma[,1:num.lv.c] <- gamma[,1:num.lv.c,drop=F]%*%gammascale
sdi <- sqrt(diag(cov(index[,(num.lv.c+1):ncol(index),drop=F])))
sdt <- sqrt(diag(cov(gamma[,(num.lv.c+1):ncol(index),drop=F])))
if( any(is.na(c(sdi,sdt)))|any(c(sdi,sdt)==0) ) { sdi<-rep(1,length(sdi)); sdt<-rep(1,length(sdt)) }
indexscale <- diag(x = 1/sdi, nrow = length(sdi))
index[,(num.lv.c+1):ncol(index)] <- index[,(num.lv.c+1):ncol(index),drop=F]%*%indexscale
gammascale <- diag(x = 0.8/sdt, nrow = length(sdi))
gamma[,(num.lv.c+1):ncol(index)] <- gamma[,(num.lv.c+1):ncol(index),drop=F]%*%gammascale
}
}
index <- index + MASS::mvrnorm(n, rep(0, num.lv+num.lv.c),diag(num.lv+num.lv.c)*jitter.var);
}else{
index <- NULL
if(num.RR==0){
gamma <- NULL
b.lv <- NULL
}else{
gamma <- RRgamma
}
}
if((num.lv.c+num.lv)>0&num.RR>0){
gammaC <- NULL
gammaU <- NULL
if(num.lv.c>0)gammaC <- gamma[,1:num.lv.c]
if(num.lv>0)gammaU <- gamma[,(ncol(gamma)-num.lv+1):ncol(gamma)]
gamma <- cbind(gammaC,RRgamma,gammaU)
}
return(list(index = index, gamma = gamma, row.params =row.params, b.lv = cbind(b.lv,RRcoef)))
}
## Calculates information matrix based on scores of VA log-likelihood
calc.infomat <- function(theta = NULL, beta0=NULL, env = NULL, row.params = NULL, vameans = NULL, lambda=NULL, phi = NULL, zeta = NULL, num.lv, family, Lambda.struc, row.eff, y, X = NULL,TR = NULL, Xd=NULL, offset=0, B=NULL) {
n <- nrow(y); p <- ncol(y)
if(!is.null(X)) num.X <- ncol(X)
if(!is.null(TR)) num.T <- ncol(TR)
trial.size <- 1
if(family == "ordinal") {
max.levels <- max(y)
if(any(max.levels == 1) || all(max.levels == 2)) stop("Ordinal data requires all columns to have at least has two levels. If all columns only have two levels, please use family == binomial instead. Thanks")
}
if (is.null(offset)) offset <- matrix(0, nrow = n, ncol = p)
if (!is.null(Xd)) nd <- ncol(Xd)
## score of VA logL wrt model parameters
grad.all.mod <- function(x,vameans=NULL,lambda=NULL) {
x2 <- x
new.vameans <- new.lambda <- new.theta <- NULL
if(num.lv > 0) {
new.lambda <- lambda
new.vameans <- vameans
new.theta <- matrix(c(x2[1:(p * num.lv)]),p,num.lv); x2 <- x2[-(1:(p * num.lv))]
new.theta[upper.tri(new.theta)] <- 0
}
new.beta0 <- x2[1:p]; x2 <- x2[-(1:p)]
new.env <- NULL
if(!is.null(X)) {
if(is.null(TR)) { new.env <- matrix(x2[1:(p * num.X)],p,num.X); x2 <- x2[-(1:(p * num.X))]
} else { B=x2[1:nd]; x2 <- x2[-(1:nd)]}
}
new.row.params <- NULL; if(row.eff) { new.row.params <- x2[1:n]; x2 <- x2[-(1:n)] }
new.phi <- NULL; if(family == "negative.binomial") { new.phi <- x2[1:p]; x2 <- x2[-(1:p)] }
new.zeta <- NULL; if(family == "ordinal") {
new.zeta <- matrix(NA,p,max.levels - 2)
for(j in 1:p) { if(max(y[,j]) > 2) { new.zeta[j,1:(max(y[,j])-2)] <- x2[1:(max(y[,j])-2)]; x2 <- x2[-(1:(max(y[,j])-2))] } }
new.zeta <- cbind(0,new.zeta);
}
eta.mat <- matrix(new.beta0,n,p,byrow=TRUE) + offset
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + X %*% t(new.env)
if(!is.null(TR)) eta.mat <- eta.mat + matrix((Xd %*% B),n,p)
if(num.lv > 0) eta.mat <- eta.mat + new.vameans %*% t(new.theta)
if(row.eff) eta.mat <- eta.mat + matrix(new.row.params,n,p)
grad.theta <- grad.beta0 <- grad.env <- grad.B <- grad.phi <- grad.row.params <- grad.zeta <- NULL;
if(Lambda.struc == "unstructured" & num.lv > 0) {
Lambda.theta <- array(NA,dim=c(p,n,num.lv))
for(i2 in 1:n) { Lambda.theta[,i2,] <- new.theta%*%new.lambda[i2,,] }
}
if(family=="poisson") {
if(num.lv > 0) eta.mat <- eta.mat + calc.quad(new.lambda,new.theta,Lambda.struc)$mat
grad.beta0 <- colSums(y-exp(eta.mat))
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(y,1,new.vameans[,l],"*") - sweep(t(Lambda.theta[,,l]),1,new.vameans[,l],"+") * exp(eta.mat)
if(Lambda.struc == "diagonal") sum1 <- sweep(y,1,new.vameans[,l],"*") - sweep((new.lambda[,l]%*%t(new.theta[,l])),1,new.vameans[,l],"+")*exp(eta.mat)
grad.theta <- c(grad.theta,colSums(sum1))
}
}
if(!is.null(X) && is.null(TR)) {
for (l in 1:num.X) {
sum1 <- sweep(y-exp(eta.mat),1,X[,l],"*")
grad.env <- c(grad.env,colSums(sum1))
}
}
if(!is.null(TR)) {
sum1 <- c(y-exp(eta.mat)) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) { grad.row.params <- rowSums(y - exp(eta.mat)) }
}
if(family=="negative.binomial") {
eta.mat <- matrix(new.beta0,n,p,byrow=TRUE) + offset
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + X %*% t(new.env)
if(!is.null(TR)) eta.mat <- eta.mat + matrix((Xd %*% B),n,p)
if(row.eff) eta.mat <- eta.mat + matrix(new.row.params,n,p)
if(num.lv > 0) eta.mat <- eta.mat + new.vameans %*% t(new.theta)
phi.mat <- matrix(new.phi,n,p,byrow=TRUE)
mu.mat.noquad <- eta.mat
if(num.lv > 0) eta.mat <- eta.mat - calc.quad(new.lambda,new.theta,Lambda.struc)$mat
grad.beta0 <- colSums(-phi.mat - (y + phi.mat) * (-phi.mat / (exp(eta.mat) + phi.mat)))
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(-phi.mat,1,new.vameans[,l],"*") + sweep(t(Lambda.theta[,,l]),1,-new.vameans[,l],"+") * (-(y + phi.mat) * (phi.mat / (exp(eta.mat) + phi.mat)))
if(Lambda.struc == "diagonal") sum1 <- sweep(-phi.mat,1,new.vameans[,l],"*") + sweep((new.lambda[,l] %*% t(new.theta[,l])),1,-new.vameans[,l],"+") * (-(y + phi.mat) * (phi.mat / (exp(eta.mat) + phi.mat)))
grad.theta <- c(grad.theta,colSums(sum1))
}
}
if(!is.null(X) && is.null(TR)) {
for(l in 1:num.X) {
sum1 <- sweep(-phi.mat - (y + phi.mat) * (-phi.mat / (exp(eta.mat) + phi.mat)),1,X[,l],"*")
grad.env <- c(grad.env,colSums(sum1))
}
}
if(!is.null(TR)) {
sum1 <- c(-phi.mat - (y + phi.mat) * (-phi.mat / (exp(eta.mat) + phi.mat))) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) grad.row.params <- rowSums(-phi.mat - (y + phi.mat) * (-phi.mat/(exp(eta.mat) + phi.mat)))
grad.phi <- colSums(-mu.mat.noquad + 1 + log(phi.mat) - digamma(phi.mat) - log(1 + phi.mat / exp(eta.mat)) - (y + phi.mat)/(phi.mat + exp(eta.mat)) + digamma(y + phi.mat))
}
if(family=="binomial") {
grad.beta0 <- colSums(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat))),na.rm=TRUE)
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat))/(pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10),1,new.vameans[,l],"*") - t(Lambda.theta[,,l])
if(Lambda.struc == "diagonal") sum1 <- sweep(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10),1,new.vameans[,l],"*") - (new.lambda[,l] %*% t(new.theta[,l]))
grad.theta <- c(grad.theta,colSums(sum1,na.rm=TRUE)) }
}
if(!is.null(X) && is.null(TR)) {
for(l in 1:num.X) {
sum1 <- sweep(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10),1,X[,l],"*")
grad.env <- c(grad.env,colSums(sum1))
}
}
if(!is.null(TR)) {
sum1 <- c(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10)) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) { grad.row.params <- rowSums((dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)))),na.rm=TRUE) }
}
if(family=="ordinal") {
grad.zeta <- matrix(0,p,ncol(new.zeta))
deriv.trunnorm <- matrix(0,n,p)
for(j in 1:p) {
deriv.trunnorm[y[,j] == 1,j] <- -dnorm(new.zeta[j,1] - eta.mat[y[,j] == 1,j]) / pnorm(new.zeta[j,1] - eta.mat[y[,j] == 1,j])
deriv.trunnorm[y[,j] == max(y[,j]),j] <- dnorm(new.zeta[j,max(y[,j]) - 1] - eta.mat[y[,j] == max(y[,j]),j]) / (1 - pnorm(new.zeta[j,max(y[,j]) - 1] - eta.mat[y[,j] == max(y[,j]),j]))
j.levels <- (1:max(y[,j]))[-c(1,max(y[,j]))]
for(k in j.levels) {
deriv.trunnorm[y[,j] == k,j] <- (-dnorm(new.zeta[j,k] - eta.mat[y[,j] == k,j]) + dnorm(new.zeta[j,k - 1] - eta.mat[y[,j] == k,j])) / (pnorm(new.zeta[j,k] - eta.mat[y[,j] == k,j]) - pnorm(new.zeta[j,k-1] - eta.mat[y[,j] == k,j])) }
}
deriv.trunnorm[!is.finite(deriv.trunnorm)] <- 0
grad.beta0 <- colSums(deriv.trunnorm,na.rm=TRUE)
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(deriv.trunnorm,1,new.vameans[,l],"*") - t(Lambda.theta[,,l])
if(Lambda.struc == "diagonal") sum1 <- sweep(deriv.trunnorm,1,new.vameans[,l],"*") - (new.lambda[,l] %*% t(new.theta[,l]))
grad.theta <- c(grad.theta,colSums(sum1,na.rm=TRUE))
}
}
if(!is.null(X) && is.null(TR)) {
for(l in 1:num.X) { sum1 <- sweep(deriv.trunnorm,1,X[,l],"*");
grad.env <- c(grad.env,colSums(sum1,na.rm=TRUE))
}
}
if(!is.null(TR)) {
sum1 <- c(deriv.trunnorm) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) { grad.row.params <- rowSums(deriv.trunnorm,na.rm=TRUE) }
for(j in 1:p) {
zeta0 <- new.zeta[j,]
grad.zeta[j,max(y[,j]) - 1] <- grad.zeta[j,max(y[,j]) - 1] - sum(dnorm(zeta0[max(y[,j]) - 1] - eta.mat[which(y[,j] == max(y[,j])) ,j]) / (1 - pnorm(zeta0[max(y[,j]) - 1] - eta.mat[which(y[,j] == max(y[,j])),j])))
j.levels <- (1:max(y[,j]))[-c(1,max(y[,j]))]
for(k in j.levels) {
grad.zeta[j,k] <- grad.zeta[j,k] + sum(dnorm(zeta0[k] - eta.mat[y[,j] == k,j]) / (pnorm(zeta0[k] - eta.mat[y[,j] == k,j]) - pnorm(zeta0[k - 1] - eta.mat[y[,j] == k,j])))
grad.zeta[j,k - 1] <- grad.zeta[j,k - 1] - sum(dnorm(zeta0[k - 1] - eta.mat[y[,j] == k,j]) / (pnorm(zeta0[k] - eta.mat[y[,j] == k,j]) - pnorm(zeta0[k - 1] - eta.mat[y[,j] == k,j])))
}
}
grad.zeta <- grad.zeta[,-1]
if(length(unique(apply(y,2,max))) > 1) { grad.zeta <- c(grad.zeta[-which(grad.zeta == 0)]) }
}
if(num.lv>0){
grad.theta=matrix(grad.theta,nrow=p,ncol=num.lv)
grad.theta[upper.tri(grad.theta)] <- 0}
if(row.eff) grad.row.params[1]=0;
return(c(grad.theta, grad.beta0, grad.env, grad.B, grad.row.params, grad.phi, grad.zeta))
}
## score of VA logL wrt vameans_sel.i and lambda_sel.i
grad.all.var <- function(x,vameans=NULL,lambda=NULL,mod.x,sel.i) {
x2 <- mod.x
if(num.lv > 0) {
new.vameans <- vameans; new.lambda <- lambda
new.theta <- matrix(c(x2[1:(p * num.lv)]),p,num.lv); x2 <- x2[-(1:(p * num.lv))]
new.theta[upper.tri(new.theta)] <- 0
}
new.beta0 <- x2[1:p]; x2 <- x2[-(1:p)]; new.env <- NULL
if(!is.null(X)) {
if(is.null(TR)) { new.env <- matrix(x2[1:(p * num.X)],p,num.X); x2 <- x2[-(1:(p * num.X))]
} else { B=x2[1:nd]; x2 <- x2[-(1:nd)]}
}
new.row.params <- NULL; if(row.eff) { new.row.params <- x2[1:n]; x2 <- x2[-(1:n)] }
new.phi <- NULL; if(family == "negative.binomial") { new.phi <- x2[1:p]; x2 <- x2[-(1:p)] }
new.zeta <- NULL
if(family == "ordinal") {
new.zeta <- matrix(NA,p,max.levels-2);
for(j in 1:p) { if(max(y[,j]) > 2) { new.zeta[j,1:(max(y[,j])-2)] <- x2[1:(max(y[,j])-2)]; x2 <- x2[-(1:(max(y[,j])-2))] } }
new.zeta <- cbind(0,new.zeta);
}
## Replace vameans and lambda with elements in x
if(num.lv > 0) {
new.vameans[sel.i,] <- x[1:num.lv]; va.x <- x[-(1:num.lv)]
if(Lambda.struc == "diagonal") { new.lambda[sel.i,] <- va.x[1:(num.lv)]; va.x <- va.x[-(1:(num.lv))] }
if(Lambda.struc == "unstructured") { ## Rebuilt lambda array from unique elements
new.lambda[sel.i,,][lower.tri(new.lambda[sel.i,,],diag=TRUE)] <- va.x[1:(num.lv + num.lv * (num.lv - 1) / 2)]; va.x <- va.x[-(1:(num.lv + num.lv * (num.lv - 1) / 2))]
new.lambda[sel.i,,][upper.tri(new.lambda[sel.i,,])] <- new.lambda[sel.i,,][lower.tri(new.lambda[sel.i,,])]
}
}
grad.vameans <- NULL
grad.lambda <- matrix(0,num.lv,num.lv)
eta.mat <- matrix(new.beta0,n,p,byrow=TRUE) + new.vameans %*% t(new.theta) + offset
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + X %*% t(new.env)
if(!is.null(TR)) eta.mat <- eta.mat + matrix((Xd %*% B),n,p)
if(row.eff) eta.mat <- eta.mat + matrix(new.row.params,n,p)
if(family=="poisson") {
eta.mat <- eta.mat + calc.quad(new.lambda,new.theta,Lambda.struc)$mat
sum1 <- (y[sel.i,]-exp(eta.mat[sel.i,]))
for(l in 1:num.lv) { grad.vameans <- c(grad.vameans,sum(sum1 * new.theta[,l]) - new.vameans[sel.i,l]) }
}
if(family=="negative.binomial") {
eta.mat <- eta.mat - calc.quad(new.lambda,new.theta,Lambda.struc)$mat
sum1 <- -new.phi - (y[sel.i,] + new.phi) * (-new.phi / (exp(eta.mat[sel.i,]) + new.phi))
for(l in 1:num.lv) { grad.vameans <- c(grad.vameans,sum(sum1 * new.theta[,l]) - new.vameans[sel.i,l]) }
}
if(family=="binomial") {
sum1 <- dnorm(eta.mat[sel.i,]) * (y[sel.i,]-trial.size * pnorm(eta.mat[sel.i,])) / (pnorm(eta.mat[sel.i,]) * (1 - pnorm(eta.mat[sel.i,])) + 1e-10)
for(l in 1:num.lv) { grad.vameans <- c(grad.vameans,sum(sum1 * new.theta[,l]) - new.vameans[sel.i,l]) }
}
if(family=="ordinal") {
deriv.trunnorm <- matrix(NA,n,p)
for(j in 1:p) {
deriv.trunnorm[y[,j] == 1,j] <- -dnorm(new.zeta[j,1] - eta.mat[y[,j] == 1,j]) / pnorm(new.zeta[j,1] - eta.mat[y[,j] == 1,j])
deriv.trunnorm[y[,j] == max(y[,j]),j] <- dnorm(new.zeta[j,max(y[,j]) - 1] - eta.mat[y[,j] == max(y[,j]),j]) / (1 - pnorm(new.zeta[j,max(y[,j]) - 1] - eta.mat[y[,j] == max(y[,j]),j]))
if(max(y[,j])>2) {
j.levels <- 2:(max(y[,j])-1)
for(k in j.levels) { deriv.trunnorm[y[,j] == k,j] <- (-dnorm(new.zeta[j,k] - eta.mat[y[,j] == k,j]) + dnorm(new.zeta[j,k-1] - eta.mat[y[,j] == k,j])) / (pnorm(new.zeta[j,k] - eta.mat[y[,j] == k,j]) - pnorm(new.zeta[j,k - 1] - eta.mat[y[,j] == k,j])) }
}
}
deriv.trunnorm[!is.finite(deriv.trunnorm)] <- 0
for(l in 1:num.lv) { grad.vameans <- c(grad.vameans,sum(deriv.trunnorm[sel.i,] * new.theta[,l]) - new.vameans[sel.i,l]) }
}
if(family %in% c("poisson","negative.binomial")) {
if(family == "poisson") { deriv1 <- exp(eta.mat[sel.i,]) }
if(family == "negative.binomial") {
eta.mat <- matrix(new.beta0,n,p,byrow=TRUE) + offset + new.vameans %*% t(new.theta) - calc.quad(new.lambda,new.theta,Lambda.struc)$mat
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + X %*% t(new.env)
if(!is.null(TR)) eta.mat <- eta.mat + matrix((Xd %*% B),n,p)
if(row.eff) eta.mat <- eta.mat + matrix(new.row.params,n,p)
deriv1 <- (y[sel.i,] + new.phi) * (new.phi / (exp(eta.mat[sel.i,]) + new.phi))
}
if(Lambda.struc == "unstructured") {
theta2 <- sapply(1:p,function(j,theta) theta[j,]%*%t(theta[j,]), theta = new.theta)
theta2 <- t(theta2)
grad.lambda <- -0.5 * matrix(apply(deriv1 * theta2,2,sum),nrow=num.lv) + 0.5 * solve(new.lambda[sel.i,,]) - 0.5 * diag(nrow=num.lv)
}
if(Lambda.struc == "diagonal") {
theta2 <- new.theta^2
grad.lambda <- -0.5 * diag(apply(deriv1 * theta2,2,sum),nrow=num.lv) + 0.5 * solve(diag(x = new.lambda[sel.i,],nrow=num.lv)) - 0.5 * diag(nrow=num.lv)
}
}
if(family %in% c("binomial","ordinal")) {
if(Lambda.struc == "unstructured") {
theta2 <- sapply(1:p,function(j,theta) theta[j,] %*% t(theta[j,]), theta = new.theta)
theta2 <- t(theta2)
grad.lambda <- -0.5 * matrix(apply(theta2,2,sum),nrow=num.lv) + 0.5 * solve(new.lambda[sel.i,,]) - 0.5 * diag(nrow=num.lv)
}
if(Lambda.struc == "diagonal") {
theta2 <- new.theta^2
grad.lambda <- -0.5 * diag(apply(theta2,2,sum),nrow=num.lv) + 0.5 * solve(diag(x=new.lambda[sel.i,],nrow=num.lv)) - 0.5 * diag(nrow=num.lv)
}
}
grad.lambda2 <- NULL
if(Lambda.struc == "diagonal") grad.lambda2 <- diag(x=as.matrix(grad.lambda)) ## Only extract diagonal elements
if(Lambda.struc == "unstructured") grad.lambda2 <- grad.lambda[lower.tri(grad.lambda,diag=TRUE)]
return(c(grad.vameans,grad.lambda2))
}
## cross derivatives of VA logL - taking VA parameters vameans_sel.i and lambda_sel.i and diffing wrt to model parameters
grad.all.cross <- function(x,vameans=NULL,lambda=NULL,mod.x,sel.i) {
x2 <- mod.x
if(num.lv > 0) {
new.vameans <- vameans; new.lambda <- lambda
new.theta <- matrix(c(x2[1:(p * num.lv)]),p,num.lv); x2 <- x2[-(1:(p * num.lv))]
new.theta[upper.tri(new.theta)] <- 0
}
new.beta0 <- x2[1:p]; x2 <- x2[-(1:p)]
new.env <- NULL
if(!is.null(X)) {
if(is.null(TR)) { new.env <- matrix(x2[1:(p * num.X)],p,num.X); x2 <- x2[-(1:(p * num.X))]
} else { B=x2[1:nd]; x2 <- x2[-(1:nd)]}
}
new.row.params <- NULL; if(row.eff) { new.row.params <- x2[1:n]; x2 <- x2[-(1:n)] }
new.phi <- NULL; if(family == "negative.binomial") { new.phi <- x2[1:p]; x2 <- x2[-(1:p)] }
new.zeta <- NULL
if(family == "ordinal") {
new.zeta <- matrix(NA,p,max.levels - 2);
for(j in 1:p) { if(max(y[,j]) > 2) { new.zeta[j,1:(max(y[,j])-2)] <- x2[1:(max(y[,j])-2)]; x2 <- x2[-(1:(max(y[,j]) - 2))] } }
new.zeta <- cbind(0,new.zeta) }
## Replace vameans and lambda with elements in x
if(num.lv > 0) {
new.vameans[sel.i,] <- x[1:num.lv]; va.x <- x[-(1:num.lv)]
if(Lambda.struc == "diagonal") { new.lambda[sel.i,] <- va.x[1:(num.lv)]; va.x <- va.x[-(1:(num.lv))] }
if(Lambda.struc == "unstructured") { ## Rebuilt lambda array from unique elements
new.lambda[sel.i,,][lower.tri(new.lambda[sel.i,,],diag=TRUE)] <- va.x[1:(num.lv+num.lv * (num.lv - 1) / 2)]; va.x <- va.x[-(1:(num.lv + num.lv * (num.lv - 1) / 2))]
new.lambda[sel.i,,][upper.tri(new.lambda[sel.i,,])] <- new.lambda[sel.i,,][lower.tri(new.lambda[sel.i,,])]
}
}
eta.mat <- matrix(new.beta0,n,p,byrow=TRUE) + offset
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + X %*% t(new.env)
if(!is.null(TR)) eta.mat <- eta.mat + matrix((Xd %*% B),n,p)
if(row.eff) eta.mat <- eta.mat + matrix(new.row.params,n,p)
if(num.lv > 0) eta.mat <- eta.mat + new.vameans %*% t(new.theta)
grad.theta <- grad.beta0 <- grad.env <- grad.B <- grad.phi <- grad.row.params <- grad.zeta <- NULL
if(row.eff) grad.row.params <- numeric(n)
if(Lambda.struc == "unstructured" & num.lv > 0) {
Lambda.theta <- array(NA,dim=c(p,n,num.lv))
for(i2 in 1:n) { Lambda.theta[,i2,] <- new.theta %*% new.lambda[i2,,] }
}
if(family=="poisson") {
if(num.lv > 0) eta.mat <- eta.mat + calc.quad(new.lambda,new.theta,Lambda.struc)$mat
grad.beta0 <- colSums(y-exp(eta.mat))
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(y,1,new.vameans[,l],"*") - sweep(t(Lambda.theta[,,l]),1,new.vameans[,l],"+") * exp(eta.mat)
if(Lambda.struc == "diagonal") sum1 <- sweep(y,1,new.vameans[,l],"*") - sweep((new.lambda[,l]%*%t(new.theta[,l])),1,new.vameans[,l],"+")*exp(eta.mat)
grad.theta <- c(grad.theta,colSums(sum1))
}
}
if(!is.null(X) && is.null(TR)) {
for (l in 1:num.X) {
sum1 <- sweep(y-exp(eta.mat),1,X[,l],"*")
grad.env <- c(grad.env,colSums(sum1)) # else{grad.env <- c(grad.env,sum(sum1)) }
}
}
if(!is.null(TR)) {
sum1 <- c(y-exp(eta.mat)) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) { grad.row.params <- rowSums(y - exp(eta.mat)) }
}
if(family=="negative.binomial") {
eta.mat <- matrix(new.beta0,n,p,byrow=TRUE) + offset
if(!is.null(X) && is.null(TR)) eta.mat <- eta.mat + X %*% t(new.env)
if(!is.null(TR)) eta.mat <- eta.mat + matrix((Xd %*% B),n,p)
if(row.eff) eta.mat <- eta.mat + matrix(new.row.params,n,p)
if(num.lv > 0) eta.mat <- eta.mat + new.vameans %*% t(new.theta)
phi.mat <- matrix(new.phi,n,p,byrow=TRUE)
mu.mat.noquad <- eta.mat
if(num.lv > 0) eta.mat <- eta.mat - calc.quad(new.lambda,new.theta,Lambda.struc)$mat
grad.beta0 <- colSums(-phi.mat - (y + phi.mat) * (-phi.mat / (exp(eta.mat) + phi.mat)))
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(-phi.mat,1,new.vameans[,l],"*") + sweep(t(Lambda.theta[,,l]),1,-new.vameans[,l],"+") * (-(y + phi.mat) * (phi.mat / (exp(eta.mat) + phi.mat)))
if(Lambda.struc == "diagonal") sum1 <- sweep(-phi.mat,1,new.vameans[,l],"*") + sweep((new.lambda[,l] %*% t(new.theta[,l])),1,-new.vameans[,l],"+") * (-(y + phi.mat) * (phi.mat / (exp(eta.mat) + phi.mat)))
grad.theta <- c(grad.theta,colSums(sum1))
}
}
if(!is.null(X) && is.null(TR)) {
for(l in 1:num.X) {
sum1 <- sweep(-phi.mat - (y + phi.mat) * (-phi.mat / (exp(eta.mat) + phi.mat)),1,X[,l],"*")
grad.env <- c(grad.env,colSums(sum1))
}
}
if(!is.null(TR)) {
sum1 <- c(-phi.mat - (y + phi.mat) * (-phi.mat / (exp(eta.mat) + phi.mat))) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) grad.row.params <- rowSums(-phi.mat - (y + phi.mat) * (-phi.mat/(exp(eta.mat) + phi.mat)))
grad.phi <- colSums(-mu.mat.noquad + 1 + log(phi.mat) - digamma(phi.mat) - log(1 + phi.mat / exp(eta.mat)) - (y + phi.mat)/(phi.mat + exp(eta.mat)) + digamma(y + phi.mat))
}
if(family=="binomial") {
grad.beta0 <- colSums(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat))),na.rm=TRUE)
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat))/(pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10),1,new.vameans[,l],"*") - t(Lambda.theta[,,l])
if(Lambda.struc == "diagonal") sum1 <- sweep(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10),1,new.vameans[,l],"*") - (new.lambda[,l] %*% t(new.theta[,l]))
grad.theta <- c(grad.theta,colSums(sum1,na.rm=TRUE)) }
}
if(!is.null(X) && is.null(TR)) {
for(l in 1:num.X) {
sum1 <- sweep(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10),1,X[,l],"*")
grad.env <- c(grad.env,colSums(sum1))
}
}
if(!is.null(TR)) {
sum1 <- c(dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)) + 1e-10)) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) { grad.row.params <- rowSums((dnorm(eta.mat) * (y - trial.size * pnorm(eta.mat)) / (pnorm(eta.mat) * (1 - pnorm(eta.mat)))),na.rm=TRUE) }
}
if(family=="ordinal") {
grad.zeta <- matrix(0,p,ncol(new.zeta))
deriv.trunnorm <- matrix(0,n,p)
for(j in 1:p) {
deriv.trunnorm[y[,j] == 1,j] <- -dnorm(new.zeta[j,1] - eta.mat[y[,j] == 1,j]) / pnorm(new.zeta[j,1] - eta.mat[y[,j] == 1,j])
deriv.trunnorm[y[,j] == max(y[,j]),j] <- dnorm(new.zeta[j,max(y[,j]) - 1] - eta.mat[y[,j] == max(y[,j]),j]) / (1 - pnorm(new.zeta[j,max(y[,j]) - 1] - eta.mat[y[,j] == max(y[,j]),j]))
j.levels <- (1:max(y[,j]))[-c(1,max(y[,j]))]
for(k in j.levels) {
deriv.trunnorm[y[,j] == k,j] <- (-dnorm(new.zeta[j,k] - eta.mat[y[,j] == k,j]) + dnorm(new.zeta[j,k - 1] - eta.mat[y[,j] == k,j])) / (pnorm(new.zeta[j,k] - eta.mat[y[,j] == k,j]) - pnorm(new.zeta[j,k-1] - eta.mat[y[,j] == k,j])) }
}
deriv.trunnorm[!is.finite(deriv.trunnorm)] <- 0
grad.beta0 <- colSums(deriv.trunnorm,na.rm=TRUE)
if(num.lv > 0) {
for(l in 1:num.lv) {
if(Lambda.struc == "unstructured") sum1 <- sweep(deriv.trunnorm,1,new.vameans[,l],"*") - t(Lambda.theta[,,l])
if(Lambda.struc == "diagonal") sum1 <- sweep(deriv.trunnorm,1,new.vameans[,l],"*") - (new.lambda[,l] %*% t(new.theta[,l]))
grad.theta <- c(grad.theta,colSums(sum1,na.rm=TRUE))
}
}
if(!is.null(X) && is.null(TR)) {
for(l in 1:num.X) { sum1 <- sweep(deriv.trunnorm,1,X[,l],"*");
grad.env <- c(grad.env,colSums(sum1,na.rm=TRUE))
}
}
if(!is.null(TR)) {
sum1 <- c(deriv.trunnorm) * Xd
grad.B <- c(grad.B,colSums(sum1))
}
if(row.eff) { grad.row.params <- rowSums(deriv.trunnorm,na.rm=TRUE) }
for(j in 1:p) {
zeta0 <- new.zeta[j,]
grad.zeta[j,max(y[,j]) - 1] <- grad.zeta[j,max(y[,j]) - 1] - sum(dnorm(zeta0[max(y[,j]) - 1] - eta.mat[which(y[,j] == max(y[,j])) ,j]) / (1 - pnorm(zeta0[max(y[,j]) - 1] - eta.mat[which(y[,j] == max(y[,j])),j])))
j.levels <- (1:max(y[,j]))[-c(1,max(y[,j]))]
for(k in j.levels) {
grad.zeta[j,k] <- grad.zeta[j,k] + sum(dnorm(zeta0[k] - eta.mat[y[,j] == k,j]) / (pnorm(zeta0[k] - eta.mat[y[,j] == k,j]) - pnorm(zeta0[k - 1] - eta.mat[y[,j] == k,j])))
grad.zeta[j,k - 1] <- grad.zeta[j,k - 1] - sum(dnorm(zeta0[k - 1] - eta.mat[y[,j] == k,j]) / (pnorm(zeta0[k] - eta.mat[y[,j] == k,j]) - pnorm(zeta0[k - 1] - eta.mat[y[,j] == k,j])))
}
}
grad.zeta <- grad.zeta[,-1]
if(length(unique(apply(y,2,max))) > 1) { grad.zeta <- c(grad.zeta[-which(grad.zeta == 0)]) }
}
if(num.lv>0){
grad.theta=matrix(grad.theta,nrow=p,ncol=num.lv)
grad.theta[upper.tri(grad.theta)] <- 0}
if(row.eff) grad.row.params[1]=0;
return(c(grad.theta, grad.beta0, grad.env, grad.B, grad.row.params, grad.phi, grad.zeta))
}
zero.cons <- which(theta==0)
if(family == "ordinal") {
zeta2 <- as.vector(t(zeta[,-1]))
find.na.zeta <- which(is.na(zeta2)); if(length(find.na.zeta) > 0) zeta2 <- zeta2[-find.na.zeta]
}
if(family != "ordinal") { zeta2 <- NULL }
A.mat <- -nd2(x0=c(theta,beta0,env,B,row.params,phi,zeta2), f = grad.all.mod, vameans = vameans, lambda = lambda)
r1.off=NULL; if(row.eff) {r1.off=length(c(theta,beta0,env,B))+1}
if(length(zero.cons) > 0) { A.mat <- A.mat[-c(zero.cons,r1.off),]; A.mat <- A.mat[,-c(zero.cons,r1.off)]
} else if(row.eff) {A.mat <- A.mat[-r1.off,]; A.mat <- A.mat[,-r1.off] } #!!!!!
w1=FALSE
if(num.lv > 0) {
D.mat.fun <- function(i3) {
if(Lambda.struc == "unstructured") { lambda.i <- lambda[i3,,][lower.tri(lambda[i3,,],diag=TRUE)] }
if(Lambda.struc == "diagonal") { lambda.i <- lambda[i3,]; }
return(-nd2(f = grad.all.var, x0=c(vameans[i3,],lambda.i), vameans = vameans, lambda = lambda, mod.x = c(theta,beta0,env,B,row.params,phi,zeta2), sel.i = i3))
}
D.mat <- vector("list",n);
if(p>1){
unique.ind <- which(!duplicated(y))
for(k in 1:length(unique.ind)) {
subD <- D.mat.fun(i3=unique.ind[k]); match.seq <- which(apply(y,1,identical,y[unique.ind[k],]) == 1) ## Unfortunately replace() only works if you want to replace elements with numerics
if( any(is.nan(subD)) ){ logi<-is.nan(subD); subD[logi]=diag(dim(subD)[1])[logi]; w1=TRUE}
for(k2 in match.seq) { D.mat[[k2]] <- subD }
}
} else {
unique.ind=1:n
for(k in 1:length(unique.ind)) {
subD <- D.mat.fun(i3=unique.ind[k]);
if( any(is.nan(subD)) ){ logi<-is.nan(subD); subD[logi]=diag(dim(subD)[1])[logi]; w1=TRUE}
for(k2 in 1:n) { D.mat[[k2]] <- subD }
}
}
D.mat <- as.matrix(Matrix::bdiag(D.mat)); rm(subD)
B.mat <- vector("list",n)
for(i3 in 1:length(unique.ind)) {
#cat("Onto row",i3,"\n")
if(Lambda.struc == "unstructured") { lambda.i <- lambda[unique.ind[i3],,][lower.tri(lambda[unique.ind[i3],,],diag=TRUE)] }
if(Lambda.struc == "diagonal") { lambda.i <- lambda[unique.ind[i3],]; }
subB.mat <- -nd2(f = grad.all.cross, x0=c(vameans[unique.ind[i3],],lambda.i), vameans = vameans, lambda = lambda, mod.x = c(theta,beta0,env,B,row.params,phi,zeta2), sel.i = unique.ind[i3])
if(length(zero.cons) > 0) { subB.mat <- subB.mat[-c(zero.cons,r1.off),]
}else if(row.eff) {subB.mat <- subB.mat[-c(r1.off),] }
if( any(is.nan(subB.mat)) ){ logi<-is.nan(subB.mat); subB.mat[logi]=0; w1=TRUE}
if(p>1){
match.seq <- which(apply(y,1,identical,y[unique.ind[i3],]) == 1)
for(k2 in match.seq) { B.mat[[k2]] <- subB.mat }
} else { for(k2 in 1:n) { B.mat[[k2]] <- subB.mat }}
}
B.mat <- do.call(cbind,B.mat)
cov.mat.mod <- MASS::ginv(A.mat-B.mat%*%solve(D.mat)%*%t(B.mat))
}
if(num.lv == 0) cov.mat.mod <- MASS::ginv(A.mat)
if(w1) warning("Standard errors of parameters might not be reliable, due to the technical remedy.\n");
sd.errors <- sqrt(diag(abs(cov.mat.mod)))
theta.errors <- NULL;
sdout=NULL
if(num.lv > 0) {
theta.vec <- sd.errors[1:(p*num.lv-length(zero.cons))]; theta.errors <- matrix(0,p,num.lv);
if(length(zero.cons) > 0) theta.errors[-zero.cons] <- theta.vec
if(length(zero.cons) == 0) theta.errors <- matrix(theta.vec,p,num.lv)
sd.errors <- sd.errors[-(1:(p*num.lv-length(zero.cons)))]
if(num.lv > 1) theta.errors[upper.tri(theta.errors)] <- 0
rownames(theta.errors) <- colnames(y);
colnames(theta.errors) <- paste("LV", 1:num.lv, sep="");
sdout$theta <- theta.errors
}
beta.errors <- sd.errors[1:p]; sd.errors <- sd.errors[-(1:p)]
names(beta.errors) <- colnames(y)
sdout$beta0 <- beta.errors
env.errors=NULL
if(!is.null(X)) {
if(is.null(TR)){
beta.errors <- matrix(sd.errors[1:(p*ncol(X))],p,ncol(X));
rownames(beta.errors) <- colnames(y); colnames(beta.errors) <- colnames(X);
sd.errors <- sd.errors[-(1:(p*ncol(X)))]
sdout$Xcoef <- beta.errors
} else {
coef.errors <- sd.errors[1:nd]
names(coef.errors) <- colnames(Xd)
sd.errors <- sd.errors[-(1:(nd))]
sdout$B <- coef.errors
}
}
row.params.errors <- NULL;
if(row.eff) {
row.params.errors <- c(0,sd.errors[1:(n-1)]); sd.errors <- sd.errors[-(1:(n-1))];
names(row.params.errors) <- rownames(y)
sdout$row.params <- row.params.errors
}
phi.errors <- NULL;
if(family == "negative.binomial") {
phi.errors <- sd.errors[1:p]; sd.errors <- sd.errors[-(1:p)];
names(phi.errors) <- colnames(y)
sdout$inv.phi <- phi.errors
}
zeta.errors <- NULL;
if(family == "ordinal") {
zeta.errors <- matrix(NA,p,max.levels - 2)
for(j in 1:p) { if(max(y[,j]) > 2) { zeta.errors[j,1:(max(y[,j]) - 2)] <- sd.errors[1:(max(y[,j]) - 2)]; sd.errors <- sd.errors[-(1:(max(y[,j]) - 2))] } }
zeta.errors <- cbind(0,zeta.errors)
rownames(zeta.errors) <- colnames(y); colnames(zeta.errors) <- paste(1:(max(y)-1),"|",2:max(y),sep="")
sdout$zeta <- zeta.errors
}
return(sdout)
}
## A function to compute highly accurate first-order derivatives
## From Fornberg and Sloan (Acta Numerica, 1994, p. 203-267; Table 1, page 213)
## Adapted by Scott Foster from code nicked off the net 2007
nd2 <- function(x0, f, m = NULL, D.accur = 2, ...) {
D.n <- length(x0)
if(is.null(m)) { D.f0 <- f(x0, ...); m <- length(D.f0) }
if(D.accur == 2) { D.w <- tcrossprod(rep(1,m), c(-1/2,1/2)); D.co <- c(-1, 1) }
else {
D.w <- tcrossprod(rep(1,m), c(1/12, -2/3, 2/3, -1/12))
D.co <- c(-2, -1, 1, 2)
}
D.n.c <- length(D.co)
macheps <- .Machine$double.eps
D.h <- macheps^(1/3) * abs(x0)
D.deriv <- matrix(NA, nrow = m, ncol = D.n)
for (ii in 1:D.n) {
D.temp.f <- matrix(0, m, D.n.c)
for (jj in 1:D.n.c) {
D.xd <- x0 + D.h[ii] * D.co[jj] * (1:D.n == ii)
D.temp.f[,jj] <- f( D.xd, ...)
}
D.deriv[,ii] <- rowSums(D.w * D.temp.f) / D.h[ii]
}
return(D.deriv)
}
calc.quad <- function(lambda,theta,Lambda.struc) {
if(Lambda.struc == "diagonal") out <- 0.5 * (lambda) %*% t(theta^2)
if(Lambda.struc == "unstructured") {
if(inherits(lambda,"array")) { n <- dim(lambda)[1]; num.lv <- dim(lambda)[2] }
if(inherits(lambda,"matrix")) { num.lv <- dim(lambda)[2]; n <- 1 }
if(inherits(theta,"matrix")) { p <- dim(theta)[1] }
if(inherits(theta,"numeric")) { p <- 1; theta <- matrix(theta,1) }
out <- matrix(NA,n,p)
if(n == 1) {
for(j in 1:p) { out[1,j] <- 0.5 * t(theta[j,]) %*% lambda %*% theta[j,] }
}
if(n > 1) {
if(n <= p) out <- t(sapply(1:n, function(x) 0.5 * rowSums(theta * (theta %*% lambda[x,,]))))
if(n > p) {
lambda.mat <- aperm(lambda,c(3,2,1)); dim(lambda.mat) <- c(num.lv,num.lv * n)
f <- function(x) 0.5 * rowSums(matrix((t(lambda.mat) * theta[x,]) %*% theta[x,],ncol=num.lv,byrow=TRUE))
out <- sapply(1:p,f)
}
}
}
return(list(mat = out, mat.sum = sum(out)))
}
lambda.convert <- function(lambda,type=1) {
if(type == 1) { ## Convert unstructured to diagonal
n <- dim(lambda)[1]
lambda2 <- matrix(1,dim(lambda)[1],dim(lambda)[3])
for(i in 1:n) { lambda2[i,] <- diag(lambda[i,,]) }
return(lambda2)
}
if(type == 2) { ## Convert diagonal to unstructured
n <- nrow(lambda)
lambda2 <- array(0,dim=c(nrow(lambda),ncol(lambda),ncol(lambda)))
for(i in 1:n) { diag(lambda2[i,,]) <- lambda[i,] }
return(lambda2)
}
}
inf.criteria <- function(fit)
{
family=fit$family
abund=fit$y
num.lv=fit$num.lv
n <- dim(abund)[1]
p <- dim(abund)[2]
k<-attributes(logLik.gllvm(fit))$df
BIC <- -2*fit$logL + (k) * log(n*p)
# AIC
AIC <- -2*fit$logL + (k) * 2
# AICc
AICc <- AIC + 2*k*(k+1)/(n*p-k-1)
list(BIC = BIC, AIC = AIC, AICc = AICc, k = k)
}
# Creates matrix of fourth corner terms from a vector
getFourthCorner<- function(object){
if(is.null(object$X) || is.null(object$TR)) stop();
n1=colnames(object$X)
n2=colnames(object$TR)
nams=names(object$params$B)
fx<-cbind(apply(sapply(n1,function(x){grepl(x, nams)}),1,any), apply(sapply(n2,function(x){grepl(x, nams)}),1,any))
fourth.index<-rowSums(fx)>1
nams2=nams[fourth.index]
fourth.corner=object$params$B[fourth.index]
splitsnam <- sapply(nams2,strsplit, split = ":")
isinvec <- function(vec, ob =""){
ob %in% vec
}
n1 <- n1[(n1 %in% unlist(splitsnam))]
n2 <- n2[(n2 %in% unlist(splitsnam))]
i=1; j=1;
fourth<-matrix(0,length(n1),length(n2))
for (i in 1:length(n1)) {
for (j in 1:length(n2)) {
fur=(unlist(lapply(splitsnam, isinvec, ob = n1[i])) + unlist(lapply(splitsnam, isinvec, ob = n2[j]))) >1
if(any(fur)){ fourth[i,j]=fourth.corner[fur]}
}
}
colnames(fourth)=n2
rownames(fourth)=n1
return(fourth)
}
# Calculates standard errors for random effects
sdrandom<-function(obj, Vtheta, incl, ignore.u = FALSE,return.covb = FALSE, type = NULL){
#For num.RR we treat the LV as zero
if(!is.null(type)){
if(type=="marginal"){
ignore.u <- TRUE
}
}
r <- obj$env$random
par = obj$env$last.par.best
num.lv <- obj$env$data$num_lv
num.lv.c <- obj$env$data$num_lv_c
num.RR <- obj$env$data$num_RR
nu<-nrow(obj$env$data$y)
xb<- obj$env$data$xb
random <- obj$env$data$random
radidx <- sum(names(obj$env$last.par.best[r])%in%c("r0","Br","b_lv"))
if(is.null(type)){
if((num.lv.c+num.RR)==0){
type <- "residual"
}else{
type <- "conditional"
}
}
if((num.lv+num.lv.c+radidx)>0|num.RR>0&random[3]>0){
hessian.random <- obj$env$spHess(par, random = TRUE)
L <- obj$env$L.created.by.newton
}
n <- nrow(obj$env$data$y)
p <- ncol(obj$env$data$y)
lv.X <- obj$env$data$x_lv
if((num.lv+num.lv.c)>0){
sigma.lv <- abs(obj$par[names(obj$par)=="sigmaLV"])
}
diag.cov.random <- array(0,dim=c(n,num.lv.c+num.lv+num.RR,num.lv.c+num.lv+num.RR))
if (ignore.u&num.RR>0&random[3]==0) {
diag.term2 <- 0
Q <- matrix(0,nrow=(num.lv+num.lv.c+num.RR)*n+radidx,ncol=dim(Vtheta)[1])
if((num.lv.c+num.lv+radidx)==0)A<-Q
if((num.lv.c+num.RR)>0){
for(q in 1:(num.lv.c+num.RR)){
Q[(1:n)+n*(q-1)+radidx,which(names(obj$par[incl])=="b_lv")[(1:(ncol(lv.X)-(q-1)))+(q-1)*ncol(lv.X)-(q-1)*(q-2)/2]] <- lv.X[,1:ncol(lv.X)-(q-1),drop=F]#/fit$params$sigma.lv[q] #divide here to multiply later in ordiplot
}
}
} else {
if((num.lv.c+num.lv+radidx)>0){
f <- obj$env$f
w <- rep(0, length(par))
reverse.sweep <- function(i) {
w[i] <- 1
f(par, order = 1, type = "ADGrad", rangeweight = w, doforward = 0)[r]
}
nonr <- setdiff(seq_along(par), r)
tmp <- sapply(nonr, reverse.sweep)
if (!is.matrix(tmp))
tmp <- matrix(tmp, ncol = length(nonr))
A <- solve(hessian.random, tmp[, incl]) #n*d by N_psi
row.names(A) <- names(obj$env$last.par.best[r])
}
if(radidx==0){
if((num.lv.c)>0&num.RR>0){
A <- rbind(A[1:(num.lv.c*nu),],matrix(0,nrow=num.RR*nu,ncol=ncol(A)),A[-c(1:(num.lv.c*nu)),])
}else if(num.lv>0&num.RR>0){
A <- rbind(matrix(0,nrow=num.RR*nu,ncol=ncol(A)),A)
}
}else{
if(num.RR>0 & random[3]==0)A <- rbind(A[1:(num.lv.c*nu+radidx),],matrix(0,nrow=num.RR*nu,ncol=ncol(A)),A[-c(1:(num.lv.c*nu+radidx)),])
}
#we never scale the prediction regions for num.lv by sigma.lv.
if(num.lv.c>0&num.lv>0){
sigma.lv <- c(sigma.lv[1:num.lv.c],rep(1,num.RR*(1-random[3])),rep(1,num.lv))
}else if(num.lv>0&num.lv.c==0){
sigma.lv <- c(rep(1,num.RR*(1-random[3])),rep(1,num.lv))
}else if(num.lv.c>0&num.lv==0){
sigma.lv <- c(sigma.lv,rep(1,num.RR*(1-random[3])))
}else if(num.RR>0){
sigma.lv <- rep(1,num.RR*(1-random[3]))
}
if(num.RR>0&random[3]==0){
row.names(A)[row.names(A)==""] <- rep("XB",num.RR*nu)
}
#Matrix Q for predictors
Q <- matrix(0,nrow=(num.lv+num.lv.c+num.RR*ifelse(random[3]==0,1,0))*n+radidx,ncol=dim(Vtheta)[1])
if((num.lv.c+num.lv+radidx)==0)A<-Q
if((num.lv.c+num.RR)>0&random[3]==0){
for(q in 1:(num.lv.c+num.RR)){
Q[(1:n)+n*(q-1)+radidx,which(names(obj$par[incl])=="b_lv")[(1:(ncol(lv.X)-(q-1)))+(q-1)*ncol(lv.X)-(q-1)*(q-2)/2]] <- lv.X[,1:ncol(lv.X)-(q-1),drop=F]#/fit$params$sigma.lv[q] #divide here to multiply later in ordiplot
}
}
# if(is.null(type)&(num.lv.c+num.RR)==0){
# type <- "residual"
# }else{
# type <- "conditional"
# }
# if((num.lv+num.lv.c+radidx)>0){
if((num.lv+num.lv.c)>0|any(random>0)){
if(type=="conditional"|type=="marginal"&random[3]>0){
S <- diag(c(rep(1,radidx),rep(sigma.lv,each=nu)))
diag.term2 <- (Q+S%*%A)%*%Vtheta%*%t(Q+S%*%A)
colnames(diag.term2)<-row.names(diag.term2)<-row.names(A)
}else if(type=="residual"){
diag.term2 <- (A)%*%Vtheta%*%t(A)
colnames(diag.term2)<-row.names(diag.term2)<-row.names(A)
}else if(type=="marginal"&random[3]==0){
diag.term2 <- Q%*%Vtheta%*%t(Q)
}
diag.term2 <- as.matrix(diag.term2)
}
}
if((num.lv+num.lv.c+radidx)>0){
diag.term1 <- Matrix::chol2inv(L)
if(radidx>0&num.RR>0&random[3]==0){
diag.term1 <- rbind(diag.term1[1:(num.lv.c*n+radidx),],matrix(0,nrow=num.RR*n,ncol=ncol(diag.term1)),diag.term1[-c(1:(num.lv.c*n+radidx)),])
diag.term1 <- cbind(diag.term1[,1:(num.lv.c*n+radidx)],matrix(0,nrow=nrow(diag.term1),ncol=num.RR*n),diag.term1[,-c(1:(num.lv.c*n+radidx))])
}else if(num.RR>0&random[3]==0){
if(num.lv.c>0){
diag.term1 <- rbind(diag.term1[1:(num.lv.c*n),],matrix(0,nrow=num.RR*n,ncol=ncol(diag.term1)),diag.term1[-c(1:(num.lv.c*n)),])
diag.term1 <- cbind(diag.term1[,1:(num.lv.c*n)],matrix(0,nrow=nrow(diag.term1),ncol=num.RR*n),diag.term1[,-c(1:(num.lv.c*n))])
}else if(num.lv>0){
diag.term1 <- rbind(matrix(0,nrow=num.RR*n,ncol=ncol(diag.term1)),diag.term1)
diag.term1 <- cbind(matrix(0,nrow=nrow(diag.term1),ncol=num.RR*n),diag.term1)
}
}
diag.term1 <- as.matrix(diag.term1)
}else if((num.lv.c+num.lv+radidx)==0){
diag.term1<-0
}
if((num.lv+num.lv.c+radidx)>0&type!="marginal"){
diag.term1 <- Matrix::chol2inv(L)
if(radidx>0 & num.RR>0 && random[3] == 0){
diag.term1 <- rbind(diag.term1[1:(num.lv.c*nu+radidx),],matrix(0,nrow=num.RR*nu,ncol=ncol(diag.term1)),diag.term1[-c(1:(num.lv.c*nu+radidx)),])
diag.term1 <- cbind(diag.term1[,1:(num.lv.c*nu+radidx)],matrix(0,nrow=nrow(diag.term1),ncol=num.RR*nu),diag.term1[,-c(1:(num.lv.c*nu+radidx))])
}else if(num.RR>0 && random[3] == 0){
if(num.lv.c>0){
diag.term1 <- rbind(diag.term1[1:(num.lv.c*nu),],matrix(0,nrow=num.RR*nu,ncol=ncol(diag.term1)),diag.term1[-c(1:(num.lv.c*nu)),])
diag.term1 <- cbind(diag.term1[,1:(num.lv.c*nu)],matrix(0,nrow=nrow(diag.term1),ncol=num.RR*nu),diag.term1[,-c(1:(num.lv.c*nu))])
}else if(num.lv>0){
diag.term1 <- rbind(matrix(0,nrow=num.RR*nu,ncol=ncol(diag.term1)),diag.term1)
diag.term1 <- cbind(matrix(0,nrow=nrow(diag.term1),ncol=num.RR*nu),diag.term1)
}
}
diag.term1 <- as.matrix(diag.term1)
}else if(type=="marginal"|(num.lv.c+num.lv)==0){
diag.term2 <- Q%*%Vtheta%*%t(Q)
diag.term1<-0
colnames(diag.term2)<-rep("XB",n*(num.RR+num.lv.c))
}
if(type%in%c("conditional")){
S <- diag(c(rep(1,radidx),rep(sigma.lv,each=nu)))
diag.term1 <- S%*%diag.term1%*%S
colnames(diag.term2)<-row.names(diag.term2)<-row.names(A)
}
covb <- diag.term1 + diag.term2
row.names(covb)<-colnames(covb) <- colnames(diag.term2)
out <- list()
#separate errors row-effects
ser0 <- diag(covb[colnames(covb)=="r0",colnames(covb)=="r0"])
covb <- covb[colnames(covb)!="r0",colnames(covb)!="r0"]
if(random[1] >0) {
CovRow <- matrix(ser0);
out$row <- CovRow
}
seb_lv <- diag(covb[colnames(covb)=="b_lv",colnames(covb)=="b_lv"])
if(random[3]>0){
seb_lv <- diag(covb[colnames(covb)=="b_lv",colnames(covb)=="b_lv"])
covb_lvErr <- matrix(seb_lv,ncol=num.lv.c+num.RR)
out$Ab_lv <- covb_lvErr
covsB <- as.matrix(covb[colnames(covb)=="b_lv",colnames(covb)=="b_lv"])
}
if(!return.covb)covb <- covb[colnames(covb)!="b_lv",colnames(covb)!="b_lv"]
#separate errors AB
seBr <- diag(covb[colnames(covb)=="Br",colnames(covb)=="Br"])
covb <- covb[colnames(covb)!="Br",colnames(covb)!="Br"]
if(random[2]>0){
CovABerr<-array(0, c(p,ncol(xb),ncol(xb)))
for (i in 1:ncol(xb)) {
CovABerr[,i,i] <- seBr[1:p]; seBr<-seBr[-(1:p)]
}
out$Ab <- CovABerr
}
# if((num.RR+num.lv+num.lv.c)>0){
if(!return.covb){
covb <- as.matrix(covb)
# se <- simplify2array(sapply(1:nu,function(i)covb[seq(i,nu*(num.lv+num.lv.c+num.RR),by=nu),seq(i,nu*(num.lv+num.lv.c+num.RR),by=nu)],simplify=F))
#re-order, select submatrices
try({
if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(random[3]==0&type!="residual",1,0))>0)
se <- simplify2array(sapply(1:nu,function(i)covb[seq(i,nu*(num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(random[3]==0&type!="residual",1,0)),by=nu),seq(i,nu*(num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(random[3]==0&type!="residual",1,0)),by=nu)],simplify=F))
if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(random[3]==0&type!="residual",1,0))>1){
se <- aperm(se,c(3,2,1))
}else if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(random[3]==0&type!="residual",1,0))>0){
se <- as.matrix(se)
}
#add error for Bs if random
if((num.RR+num.lv.c)>0&type!="residual"&random[3]>0){
if((num.lv+num.lv.c+num.RR)>0){
se2 <- array(0,dim=c(n,num.lv.c+num.RR+num.lv*ifelse(type=="marginal",0,1),num.lv.c+num.RR+num.lv*ifelse(type=="marginal",0,1)))
if(type=="conditional"&(num.lv+num.lv.c)>0&num.RR>0){
se2[,-c((num.lv.c+1):(num.lv.c+num.RR)),-c((num.lv.c+1):(num.lv.c+num.RR))] <- se #0s for RRR
}else if(type=="conditional"&(num.lv+num.lv.c)>0){
se2 <- se
}
se <- se2
}
for(i in 1:n){
Q <- as.matrix(Matrix::bdiag(replicate(num.RR+num.lv.c,lv.X[i,,drop=F],simplify=F)))
temp <- Q%*%covsB%*%t(Q)
temp[col(temp)!=row(temp)] <- 2*temp[col(temp)!=row(temp)] ##should be double the covariance
se[i,1:(num.RR+num.lv.c),1:(num.RR+num.lv.c)] <- se[i,1:(num.RR+num.lv.c),1:(num.RR+num.lv.c)] + temp
}
}
if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(type!="residual",1,0)) > 0) {
out$A <- se
}
},silent=T)
}else{
out <- covb
}
# }
return(out)
}
start.values.randomX <- function(y, X, family, formula =NULL, starting.val, Power = NULL, link=NULL, method="VA", Ntrials = 1) {
y <- as.matrix(y)
Xb <- as.matrix(model.matrix(formula, data = data.frame(X)))
rnam <- colnames(Xb)[!(colnames(Xb) %in% c("(Intercept)"))]
Xb <- as.matrix(Xb[, rnam]); #as.matrix(X.new[, rnam])
if(NCOL(Xb) == 1) colnames(Xb) <- rnam
Xb <- as.matrix(Xb)
n <- NROW(y)
p <- NCOL(y)
tr0 <- try({
if(starting.val %in% c("res", "random")){
if(family %in% c("poisson", "negative.binomial", "binomial", "ZIP")){
if(family == "ZIP") family <- "poisson"
f1 <- gllvm.TMB(y=y, X=X, family = family, formula=formula, num.lv=0, starting.val = "zero", link =link, Ntrials = Ntrials) #, method=method
coefs0 <- as.matrix(scale((f1$params$Xcoef), scale = FALSE))
Br <- coefs0/max(apply(coefs0, 2, sd))
sigmaB <- cov(Br)
Br <- t(Br)
} else if(family == "tweedie"){
coefs0 <- NULL
for(j in 1:p){
fitj <- try(glm( y[,j] ~ Xb, family = statmod::tweedie(var.power=Power, link.power=0) ),silent = TRUE)
if(!inherits(fitj, "try-error")){
coefs0 <- rbind(coefs0, fitj$coefficients[-1])
} else { coefs0 <- rbind(coefs0,rnorm(dim(Xb)[2])); }
}
Br <- coefs0/max(apply(coefs0, 2, sd))
sigmaB <- cov(Br)
Br <- t(Br)
} else {
Br <- matrix(0, ncol(Xb), p)
sigmaB <- diag(ncol(Xb))
}
} else {
Br <- matrix(0, ncol(Xb), p)
sigmaB <- diag(ncol(Xb))
}
}, silent = TRUE)
if(inherits(tr0, "try-error")){
Br <- matrix(0, ncol(Xb), p)
sigmaB <- diag(ncol(Xb))
}
return(list(Br = Br, sigmaB = sigmaB))
}
# Calculates adjusted prediction errors for random effects
# sdA<-function(fit){
# n<-nrow(fit$y)
# A<- fit$Hess$cov.mat.mod #
# B<- fit$Hess$Hess.full[fit$Hess$incl, fit$Hess$incla]
# C<- fit$Hess$Hess.full[fit$Hess$incla, fit$Hess$incl]
# D<- solve(fit$Hess$Hess.full[fit$Hess$incla, fit$Hess$incla])
# covb <- (D%*%C)%*%(A)%*%(B%*%t(D))
# se <- (diag(abs(covb)))
# CovAerr<-array(0, dim(fit$A))
# for (i in 1:ncol(fit$A)) {
# CovAerr[,i,i] <- se[1:n]; se<-se[-(1:n)]
# }
# CovAerr
# }
#
# sdB<-function(fit){
# n<-nrow(fit$y)
# p<-ncol(fit$y)
# incla<-rep(FALSE, length(fit$Hess$incl))
# incla[names(fit$TMBfn$par)%in%c("u", "Br")] <- TRUE
# if(fit$row.eff == "random") incla[names(fit$TMBfn$par)%in%c("r0")] <- TRUE
# fit$Hess$incla <- incla
#
# A<- fit$Hess$cov.mat.mod #
# B<- fit$Hess$Hess.full[fit$Hess$incl, fit$Hess$incla]
# C<- fit$Hess$Hess.full[fit$Hess$incla, fit$Hess$incl]
# D<- solve(fit$Hess$Hess.full[fit$Hess$incla, fit$Hess$incla])
# covb <- (D%*%C)%*%(A)%*%(B%*%t(D))
# se <- (diag(abs(covb)))
# if(fit$row.eff == "random") {
# CovArerr <- matrix(se[1:length(fit$params$row.params)]);
# se<-se[-(1:length(fit$params$row.params))]
# }
# CovABerr<-array(0, dim(fit$Ab))
# for (i in 1:ncol(fit$Ab)) {
# CovABerr[,i,i] <- se[1:p]; se<-se[-(1:p)]
# }
# if((fit$num.lv+fit$num.lv.c) > 0) {
# CovAerr<-array(0, dim(fit$A))
# for (i in 1:ncol(fit$A)) {
# CovAerr[,i,i] <- se[1:n]; se<-se[-(1:n)]
# }
# }
# CovABerr
# }
CMSEPf <- function(fit, return.covb = F, type = NULL){
#for num.RR we are treating the LV as zero
n<-nrow(fit$y)
p<-ncol(fit$y)
num.lv <- fit$num.lv
num.lv.c <- fit$num.lv.c
num.RR <- fit$num.RR
num.lv.cor <- fit$num.lvcor
randomB <- fit$randomB
incla<-rep(FALSE, length(fit$Hess$incl))
if(fit$row.eff == "random") incla[names(fit$TMBfn$par)%in%c("r0")] <- TRUE
if(!is.null(fit$randomX)) incla[names(fit$TMBfn$par)%in%c("Br")] <- TRUE
if((num.lv+num.lv.c)>0) incla[names(fit$TMBfn$par)%in%c("u")] <- TRUE
if(randomB!=FALSE)incla[names(fit$TMBfn$par)%in%c("b_lv")] <- TRUE
fit$Hess$incla <- incla
A<- fit$Hess$cov.mat.mod #
B<- fit$Hess$Hess.full[fit$Hess$incl, fit$Hess$incla]
C<- fit$Hess$Hess.full[fit$Hess$incla, fit$Hess$incl]
D<- fit$Hess$Hess.full[fit$Hess$incla, fit$Hess$incla]
num.lv.cor <- fit$num.lvcor
radidx <- 0
if(is.null(type)){
if((num.lv.c+num.RR)==0){
type <- "residual"
}else{
type <- "conditional"
}
}
#### add errors for Bs if fixed ####
if(!fit$row.eff%in%c(FALSE,"fixed",NULL)){
radidx <- radidx +length(fit$TMBfn$par[names(fit$TMBfn$par)=="r0"])
}
if(!is.null(fit$randomX)){
radidx <-radidx+ length(fit$TMBfn$par[names(fit$TMBfn$par)=="Br"])
}
if(!is.null(fit$lv.X)&randomB!=FALSE){
radidx <-radidx+ length(fit$TMBfn$par[names(fit$TMBfn$par)=="b_lv"])
}
if((num.lv+num.lv.c+radidx)>0){
D <- solve(D)
}
sigma.lv <- fit$params$sigma.lv
#we never scale the prediction regions for num.lv by sigma.lv.
if(num.lv.c>0&num.lv>0){
sigma.lv <- c(sigma.lv[1:num.lv.c],rep(1,num.RR*ifelse(randomB==FALSE,1,0)),rep(1,num.lv))
}else if(num.lv>0&num.lv.c==0){
sigma.lv <- c(rep(1,num.RR*ifelse(randomB==FALSE,1,0)),rep(1,num.lv))
}else if(num.lv.c>0&num.lv==0){
sigma.lv <- c(sigma.lv,rep(1,num.RR*ifelse(randomB==FALSE,1,0)))
}else if(num.RR>0){
sigma.lv <- rep(1,num.RR*ifelse(randomB==FALSE,1,0))
}
if(prod(dim(D))!=0){colnames(D)<-row.names(D)<-names(fit$TMBfn$par[fit$Hess$incla])}
if(radidx>0){
if(prod(dim(D))==0&num.RR>0){
D <- matrix(0,ncol=num.RR*n,nrow=num.RR*n)
B <- matrix(0,nrow=sum(fit$Hess$incl),ncol=num.RR*n)
C <- t(B)
}else if(num.RR>0&randomB==FALSE){
D <- cbind(D[,1:(num.lv.c*n+radidx)],matrix(0,nrow=nrow(D),ncol=num.RR*n),D[,-c(1:(num.lv.c*n+radidx))])
D <- rbind(D[1:(num.lv.c*n+radidx),],matrix(0,nrow=num.RR*n,ncol=ncol(D)),D[-c(1:(num.lv.c*n+radidx)),])
B <- cbind(B[,1:(num.lv.c*n+radidx)],matrix(0,nrow=sum(fit$Hess$incl),ncol=num.RR*n),B[,-c(1:(num.lv.c*n+radidx))])
C <- t(B)
}
}else{
if(prod(dim(D))==0&num.RR>0){
D <- matrix(0,ncol=num.RR*n,nrow=num.RR*n)
B <- matrix(0,nrow=sum(fit$Hess$incl),ncol=num.RR*n)
C <- t(B)
}else if(num.lv.c>0&num.RR>0&randomB==FALSE){
D <- cbind(D[,1:(num.lv.c*n)],matrix(0,nrow=nrow(D),ncol=num.RR*n),D[,-c(1:(num.lv.c*n))])
D <- rbind(D[1:(num.lv.c*n),],matrix(0,nrow=num.RR*n,ncol=ncol(D)),D[-c(1:(num.lv.c*n)),])
B <- cbind(B[,1:(num.lv.c*n)],matrix(0,nrow=sum(fit$Hess$incl),ncol=num.RR*n),B[,-c(1:(num.lv.c*n))])
C <- t(B)
}else if(num.lv>0&num.RR>0&randomB==FALSE){
D <- cbind(matrix(0,nrow=nrow(D),ncol=num.RR*n),D)
D <- rbind(matrix(0,nrow=num.RR*n,ncol=ncol(D)),D)
B <- cbind(matrix(0,nrow=sum(fit$Hess$incl),ncol=num.RR*n),B)
C <- t(B)
}
}
if(num.RR>0&randomB==FALSE){
if((num.lv+num.lv.c+radidx)==0){colnames(D)<-rep("",ncol(D))}
colnames(D)[colnames(D)==""]<-"XB"
row.names(D)<-colnames(D)
}
if(num.lv.cor>0)
Q <- matrix(0,nrow=sum(names(fit$TMBfn$par)%in%c("u"))+radidx,ncol=dim(A)[1])
else
Q <- matrix(0,nrow=(num.lv+num.lv.c+num.RR*ifelse(randomB!=FALSE,0,1))*n+radidx,ncol=dim(A)[1])
# Q <- matrix(0,nrow=(num.lv+num.lv.c+num.RR)*n+radidx,ncol=dim(A)[1])
if((num.lv.c+num.RR)>0&randomB==FALSE){
for(q in 1:(num.lv.c+num.RR)){
Q[(1:n)+n*(q-1)+radidx,which(names(fit$TMBfn$par[fit$Hess$incl])=="b_lv")[(1:ncol(fit$lv.X))+(ncol(fit$lv.X)*(q-1))]] <- fit$lv.X#/fit$params$sigma.lv[q] #divide here to multiply later in ordiplot
}
}
####################################
if(type=="conditional"|type=="marginal"&randomB!=FALSE){
S <- diag(c(rep(1,radidx),rep(sigma.lv,each=n)))
covb <- (Q+S%*%D%*%C)%*%(A)%*%(t(Q)+B%*%t(D)%*%S)
}else if(type=="residual"){
covb <- (D%*%C)%*%(A)%*%(B%*%t(D))
}else if(type=="marginal"&randomB==FALSE){
covb <- Q%*%(A)%*%t(Q)
}
colnames(covb)<-row.names(covb)<-colnames(D)
#separate errors row-effects
ser0 <- diag(covb[colnames(covb)=="r0",colnames(covb)=="r0"])
covb <- covb[colnames(covb)!="r0",colnames(covb)!="r0"]
out <- list()
if(fit$row.eff == "random") {
CovArerr <- matrix(ser0);
out$Ar <- CovArerr
}
#separate errors b_lv
if(fit$randomB!=FALSE){
seb_lv <- diag(covb[colnames(covb)=="b_lv",colnames(covb)=="b_lv"])
covb_lvErr <- matrix(seb_lv,ncol=num.lv.c+num.RR)
out$Ab_lv <- covb_lvErr
covsB <- covb[colnames(covb)=="b_lv",colnames(covb)=="b_lv"]
}
if(!return.covb)covb <- covb[colnames(covb)!="b_lv",colnames(covb)!="b_lv"]
# separate errors AB
seAb <- diag(covb[colnames(covb)=="Br",colnames(covb)=="Br"])
covb <- covb[colnames(covb)!="Br",colnames(covb)!="Br"]
if(!is.null(fit$randomX)){
CovABerr<-array(0, dim(fit$Ab))
for (i in 1:ncol(fit$Ab)) {
CovABerr[,i,i] <- seAb[1:p]; seAb<-seAb[-(1:p)]
}
out$Ab <- CovABerr
}
if(!return.covb){
#re-order, select submatrices
try({
if(num.lv.cor>0){
# nS<- nrow(fit$TMBfn$env$parameters$u)
nS<- nrow(fit$A)
se <- simplify2array(sapply(1:nS,function(i)covb[seq(i,nS*(num.lv+num.lv.c+num.RR),by=nS),seq(i,nS*(num.lv+num.lv.c+num.RR),by=nS)],simplify=F))
} else if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(randomB==FALSE&type!="residual",1,0))>0){
se <- simplify2array(sapply(1:n,function(i)covb[seq(i,n*(num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(randomB==FALSE&type!="residual",1,0)),by=n),seq(i,n*(num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(randomB==FALSE&type!="residual",1,0)),by=n)],simplify=F))
}
if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(randomB==FALSE&type!="residual",1,0))>1){
se <- aperm(se,c(3,2,1))
}else if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(randomB==FALSE&type!="residual",1,0))>0 ){
se <- as.matrix(se)
}
# else if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(randomB==FALSE&type!="residual",1,0))>0 & (num.lv.cor >0 & (fit$corP$cstruc[2] !="diag"))){
# se <- diag(se)
# }
#add error for Bs if random
if((num.RR+num.lv.c)>0&type!="residual"&randomB!=FALSE){
if((num.lv+num.lv.c+num.RR)>0){
se2 <- array(0,dim=c(n,num.lv.c+num.RR+num.lv*ifelse(type=="marginal",0,1),num.lv.c+num.RR+num.lv*ifelse(type=="marginal",0,1)))
if(type=="conditional"&(num.lv+num.lv.c)>0&num.RR>0){
se2[,-c((num.lv.c+1):(num.lv.c+num.RR)),-c((num.lv.c+1):(num.lv.c+num.RR))] <- se #0s for RRR
}else if(type=="conditional"&(num.lv+num.lv.c)>0){
se2 <- se
}
se <- se2
}
for(i in 1:n){
Q <- as.matrix(Matrix::bdiag(replicate(num.RR+num.lv.c,fit$lv.X[i,,drop=F],simplify=F)))
temp <- Q%*%covsB%*%t(Q)
temp[col(temp)!=row(temp)] <- 2*temp[col(temp)!=row(temp)] ##should be double the covariance
se[i,1:(num.RR+num.lv.c),1:(num.RR+num.lv.c)] <- se[i,1:(num.RR+num.lv.c),1:(num.RR+num.lv.c)] + temp
}
}
if((num.lv*ifelse(type=="marginal",0,1)+num.lv.c+num.RR*ifelse(type!="residual",1,0)) > 0) {
out$A <- se
}
},silent=T)
}else{
out <- covb
}
return(out)
}
# draw an ellipse
ellipse<-function(center, covM, rad, col=4, lwd = 1, lty = 1){
seg <- 51
Qc <- chol(covM, pivot = TRUE)
angles <- (0:seg) * 2 * pi / seg
unit.circ <- cbind(cos(angles), sin(angles))
order <- order(attr(Qc, "pivot"))
ellips <- t(center + rad * t(unit.circ %*% Qc[, order]))
lines(ellips, col = col, lwd = lwd, lty = lty)
}
gamEnvelope <- function(x, y,line.col = "red", envelope.col = c("blue","lightblue"), col = 1, envelopes = TRUE, ...){
xSort <- sort(x, index.return = TRUE)
gam.yx <- gam(y[xSort$ix] ~ xSort$x)
pr.y <- predict.gam(gam.yx, se.fit = TRUE)
n.obs <- length(xSort$ix)
prHi <- pr.y$fit + 1.96*pr.y$se.fit
prLow <- pr.y$fit - 1.96*pr.y$se.fit
if(envelopes) polygon(xSort$x[c(1:n.obs,n.obs:1)], c(prHi,prLow[n.obs:1]), col = envelope.col[2], border = NA)
lines(xSort$x, pr.y$fit, col = envelope.col[1])
abline(h = 0, col = 1)
points(x, y, col = col, ...)
}
mlm <- function(y, X = NULL, index = NULL){
out <- list(coefficients=NULL, residuals = NULL)
coefficients <- residuals <- NULL
if(!is.null(X) || !is.null(index)){
Xx<- cbind(X, index)
for (j in 1:ncol(y)) {
f0 <- lm(y[,j]~Xx)
coefficients <- rbind(coefficients, f0$coef)
residuals <- cbind(residuals, f0$residuals)
}
} else {
for (j in 1:ncol(y)) {
f0 <- lm(y[,j]~1)
coefficients <- rbind(coefficients, f0$coef)
residuals <- cbind(residuals, f0$residuals)
}
}
out$coefficients <- t(coefficients)
out$residuals <- residuals
out
}
RRse <- function(object){
if(object$quadratic!=FALSE){
warning("Coefplot for quadratic coefficients not yet implemented.")
}
if(!is.list(object$sd)){
stop("Cannot construct coefplot without standard errors in the model.")
}
K <- ncol(object$lv.X)
d <- object$num.RR+object$num.lv.c
p <- ncol(object$y)
#still add RR or num.lv.c or something to predictor names
beta <- object$params$theta[,1:d,drop=F]%*%t(object$params$LvXcoef)
betaSE<-matrix(0,ncol=K,nrow=ncol(object$y))
colnames(betaSE)<-colnames(object$lv.X)
row.names(betaSE)<-colnames(object$y)
if(isFALSE(object$randomB)){
covMat <- object$Hess$cov.mat.mod
colnames(covMat) <- row.names(covMat) <- names(object$TMBfn$par[object$Hess$incl])
covMat <- covMat[colnames(covMat)%in%c("b_lv","lambda"),colnames(covMat)%in%c("b_lv","lambda")]
#add first row and column of zeros before b_lv, for first species
covMat <- rbind(covMat[1:(d*K),, drop=FALSE],0,covMat[-c(1:(d*K)),, drop=FALSE])
covMat <- cbind(covMat[,1:(d*K), drop=FALSE],0,covMat[,-c(1:(d*K)), drop=FALSE])
if(d>1){
idx<-which(c(upper.tri(object$params$theta[,1:d],diag=T)))[-1]
#add zeros where necessary
for(q in 1:length(idx)){
covMat <- rbind(covMat[1:(d*K+idx[q]-1),],0,covMat[(d*K+idx[q]):ncol(covMat),])
covMat <- cbind(covMat[,1:(d*K+idx[q]-1)],0,covMat[,(d*K+idx[q]):ncol(covMat)])
}
}
row.names(covMat)[row.names(covMat)==""]<-colnames(covMat)[colnames(covMat)==""]<-"lambda"
covB <- covMat[colnames(covMat)=="b_lv",colnames(covMat)=="b_lv", drop=FALSE]
covL <- covMat[colnames(covMat)=="lambda",colnames(covMat)=="lambda", drop=FALSE]
covLB <- covMat[colnames(covMat)=="lambda",colnames(covMat)=="b_lv", drop=FALSE]
}else{
if(object$method=="LA"){
covMat <- object$Hess$cov.mat.mod
colnames(covMat) <- row.names(covMat) <- names(object$TMBfn$par)[object$Hess$incl]
covMat <- covMat[colnames(covMat)%in%"lambda",colnames(covMat)%in%"lambda"]
#add first row and column of zeros before b_lv, for first species
covMat <- rbind(covMat,0,covMat)
covMat <- cbind(covMat,0,covMat)
if(d>1){
idx<-which(c(upper.tri(object$params$theta[,1:d],diag=T)))[-1]
#add zeros where necessary
for(q in 1:length(idx)){
covMat <- rbind(covMat[1:(idx[q]-1),],0,covMat[(idx[q]):ncol(covMat),])
covMat <- cbind(covMat[,1:(idx[q]-1)],0,covMat[,(idx[q]):ncol(covMat)])
}
}
row.names(covMat)[row.names(covMat)==""]<-colnames(covMat)[colnames(covMat)==""]<-"lambda"
covB <- sdrandom(object$TMBfn, object$Hess$cov.mat.mod, object$Hess$incl,ignore.u = F, type = "residual", return.covb=T)
covB <- covB[colnames(covB)=="b_lv",colnames(covB)=="b_lv"]
covL <- covMat[colnames(covMat)=="lambda",colnames(covMat)=="lambda", drop=FALSE]
# getting cov(b,gamma) is more difficult for Laplace because we do not have the full hessian
# the solution below goes via the jointPrecision matrix from TMB
r <- object$TMBfn$env$random
par = object$TMBfn$env$last.par
hessian.random <- object$TMBfn$env$spHess(par, random = TRUE)
f <- object$TMBfn$env$f
w <- rep(0, length(par))
reverse.sweep <- function(i) {
w[i] <- 1
f(par, order = 1, type = "ADGrad", rangeweight = w, doforward = 0)[r]
}
nonr <- setdiff(seq_along(par), r)
tmp <- sapply(nonr, reverse.sweep)
if(!is.matrix(tmp)) tmp <- matrix(tmp, ncol=length(nonr) )
A <- solve(hessian.random, tmp[, object$Hess$incl])
G <- hessian.random %*% A
G <- as.matrix(G)
M1 <- cbind2(hessian.random,G)
Vtheta <- object$Hess$Hess.full[object$Hess$incl,object$Hess$incl]
M2 <- cbind2(t(G), as.matrix(t(A)%*%G)+Vtheta)
M <- rbind2(M1,M2)
dn <- c(names(par)[r],names(par[-r])[object$Hess$incl])
dimnames(M) <- list(dn,dn)
ip <- Matrix::invPerm(c(r,(1:length(par))[-r][object$Hess$incl]))
jointPrecision <- M[ip,ip]
# bottom-left block of covariance matrix via block inversion
incl = row.names(jointPrecision)%in%names(object$TMBfn$env$last.par.best[-r][object$Hess$incl])
incld = row.names(jointPrecision)[r]
covMat <- try(-solve(as.matrix(jointPrecision[incld,incld]),as.matrix(jointPrecision[incld,incl]))%*%object$Hess$cov.mat.mod,silent=T)
if(inherits(covMat,"try-error")){
# Via fixed-effects part of Hessian if random-effects part is singular
Ai <- solve(as.matrix(jointPrecision[incl,incl]))
B.mat <- as.matrix(jointPrecision[incld,incl])
D.mat <- as.matrix(jointPrecision[incld,incld])
covMat<- -MASS::ginv(-D.mat-B.mat%*%Ai%*%t(B.mat))%*%B.mat%*%Ai
}
row.names(covMat) <- names(object$TMBfn$env$last.par.best)[r]
colnames(covMat) <- names(object$TMBfn$env$par)[-r][object$Hess$incl]
covMat <- covMat[row.names(covMat)=="b_lv",colnames(covMat) == "lambda"]
#add first column of zeros for first species
covMat <- cbind(0,covMat)
if(d>1){
idx<-which(c(upper.tri(object$params$theta[,1:d],diag=T)))[-1]
#add zeros where necessary
for(q in 1:length(idx)){
covMat <- cbind(covMat[,1:(idx[q]-1)],0,covMat[,idx[q]:ncol(covMat)])
}
}
row.names(covMat)[row.names(covMat)==""]<-colnames(covMat)[colnames(covMat)==""]<-"lambda"
covLB <- t(covMat)
}else{
if(object$randomB=="P"|object$randomB=="single"){
covB <- as.matrix(Matrix::bdiag(lapply(seq(dim(object$Ab.lv)[1]), function(k) object$Ab.lv[k , ,])))
}else if(object$randomB=="LV"){
covB <- as.matrix(Matrix::bdiag(lapply(seq(dim(object$Ab.lv)[1]), function(q) object$Ab.lv[q , ,])))
}
covsB <- CMSEPf(object, return.covb = T)
covB = covB + covsB[row.names(covsB)=="b_lv",colnames(covsB)=="b_lv"]
# covariance matrix of loadings
covL <- object$Hess$cov.mat.mod;colnames(covL) <- names(object$TMBfn$par)[object$Hess$incl];covL <- covL[colnames(covL)=="lambda",colnames(covL)=="lambda", drop=FALSE]
# covariance of parameters via block inversion
incl=object$Hess$incl;incld=object$Hess$incld
covMat <- -solve(object$Hess$Hess.full[incld,incld], object$Hess$Hess.full[incld,incl])%*%object$Hess$cov.mat.mod
if(inherits(covMat,"try-error")){
# Via fixed-effects part of Hessian if random-effects part is singular
Ai <- solve(object$Hess$Hess.full[incl,incl])
B.mat <- object$Hess$Hess.full[incld,incl]
D.mat <- object$Hess$Hess.full[incld,incld]
covMat<- -MASS::ginv(-D.mat-B.mat%*%Ai%*%t(B.mat))%*%B.mat%*%Ai
}
row.names(covMat) <- names(object$TMBfn$par)[incld]
colnames(covMat) <- names(object$TMBfn$par)[incl]
covMat <- covMat[row.names(covMat)=="b_lv",colnames(covMat) == "lambda"]
#add first column of zeros for first species
covMat <- cbind(0,covMat)
covL <- rbind(0,cbind(0,covL))
if(d>1){
idx<-which(c(upper.tri(object$params$theta[,1:d],diag=T)))[-1]
#add zeros where necessary
for(q in 1:length(idx)){
covMat <- cbind(covMat[,1:(idx[q]-1)],0,covMat[,idx[q]:ncol(covMat)])
covL <- rbind(covL[1:(idx[q]-1),],0,covL[(idx[q]):ncol(covL),])
covL <- cbind(covL[,1:(idx[q]-1)],0,covL[,(idx[q]):ncol(covL)])
}
}
row.names(covMat)[row.names(covMat)==""]<-colnames(covMat)[colnames(covMat)==""]<-"lambda"
covLB <- t(covMat)
}
}
#all ordered by LV, so LV11 LV12 LV13 etc
for(k in 1:K){
for(j in 1:p){
for(q in 1:d){
for(r in 1:d){
betaSE[j,k] <- betaSE[j,k]+object$params$LvXcoef[k,q]*object$params$LvXcoef[k,r]*covL[(q-1)*p+j,(r-1)*p+j]+
object$params$LvXcoef[k,q]*object$params$theta[j,r]*covLB[(q-1)*p+j,(r-1)*K+k]+
object$params$LvXcoef[k,r]*object$params$theta[j,q]*covLB[(r-1)*p+j,(q-1)*K+k]+
object$params$theta[j,q]*object$params$theta[j,r]*covB[(r-1)*K+k,(q-1)*K+k]+
covB[(r-1)*K+k,(q-1)*K+k]*covL[(q-1)*p+j,(r-1)*p+j]+
covLB[(r-1)*p+j,(q-1)*K+k]*covLB[(q-1)*p+j,(r-1)*K+k]
}
}
}
}
betaSE<-sqrt(abs(betaSE))
return(betaSE)
}
#functions for optimization with equality constraints using alabama
#constraint
eval_eq_c <- function(x, obj, ...){ #alabama requires a ... argument
B <- matrix(x[names(obj$par)=="b_lv"],ncol=obj$env$data$num_lv_c+obj$env$data$num_RR)
con<-NULL
d <- obj$env$data$num_lv_c+obj$env$data$num_RR
nc <- d*(d-1)/2# number of constraints
combs <- combn(1:(obj$env$data$num_RR+obj$env$data$num_lv_c),2)
con <- colSums(B[,combs[1,],drop=F]*B[,combs[2,],drop=F])
return(con)
}
#jacobian
eval_eq_j <- function(x, obj, ...){
B <- matrix(x[names(obj$par)=="b_lv"],ncol=obj$env$data$num_lv_c+obj$env$data$num_RR)
jacob <- NULL
d <- obj$env$data$num_lv_c+obj$env$data$num_RR
nc <- d*(d-1)/2#number of constraints
#LV constraint combinations
combs <- combn(1:(obj$env$data$num_RR+obj$env$data$num_lv_c),2)
#Build jacobian
jacob.B <- matrix(0,nrow=nc,ncol=sum(names(obj$par)=="b_lv"))
#Indices
idx.i <- rep(1:nc,each=nrow(B))
idx.j <- nrow(B)*(rep(combs[1,],each=nrow(B))-1)+rep(1:nrow(B),nc)
#d(constraint)/d(b1)
jacob.B[cbind(idx.i,idx.j)] <- c(B[,combs[2,]])
idx.j <- nrow(B)*(rep(combs[2,],each=nrow(B))-1)+rep(1:nrow(B),nc)
#d(constraint)/d(b2)
jacob.B[cbind(idx.i,idx.j)] <- c(B[,combs[1,]])
#padd with zeros for all unrelated parameters
jacob <- cbind(matrix(0,nrow=nc,ncol=which(names(obj$par)=="b_lv")[1]-1),jacob.B,matrix(0,nrow=nc,ncol=length(x)-tail(which(names(obj$par)=="b_lv"),1)))
return(jacob)
}
#functions for optimization with equality constraints using nloptr
eval_f<-function(x, obj = NULL){
#nloptr requires to return likelihood and gradient simultaneously
return(list("objective" = obj$fn(x),
"gradient" = obj$gr(x)))
}
eval_g_eq <- function(x, obj = NULL){
B <- matrix(x[names(obj$par)=="b_lv"],ncol=obj$env$data$num_lv_c+obj$env$data$num_RR)
jacob <- NULL
con<-NULL
d <- obj$env$data$num_lv_c+obj$env$data$num_RR
nc <- d*(d-1)/2#number of constraints
combs <- combn(1:(obj$env$data$num_RR+obj$env$data$num_lv_c),2)
con <- colSums(B[,combs[1,],drop=F]*B[,combs[2,],drop=F])
jacob.B <- matrix(0,nrow=nc,ncol=sum(names(obj$par)=="b_lv"))
idx.i <- rep(1:nc,each=nrow(B))
idx.j <- nrow(B)*(rep(combs[1,],each=nrow(B))-1)+rep(1:nrow(B),nc)
jacob.B[cbind(idx.i,idx.j)] <- c(B[,combs[2,]])
idx.j <- nrow(B)*(rep(combs[2,],each=nrow(B))-1)+rep(1:nrow(B),nc)
jacob.B[cbind(idx.i,idx.j)] <- c(B[,combs[1,]])
jacob <- cbind(matrix(0,nrow=nc,ncol=which(names(obj$par)=="b_lv")[1]-1),jacob.B,matrix(0,nrow=nc,ncol=length(x)-tail(which(names(obj$par)=="b_lv"),1)))
#nloptr requires to return constraint and jacobian simultaneously
res <- list("constraints" = con,"jacobian" = jacob)
return(res)
}
# function to post-hoc estimate lagranian multipliers
# see https://discourse.julialang.org/t/lagrangian-function/38287/18?u=stevengj
lambda<-function(x,obj)c(-obj$gr()%*%MASS::ginv(eval_eq_j(x,obj = obj)))
# gradient of Lagranian
gradL <- function(x,lambda,objr){
res <- eval_g_eq(x,objr)
con <- res$con
con.j <- res$jacobian
objr$gr(x)-colSums(lambda(x,objr)*con.j)#+sig*drop(t(con.j)%*%con)
}
# 2nd derivative, for constraint for correction of SEs
eval_eq_j <- function(x, obj, ...){
B <- matrix(x[names(obj$par)=="b_lv"],ncol=obj$env$data$num_lv_c+obj$env$data$num_RR)
jacob <- NULL
d <- obj$env$data$num_lv_c+obj$env$data$num_RR
nc <- d*(d-1)/2#number of constraints
#LV constraint combinations
combs <- combn(1:(obj$env$data$num_RR+obj$env$data$num_lv_c),2)
#Build jacobian
jacob.B <- matrix(0,nrow=nc,ncol=sum(names(obj$par)=="b_lv"))
#Indices
idx.i <- rep(1:nc,each=nrow(B))
idx.j <- nrow(B)*(rep(combs[1,],each=nrow(B))-1)+rep(1:nrow(B),nc)
#d(constraint)/d(b1)
jacob.B[cbind(idx.i,idx.j)] <- c(B[,combs[2,]])
idx.j <- nrow(B)*(rep(combs[2,],each=nrow(B))-1)+rep(1:nrow(B),nc)
#d(constraint)/d(b2)
jacob.B[cbind(idx.i,idx.j)] <- c(B[,combs[1,]])
#padd with zeros for all unrelated parameters
jacob <- cbind(matrix(0,nrow=nc,ncol=which(names(obj$par)=="b_lv")[1]-1),jacob.B,matrix(0,nrow=nc,ncol=length(x)-tail(which(names(obj$par)=="b_lv"),1)))
return(jacob)
}
# this is dc(x)/dbdb, i.e., hessian w.r.t. the constraint function Lmult*(B^tB - I)
b_lvHEcorrect <- function(Lmult,K,d){
corHE <- matrix(0,K*d,K*d)
combs <- combn(1:d,2)
for(q in 1:ncol(combs)){
diag(corHE[(combs[1,q]*K-K+1):(combs[1,q]*K),(combs[2,q]*K-K+1):(combs[2,q]*K)])<- Lmult[q]
diag(corHE[(combs[2,q]*K-K+1):(combs[2,q]*K),(combs[1,q]*K-K+1):(combs[1,q]*K)])<- Lmult[q]
}
corHE
}
# distribution functions for ZIP andZINB
pzip <- function(y, mu, sigma)
{
pp <- NULL
if (y > -1) {
cdf <- ppois(y, lambda = mu, lower.tail = TRUE, log.p = FALSE)
cdf <- sigma + (1 - sigma) * cdf
pp <- cdf
}
if (y < 0) {
pp <- 0
}
pp
}
pzinb <- function(y, mu, sigma)
{
pp <- NULL
if (y > -1) {
cdf <- pnbinom(y, mu = mu, size = 1 / sigma)
cdf <- sigma + (1 - sigma) * cdf
pp <- cdf
}
if (y < 0) {
pp <- 0
}
pp
}
# function to get the derivative w.r.t. the squared constraint function
# eval_eq_j2 <- function(x, obj, ...){
# B <- matrix(x[names(obj$par)=="b_lv"],ncol=obj$env$data$num_lv_c+obj$env$data$num_RR)
# jacob <- NULL
#
# d <- obj$env$data$num_lv_c+obj$env$data$num_RR
# nc <- d*(d-1)/2#number of constraints
#
# #LV constraint combinations
# combs <- combn(1:(obj$env$data$num_RR+obj$env$data$num_lv_c),2)
#
# #Build jacobian
# jacob.B <- matrix(0,nrow=nc,ncol=sum(names(obj$par)=="b_lv"))
# #Indices
# idx.i <- rep(1:nc,each=nrow(B))
# idx.j <- nrow(B)*(rep(combs[1,],each=nrow(B))-1)+rep(1:nrow(B),nc)
# #d(constraint^2)/d(b1)
# jacob.B[cbind(idx.i,idx.j)] <- 2*B[,combs[2,]]^2*B[,combs[1,]]
#
# idx.j <- nrow(B)*(rep(combs[2,],each=nrow(B))-1)+rep(1:nrow(B),nc)
# #d(constraint^2)/d(b2)
# jacob.B[cbind(idx.i,idx.j)] <- 2*B[,combs[1,]]^2*B[,combs[2,]]
# #padd with zeros for all unrelated parameters
# jacob <- cbind(matrix(0,nrow=nc,ncol=which(names(obj$par)=="b_lv")[1]-1),jacob.B,matrix(0,nrow=nc,ncol=length(x)-tail(which(names(obj$par)=="b_lv"),1)))
# return(jacob)
# }
# function adapted from utils package to re-order gllvm's parameters and/or standard errors
relist.gllvm <- function (flesh, skeleton = attr(flesh, "skeleton"))
{
ind <- 1L
nam = unique(names(flesh))
skeleton <- result <- skeleton[nam]
for (i in nam) {
skel_i <- result[[i]]
size <- length(unlist(skel_i))
result[[i]] <- relist(flesh[seq.int(ind, length.out = size)],
skel_i)
ind <- ind + size
}
result
}
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.