post-hoc test

[AlainZuur]

I'm preparing an exercise in which I model 16 morphometric variables using GLLVM with a Gamma distribution...honeybee data. With one categorical covariate "Location" with 8 levels, a row-specific random effect, and 2 latent variables. All very nice....works perfectly.

I can see the following question coming: Can we test which locations differ from each other? So..that is post-hoc testing territory.

I guess you have to do that manually?

Alain

[BertvanderVeen]

See my answer of this morning in #157 . In short: yes, you need to do it manually, but I provide some code there that you can adjust to help you do it.

[AlainZuur]

Thanks....!

Alain

[jebyrnes]

Oh, can this be made a vignette? Posthocs and comparisons for individual species are still something that I struggle with - and struggle with teaching! Would be very useful!

[BertvanderVeen]

Happy to, but I'm too busy for the moment. Perhaps @AlainZuur is willing to contribute some of the code for his exercise once finished? That would make things easier.

[AlainZuur]

yes...sure. Not a problem. It is a Gamma GLLVM with as only covariate factor Location with (8!) levels plus some LVs and random effects.

The gllvm output gives me the intercepts (which is Location 1 in this specific model) and the coefplot gives me the correction of the intercept for each other level. I wanted to have "Intercept + Corrections of the intercept", plus its SE. So that I get the actual Location effect.

I managed to calculate C * beta where C is for example C = (1, 1, 0, 0, 0, 0, 0) and the corresponding standard errors....for each level of Location and each response variable. Horrible to program...and it would be horrible to write generic code for that.

Doing a post-hoc test would mean using C = (1, -1, 0, 0, 0, 0, 0) and repeating all the code...and then C = (1, 0, -1, 0, 0, 0, 0)...and then C = (1, 0, 0, -1, 0, 0, 0)...etc etc.

According to ChatGPT that gives me 8 * (8-1) / 2= 28 unique comparisons per variable. And there are 16 response variables....adding up to 28 * 16 = 448 comparisons. I hardly ever do post-hoc tests, but it tells me to use the Benjamini-Hochberg procedure. It was friendly enough to give me R code for that.

But sure...once it is all worked out I can share it. Probably later on this week.

Alain

[jebyrnes]

That would be amazing! Note, I've signed up for the seminar and am urging folk I know to do so. VERY excited.

[BertvanderVeen]

That would be amazing! Note, I've signed up for the seminar and am urging folk I know to do so. VERY excited.

Which seminar is that?

[jebyrnes]

https://www.highstat.com/index.php/livevideo2

[AlainZuur]

That is in 2025!

[BertvanderVeen]

@jebyrnes this got me thinking. So I have something for you. It is only a start, and I don't have much time to improve it at present, so please don't get too excited! But I know you have been waiting for it for some time.

recover_data.gllvm <- function(object, component = c("main", "LV"), ...){
 component <- match.arg(component)
  if(component == "main"){
    fcall <- getCall(object)
    X = object$X
    if(is.null(fcall$formula)){
      stop("Model without formula not yet implemented. Please refit your model using the formula argument.")
      # fcall$formula = as.formula(object$formula)
    }
  }else{
    fcall <- getCall(object)
    X = object$lv.X
    object$X <- X
    if(is.null(fcall$lv.formula)){
      stop("Model without lv.formula not yet implemented. Please refit your model using the formula argument.")
    }
    fcall$formula <- fcall$lv.formula
  }
  
 fcall$formula <- as.formula(paste0("~0+species+",paste0("species:", labels(terms(as.formula(fcall$formula))), collapse = "+")))
 X = cbind(do.call("rbind",replicate(ncol(object$y),data.frame(X,check.rows = FALSE),simplify=FALSE)),species = as.factor(rep(colnames(object$y),each=nrow(object$X))))
 frm = model.frame(fcall$formula, X)
 
  emmeans:::recover_data.call(getCall(object), trms = terms(fcall$formula), na.action = NULL, data = X, params = "pi", frame = frm, pwts = NULL, addl.vars = NULL)
}


