Merge pull request #151 from BertvanderVeen/coleff

clean-up

NAMESPACE

@@ -86,20 +86,24 @@ importFrom(graphics,points)
 importFrom(graphics,polygon)
 importFrom(graphics,segments)
 importFrom(graphics,text)
+importFrom(methods,as)
 importFrom(methods,cbind2)
+importFrom(methods,is)
 importFrom(methods,rbind2)
 importFrom(mgcv,gam)
 importFrom(mgcv,predict.gam)
 importFrom(mvabund,manyglm)
 importFrom(nloptr,nloptr)
 importFrom(stats,AIC)
 importFrom(stats,BIC)
+importFrom(stats,aggregate)
 importFrom(stats,as.formula)
 importFrom(stats,binomial)
 importFrom(stats,coef)
 importFrom(stats,confint)
 importFrom(stats,constrOptim)
 importFrom(stats,cov)
+importFrom(stats,cov2cor)
 importFrom(stats,dbinom)
 importFrom(stats,dnorm)
 importFrom(stats,factanal)
@@ -135,6 +139,7 @@ importFrom(stats,qqnorm)
 importFrom(stats,quantile)
 importFrom(stats,rbeta)
 importFrom(stats,rbinom)
+importFrom(stats,reformulate)
 importFrom(stats,reshape)
 importFrom(stats,residuals)
 importFrom(stats,residuals.lm)
@@ -146,8 +151,10 @@ importFrom(stats,rpois)
 importFrom(stats,runif)
 importFrom(stats,sd)
 importFrom(stats,terms)
+importFrom(stats,update.formula)
 importFrom(stats,vcov)
 importFrom(utils,combn)
+importFrom(utils,head)
 importFrom(utils,relist)
 importFrom(utils,tail)
 useDynLib(gllvm, .registration = TRUE)

R/TMBtrait.R