emm_basis.gllvm = function(object, trms, xlev, grid, component = c("main", "LV"), ...){
  # if(is.null(object$call$formula)){
  #   fcall <- getCall(object)
  #   fcall$formula = as.formula(object$formula)
  #   fcall$formula <- as.formula(paste0("~0+",paste0(labels(terms(as.formula(fcall$formula))),":species",collapse = "+")))
  #   trms = terms(as.formula(object$formula))
  #   }
  component = match.arg(component)
   xlev$species = colnames(object$y) # make sure species ordering is correctly retained

  if(component == "main"){
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
    X = model.matrix(trms, m)
    
    bhat = cbind(object$param$beta0,object$params$Xcoef)
    V = vcov(object) # per species organized
    V <- V[row.names(V)=="b", colnames(V)=="b"]
    #reorder V to per covariate
    V <- V[order(rep(1:ncol(bhat),times=ncol(object$y))),order(rep(1:ncol(bhat),times=ncol(object$y)))]
  }else if((object$num.lv.c+object$num.RR)>0 && component == "LV"){
    m = model.frame(trms, grid, na.action = na.pass, xlev = xlev)
    X = model.matrix(trms, m)
    
    bhat = cbind(object$params$beta0, object$params$theta[,1:(object$num.lv.c+object$num.RR),drop=FALSE]%*%t(object$params$LvXcoef))
    V = gllvm:::RRse(object, return.covb = TRUE) # per covariate organised
    covMat <- object$Hess$cov.mat.mod
    colnames(covMat) <- row.names(covMat) <- names(object$TMBfn$par[object$Hess$incl])
    Sb <- covMat[colnames(covMat)=="b",colnames(covMat)=="b"]
    if(!object$beta0com){
      bidx <- rep(c(TRUE,rep(FALSE,ncol(Sb)/ncol(object$y)-1)), ncol(object$y))
    }else{
      bidx <- c(TRUE,rep(FALSE, ncol(Sb)-1))
    }
    Sb <- Sb[bidx,bidx]
    
    # Add covariances for intercept.
    K = ncol(object$lv.X.design)
    d = object$num.lv.c+object$num.RR
    p = ncol(object$y)

    covMat <- covMat[colnames(covMat)%in%c("b","b_lv","lambda"),colnames(covMat)%in%c("b","b_lv","lambda")]
    covMat <- covMat[c(bidx,rep(TRUE,sum(colnames(covMat)%in%c("b_lv","lambda")))),c(bidx,rep(TRUE,sum(colnames(covMat)%in%c("b_lv","lambda"))))]
    #add first row and column of zeros before b_lv, for first species
    covMat <- rbind(covMat[1:(p+d*K),, drop=FALSE],0,covMat[-c(1:(p+d*K)),, drop=FALSE])
    covMat <- cbind(covMat[,1:(p+d*K), drop=FALSE],0,covMat[,-c(1:(p+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:(p+d*K+idx[q]-1),],0,covMat[(p+d*K+idx[q]):ncol(covMat),])
        covMat <- cbind(covMat[,1:(p+d*K+idx[q]-1)],0,covMat[,(p+d*K+idx[q]):ncol(covMat)])
      }
    }
    row.names(covMat)[row.names(covMat)==""]<-colnames(covMat)[colnames(covMat)==""]<-"lambda"
    
    covLb <- covMat[colnames(covMat) =="lambda",colnames(covMat)=="b",drop=FALSE]
    if(object$beta0com){
      covLb <- do.call(cbind,replicate(ncol(object$y),covLB, simplify=FALSE))
    }
    covLB <- covMat[colnames(covMat)=="lambda",colnames(covMat)=="b_lv", drop=FALSE]
    covb_lvb <- covMat[colnames(covMat) =="b_lv",colnames(covMat)=="b",drop=FALSE]
    if(object$beta0com){
      covb_lvb <- do.call(cbind,replicate(ncol(object$y),covb_lvb, simplify=FALSE))
    }
    covBb <- matrix(0,  K*p,ncol(object$y))

    for(k in 1:K){
        for(j in 1:p){
          for(j2 in 1:p){
            for(q in 1:d){
              covBb[j+p*(k-1),j2] <- covBb[j+p*(k-1),j2]+object$params$LvXcoef[k,d]*covLb[j+p*(q-1),j2]+object$params$theta[j,q]*covb_lvb[(q-1)*K+k,j2]
            }
        }
        }
    }
  V = rbind(cbind(Sb,t(covBb)),cbind(covBb,V))
  }else{
    stop("Invalid model.")
  }

  nbasis=matrix(NA)
  dfargs = list(df = Inf)
  dffun <- function(k, dfargs)Inf
  
  list(X=X, bhat=c(bhat), nbasis = nbasis, V=V,dffun = dffun, dfargs = dfargs)
}

emmeans::.emm_register("gllvm","gllvm")

library(gllvm)
data(spider)
X = data.frame(covariate = rep(c("one","two","three","four"),7), soil.dry = spider$x$soil.dry)
model <- gllvm(spider$abund,X,formula=~covariate+soil.dry,family="poisson",num.lv=0)

emmeans(model, ~covariate:species)

NOTE: A nesting structure was detected in the fitted model:
    covariate %in% species
 covariate species     emmean     SE df lower.CL upper.CL
 four      Alopacce   1.35571 0.1714 60    1.013    1.699
 one       Alopacce   1.87027 0.1601 60    1.550    2.191
 three     Alopacce   1.81912 0.1497 60    1.520    2.119
 two       Alopacce   0.10058 0.2930 60   -0.486    0.687
 four      Alopcune   1.50513 0.1750 60    1.155    1.855
 one       Alopcune   1.85791 0.1548 60    1.548    2.168
 three     Alopcune   1.21388 0.1993 60    0.815    1.612
 two       Alopcune   0.09861 0.3159 60   -0.533    0.730
 four      Alopfabr   0.33854 0.2788 60   -0.219    0.896
 one       Alopfabr   0.40582 0.3194 60   -0.233    1.045
 three     Alopfabr  -0.00134 0.3012 60   -0.604    0.601
 two       Alopfabr  -0.17177 0.3135 60   -0.799    0.455
 four      Arctlute  -2.77283 1.0772 60   -4.928   -0.618
 one       Arctlute  -2.45605 0.8689 60   -4.194   -0.718
 three     Arctlute  -0.55971 0.6116 60   -1.783    0.664
 two       Arctlute  -1.58347 0.7229 60   -3.029   -0.137
 four      Arctperi  -2.18649 0.7384 60   -3.663   -0.709
 one       Arctperi -11.92980 0.7384 60  -13.407  -10.453
 three     Arctperi  -1.84109 0.7037 60   -3.249   -0.434
 two       Arctperi  -1.03242 0.6552 60   -2.343    0.278
 four      Auloalbi   1.34367 0.1884 60    0.967    1.721
 one       Auloalbi   0.93478 0.2131 60    0.509    1.361
 three     Auloalbi   1.04166 0.2152 60    0.611    1.472
 two       Auloalbi   1.30856 0.1941 60    0.920    1.697
 four      Pardlugu   2.27998 0.1244 60    2.031    2.529
 one       Pardlugu   0.60111 0.2437 60    0.114    1.089
 three     Pardlugu   0.95790 0.2257 60    0.506    1.409
 two       Pardlugu  -0.07552 0.3473 60   -0.770    0.619
 four      Pardmont   2.78308 0.0942 60    2.595    2.972
 one       Pardmont   2.88896 0.0917 60    2.706    3.072
 three     Pardmont   3.12778 0.0796 60    2.969    3.287
 two       Pardmont   1.96368 0.1415 60    1.681    2.247
 four      Pardnigr   1.37521 0.1759 60    1.023    1.727
 one       Pardnigr   2.64373 0.1055 60    2.433    2.855
 three     Pardnigr   2.37218 0.1130 60    2.146    2.598
 two       Pardnigr   1.54579 0.1513 60    1.243    1.848
 four      Pardpull   2.55632 0.1018 60    2.353    2.760
 one       Pardpull   2.73715 0.0909 60    2.555    2.919
 three     Pardpull   2.81941 0.0907 60    2.638    3.001
 two       Pardpull   2.53371 0.1007 60    2.332    2.735
 four      Trocterr   3.11751 0.0770 60    2.964    3.272
 one       Trocterr   3.17553 0.0721 60    3.031    3.320
 three     Trocterr   3.26193 0.0722 60    3.117    3.406
 two       Trocterr   3.15966 0.0748 60    3.010    3.309
 four      Zoraspin   1.29279 0.1863 60    0.920    1.665
 one       Zoraspin   1.11741 0.1972 60    0.723    1.512
 three     Zoraspin   0.70183 0.2347 60    0.232    1.171
 two       Zoraspin   1.66955 0.1737 60    1.322    2.017

Confidence level used: 0.95 

Note to self: ordering of the SEs is a bit tricky, since the order of emmeans doesn't quite correspond with the order of the asymptotic covariance matrix in gllvm, so that needs to be kept an eye on.

This also works for models with reduced rank terms, but at present this requires installing the package from the development repo.

[BertvanderVeen]

Yup, all good.

[gerverska]

Thanks, Bert. It looked like dfargs = list(df = Inf) does this, so I was a little confused with the df = 60 emmGrid output above. Apologies for my ignorance, but what does setting dfargs = list(df = Inf) achieve?

[BertvanderVeen]

Thanks, Bert. It looked like dfargs = list(df = Inf) does this, so I was a little confused with the df = 60 emmGrid output above. Apologies for my ignorance, but what does setting dfargs = list(df = Inf) achieve?