@@ -1302,7 +1302,7 @@ trait.TMB <- function(
 
         parameter.list = list(r0=r1, b = b1, b_lv = matrix(0), sigmab_lv = 0, Ab_lv = 0, B=B1, Br=Br1, lambda = lambda1, lambda2 = t(lambda2), sigmaLV = sigma.lv1, u = u1, lg_phi=lg_phi1, sigmaB=sigmaB1, sigmaij=sigmaij1, log_sigma=log_sigma1, rho_lvc=rho_lvc, Au=Au, lg_Ar=lg_Ar, Abb=Abb, zeta=zeta, ePower = ePower, lg_phiZINB = lg_phiZINB1) #, scaledc=scaledc, thetaH = thetaH, bH=bH
 
-      data.list = list(y = y, x = Xd, x_lv = matrix(0), xr=xr, xb=xb, dr0 = dr, dLV = dLV, offset=offset, nr = nr, num_lv = num.lv, num_RR = 0, num_lv_c = 0, quadratic = 1, family=familyn, extra=extra, method=switch(method, VA=0, EVA=2), model=1, random=randoml, zetastruc = ifelse(zeta.struc=="species",1,0), times = times, cstruc=cstrucn, cstruclv = cstruclvn, dc=dist, dc_lv = distLv, Astruc=Astruc, NN = NN, Ntrials = Ntrials, colMatBlocksI = blocks,  Abranks = Abranks, Abstruc = Abstruc, cs = matrix(0), nncolMat = nncolMat)
+      data.list = list(y = y, x = Xd, x_lv = matrix(0), xr=xr, xb=xb, dr0 = dr, dLV = dLV, offset=offset, nr = nr, num_lv = num.lv, num_RR = 0, num_lv_c = 0, quadratic = 1, family=familyn, extra=extra, method=switch(method, VA=0, EVA=2), model=1, random=randoml, zetastruc = ifelse(zeta.struc=="species",1,0), times = times, cstruc=cstrucn, cstruclv = cstruclvn, dc=dist, dc_lv = distLV, Astruc=Astruc, NN = NN, Ntrials = Ntrials, colMatBlocksI = blocks,  Abranks = Abranks, Abstruc = Abstruc, cs = matrix(0), nncolMat = nncolMat)
 
       objr <- TMB::MakeADFun(
         data = data.list, silent=!trace,

R/getEnvironCov.gllvm.R

@@ -42,12 +42,12 @@
 #'# As proportion of variance in the model
 #'envcov$trace.randomB/sum(envcov$trace.randomB)
 #'}
-#'@aliases getEnvironCov getEnvironCor
+#'@aliases getEnvironCov getEnvironCor getEnvironCor.gllvm
 #'@method getEnvironCov gllvm
 #'@export
 #'@export getEnvironCov.gllvm
 
-getEnvironCov.gllvm <- function(object, x = NULL){
+getEnvironCov.gllvm <- function(object, x = NULL, ...){
 
   if(isFALSE(object$randomB) & isFALSE(object$col.eff$col.eff)& is.null(object$randomX)){
     stop("Canot calculate correlations without random effects for covariates in the model.")
@@ -96,7 +96,7 @@ if(object$col.eff$col.eff=="random"){
 
   if(!isFALSE(object$quadratic)){
     # add tr(D_jSigma_zD_kSigma_z) for z = B^t x
-    Sigmaz <- diag(rep(x%*%diag(model2$params$sigmaLvXcoef^2)%*%x,object$num.lv.c+object$num.RR))
+    Sigmaz <- diag(rep(x%*%diag(object$params$sigmaLvXcoef^2)%*%x,object$num.lv.c+object$num.RR))
     theta <- object$params$theta[,-c(1:(object$num.lv.c+object$num.RR+object$num.lv))][,1:(object$num.lv.c+object$num.RR)]
     cov.environ.randomB.quad <- 2*abs(theta)%*%Sigmaz%*%Sigmaz%*%t(abs(theta))
     trace.environ.randomB.quad <- lapply(cov.environ.randomB.quad, function(x)sum(diag(x)))

R/getPredictErr.gllvm.R

@@ -127,7 +127,7 @@ getPredictErr.gllvm = function(object, CMSEP = TRUE, cov = FALSE, ...)
 
       if((num.lv+num.lv.c)>0){ object$A<-sdb$A+A} else{object$A <- sdb$A}
       if(num.RR>0&object$randomB!=FALSE){
-        if(object$randomB=="P"|object$randomB=="single"|randomB=="iid"){
+        if(object$randomB=="P"|object$randomB=="single"|object$randomB=="iid"){
           covsB <- as.matrix(Matrix::bdiag(lapply(seq(dim(object$Ab.lv)[1]), function(k) object$Ab.lv[k , ,])))
         }else if(object$randomB=="LV"){
           covsB <- as.matrix(Matrix::bdiag(lapply(seq(dim(object$Ab.lv)[1]), function(q) object$Ab.lv[q , ,])))
@@ -146,7 +146,7 @@ getPredictErr.gllvm = function(object, CMSEP = TRUE, cov = FALSE, ...)
       sdb <- list(Ab_lv = 0)
 
       if(num.RR>0&object$randomB!=FALSE){
-        if(object$randomB=="P"|object$randomB=="single"|randomB=="iid"){
+        if(object$randomB=="P"|object$randomB=="single"|object$randomB=="iid"){
           covsB <- as.matrix(Matrix::bdiag(lapply(seq(dim(object$Ab.lv)[1]), function(k) object$Ab.lv[k , ,])))
         }else if(object$randomB=="LV"){
           covsB <- as.matrix(Matrix::bdiag(lapply(seq(dim(object$Ab.lv)[1]), function(q) object$Ab.lv[q , ,])))

R/gllvm.R

@@ -405,14 +405,14 @@
 #'@importFrom mvabund manyglm
 #'@importFrom graphics abline axis par plot segments text points boxplot panel.smooth lines polygon arrows 
 #'@importFrom grDevices rainbow hcl
-#'@importFrom stats dnorm pnorm qnorm rnorm dbinom pbinom rbinom pnbinom rnbinom pbeta rbeta pexp rexp pgamma rgamma ppois rpois runif pchisq qchisq qqnorm lm AIC binomial constrOptim factanal glm model.extract model.frame model.matrix model.response nlminb optim optimHess reshape residuals terms BIC qqline sd formula ppoints quantile gaussian cov princomp as.formula residuals.lm coef printCoefmat nobs predict
+#'@importFrom stats dnorm pnorm qnorm rnorm dbinom pbinom rbinom pnbinom rnbinom pbeta rbeta pexp rexp pgamma rgamma ppois rpois runif pchisq qchisq qqnorm lm AIC binomial constrOptim factanal glm model.extract model.frame model.matrix model.response nlminb optim optimHess reshape residuals terms BIC qqline sd formula ppoints quantile gaussian cov princomp as.formula residuals.lm coef printCoefmat nobs predict cov2cor reformulate update.formula aggregate
 #'@importFrom Matrix bdiag chol2inv diag
 #'@importFrom MASS ginv polr mvrnorm
 #'@importFrom mgcv gam predict.gam
 #'@importFrom nloptr nloptr
 #'@importFrom alabama auglag
-#'@importFrom utils combn tail relist
-#'@importFrom methods cbind2 rbind2
+#'@importFrom utils combn tail relist head
+#'@importFrom methods cbind2 rbind2 is as
 #'
 
 gllvm <- function(y = NULL, X = NULL, TR = NULL, data = NULL, formula = NULL, family,
@@ -718,7 +718,7 @@ gllvm <- function(y = NULL, X = NULL, TR = NULL, data = NULL, formula = NULL, fa
       }else{
         X <- cbind(X, X.col.eff)
         colnames(X)[tail(1:ncol(X),ncol(X.col.eff))] <-  make.names(paste0("RE_mean_", make.unique(colnames(X.col.eff))))
-        formula <- update(as.formula(formula), as.formula(paste0("~. + ", paste0("RE_mean_", make.unique(colnames(X.col.eff)), collapse = "+"))))
+        formula <- update.formula(as.formula(formula), as.formula(paste0("~. + ", paste0("RE_mean_", make.unique(colnames(X.col.eff)), collapse = "+"))))
       }
     }
 
@@ -1104,7 +1104,7 @@ gllvm <- function(y = NULL, X = NULL, TR = NULL, data = NULL, formula = NULL, fa
       num.lv.cor = 0
     }
     if(Lambda.struc %in% c("bdNN","UNN") & num.lv.cor>0){
-      NN<-t(apply(as.matrix(dist(distLv, upper = TRUE, diag = TRUE)),1, order)[1+(1:NN),])
+      NN<-t(apply(as.matrix(dist(distLV, upper = TRUE, diag = TRUE)),1, order)[1+(1:NN),])
       i1<-rep(1:nrow(NN), each=ncol(NN))
       i2<-c(t(NN))
       indM<-cbind(i1,i2)

R/predict.gllvm.R

@@ -314,7 +314,7 @@ predict.gllvm <- function(object, newX = NULL, newTR = NULL, newLV = NULL, type
     }
     r0 <- object$params$row.params
     if((object$row.eff %in% "random") && (level==0)) r0 = r0*0
-    eta <- eta + r0%*%rep(1,ncol(model$y))
+    eta <- eta + r0%*%rep(1,ncol(object$y))
   }
 
   if (object$col.eff$col.eff == "random" && is.null(newX)) {

R/residuals.gllvm.R

@@ -86,7 +86,7 @@ residuals.gllvm <- function(object, ...) {
       if (object$family == "gamma") {
         phis <- object$params$phi # - 1
         b <- pgamma(as.vector(y), shape = rep(phis, each = n), scale = as.vector(mu)/rep(phis, each = n))
-        a <- pmin(b, pgamma(as.vector(y)), shape = rep(phis, each = n), scale = mas.vector(mu)/rep(phis, each = n))
+        a <- pmin(b, pgamma(as.vector(y)), shape = rep(phis, each = n), scale = as.vector(mu)/rep(phis, each = n))
 
         u <- runif(n*p, min = a, max = b)
         if(any(u==1))u[u==1] <- 1-1e-16