1. it gets emmeans to recognize we are looking at confidence intervals based on a large sample approximation of the sampling distribution
2. it keeps emmeans from letting the user do t-tests of the parameters while not appropriate.

[BertvanderVeen]

Minor correction, species ordering was changed by emmeans so I have hard-coded a fix for this.

[gerverska]

Thanks, Bert!

[AlainZuur]

Here is my small contribution.

I have this model:
M2.X <- gllvm(y = SpecData,
X = CovData,
formula = ~ Location,
row.eff = "random",
num.lv = 2,
seed = 12345,
control.start = list(n.init = 5, jitter.var = 0.1),
family = "gamma")

Location is a categorical covariate with 8 levels (sorry...can't help). SpecData contains 16 morphometric variables. 240-ish rows (bees).

The output from gllvm gives the estimated intercepts:
b0 <- M2.X$params$beta0

And the slopes...which are corrections of the intercept in this specific example. You can plot them with:
coefplot(M2.X,
which.Xcoef = c("Location2", "Location3",
"Location4", "Location5",
"Location6", "Location7",
"Location8"),
order = FALSE,
cex.ylab = 0.7,
mar = c(5, 5, 1, 1),
mfrow = c(2, 4))

Or you can grab them from:
M2.X$params$Xcoef

So far so good. I then wanted to calculate the intercept + the correction of the intercept for each of the 8 locations...and for each morphometric variable. Hence...in this model:

mu_ij = exp(Intercept + Location *beta_j + row-specific random effect + LV * theta_j)

I want the combined effect of: Intercept + Location *beta_j
And ..also the corresponding SE...which allows me to calculate a 95% CI.
That would tell me the Location effect for each location and each morphometric variable.

I think it would be quite tricky to write code that can do this in any scenario (that is...add more covariates etc). So..I set myself on a mission to write an exercise that explains this step-by-step. And once you understand how that works, you can easily adjust it for your own data sets (instead of relying on emmeans and friends).

So..in a univariate lm, it simplifies to:

Bees2$fLocation <- factor(Bees2$Location)
Ma <- lm(Forewing_length ~ fLocation, data = Bees2)

And then you get things like this:

summary(Ma)$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) 9.18853333 0.05210054 176.361573 3.791235e-237
fLocation2 0.18266667 0.07368129 2.479146 1.392786e-02
fLocation3 -0.53356667 0.07368129 -7.241549 7.526368e-12
fLocation4 -0.09946667 0.07368129 -1.349958 1.784294e-01
fLocation5 -0.21596667 0.07368129 -2.931092 3.737604e-03
fLocation6 0.20896667 0.08834077 2.365461 1.888456e-02
fLocation7 0.23130000 0.07368129 3.139196 1.928569e-03
fLocation8 -0.29396667 0.07368129 -3.989706 9.034554e-05

The intercept is the fitted value for location 1.
Intercept + fLocation2 is the fitted value for location 2.
Intercept + fLocation3 is the fitted value for location 3.
Intercept + fLocation4 is the fitted value for location 4.
etc etc

But that does not give you SEs.

So..we need a little trick. First, we need these two things:
Betas <- coef(Ma)
CovM <- vcov(Ma)

The estimated betas and their covariance matrix. Instead of calculating 9.18853333 + 0.18266667 by hand, I will do it via matrix multiplication.

Define C as follows
C <- matrix(nrow = 1, ncol = 8, c(1, 1, 0, 0, 0, 0, 0, 0))

Multiply this matrix with the estimates betas:
FitLoc2 <- C %*% Betas

So that is the intercept plus the correction of the intercept for location 2. The fitted value for Location 2.

Now work towards the SE. There is this stats rule:
var(a + b) = var(a) + var(b) + 2 * cov(a,b)

In other words, to get the SE of 'Intercept + fLocation2, we need:

var(Intercept) + var(fLocation2) + 2 * cov(Intercept, fLocation2)

where:
#' The var(Intercept) is the 1,1 element of CovM.
#' The var(fLocation2) is the 2,2 element of CovM.
#' The cov(Intercept,fLocation2) is the 1,2 element of CovM.

#' We can also do it a little bit more sophisticated, and use
#' the following rule: If mu = C x Betas, then the variance of
#' mu is equal to C x var(Betas) x t(C). The t() stands for 'transpose'.

#' Utilising this, we get:
Var.12 <- C %% CovM %% t(C)
SE.12 <- sqrt ( Var.12 )

In other words, to get the fitted value for Location 2 and its SE, we need to define the C as:
c(1, 1, 0, 0, 0, 0, 0, 0)

And then it is two lines of code:
Fit.12 <- C %% Betas
SE.12 <- sqrt ( C %% CovM %*% t(C) )

This gives you the 'Intercept + fLocation2' value and also its standard error. What remains is to do this for:
Intercept + fLocation3
Intercept + fLocation4
Intercept + fLocation5
Intercept + fLocation6
etc.

But that is just a matter of changing the C 6 more times. And this is all for a linear regression model.

#' What did we learn?
#' 1. If you want to calculate the sum of two regression parameters,
#' then you define C and calculate C x beta.
#' 2. If you want its standard error, you take the square root of
#' C x cov(betas) x t(C)

Now..the good news is...in GLLVM it is a monkey-see-monkey copy.

Here are all the parameters in the model:
Betas <- M2.X$TMBfn$par[M2.X$Hess$incl]

And the covariance matrix:
CovM <- M2.X$Hess$cov.mat.mod

Now we need to do exactly the same thing...but then 16 times. And you get a headache which values of Betas to add up.

This is me doing it in a loop (j is for the 16 morphometric variables).

#' And now for all morphometric variables.
#' Initialise the matrix
LocationEffect <- matrix(nrow = 16, ncol = 8)

#' Pre-create C outside the loop and modify it in the loop
C <- matrix(0, nrow = 1, ncol = 8)

#' Iterate over morphometric variables j.
for (j in 1:16){
for (k in 1:8){
C[,] <- 0 #' Reset C
C[, k] <- 1 #' Set the kth position to 1
C[, 1] <- 1 #' Ensure the first element is always 1

#' Calculate Intercept + correction of the intercept via C x Beta
offset <- 8 * (j - 1)
LocationEffect[j, k] <- C %*% Betas[1:8 + offset]

}
}
#' Results:
rownames(LocationEffect) <- rownames(BetasMatrix)
colnames(LocationEffect) <- colnames(BetasMatrix)
LocationEffect

This is really the same thing as 16 times what I did for the univariate lm.

And the standard errors:
#' Initialise the matrices
LocationEffect <- matrix(nrow = 16, ncol = 8)
SELocationEffect <- matrix(nrow = 16, ncol = 8)

#' Pre-create C outside the loop and modify it in the loop
C <- matrix(0, nrow = 1, ncol = 8)

#' Iterate over morphometric variables j.
for (j in 1:16){
for (k in 1:8){
C[,] <- 0 #' Reset C
C[, k] <- 1 #' Set the kth position to 1
C[, 1] <- 1 #' Ensure the first element is always 1

#' Calculate Intercept + correction of the intercept via C x Beta[..]
offset <- 8 * (j - 1)
LocationEffect[j, k] <- C %*% Betas[1:8 + offset]

# Calculate C x CovM x t(C), and take the square root 
VAR <- C %*% CovM[1:8 + offset, 1:8 + offset] %*% t(C)
SELocationEffect[j, k] <- sqrt(VAR)

}
}

I'm sure that it can be done more efficiently. The end result is this graph:

This shows the 'Intercept + Location" effect for each morphometric variable. Keep in mind that there are other terms in the model, and there is also a link function.

For this specific example, the 'Intercept + Location' should probably be the same as predict(.., type = "link", and without the random stuff).

Post-hoc tests
The next question is then....which of those locations are statistically different from each other?

In a univariate case, you would do this:
Ma <- lm(Forewing_length ~ fLocation, data = Bees2)
TukeyTest <- glht(Ma, linfct = mcp(fLocation = "Tukey")) #' Or whatever test you want to use.
summary(TukeyTest)

#' And you can plot the results.
par(mar = c(5,5,2,2))
plot(TukeyTest)

All what that is doing is defining your C matrix like this:
C <- matrix(nrow = 1, ncol = 8,
c(0, 0, 0, 0, 0, 0, 1, -1))

This is testing whether Location 7 is different from 8.
Here is the difference (between fLocation7 and fLocation8) and the SE of the difference:

#' Calculate the C x B:
Diff <- as.numeric(C %*% Betas)

#' Calculate the standard error of Diff:
SE.Diff <- as.numeric(sqrt(C %% CovM %% t(C)))

It gives exactly the same value as glht. Then all that is needed is to do this for every combi of locations. And correcting for multiple tests.

And this translates directly to GLLVM. But because of the 16 morphometric variables, the number of tests is quite large.

And this is where I am right now. So..my weekend plan is to trim some trees and extend this to the GLLVM code.

I hope this helps the original question from @jebyrnes ? And hopefully, this is actually correct...:-).

Alain