Combined standard errors for subsets of coefficients?

[gerverska]

Hello @BertvanderVeen and @JenniNiku!

Following the discussion for #113, it seems based on Gordana's overview that if

$log\left(\lambda_0\right) = \beta_0 + \beta_1x_0$
and
$\lambda_1 = \lambda_0exp\left(\beta_1\right)$

then the following should be true:
$\lambda_1 = \lambda_0exp\left(\beta_1 + \beta_2 + \beta_3\right)$
assuming a base model of
$log\left(\lambda_0\right) = \beta_0 + \beta_1x_0 + \beta_2x_0 + \beta_3x_0$

So if coefficients are combined like this (with each coefficient having its own standard error) to estimate the total multiplicative effect on the mean, is there then a cumulative standard error for the coefficients involved in this combination? If so, how is this calculated?

I've already combined the estimates of the term:species coefficients output from summary.gllvm() and exponentiated them as above to get the total multiplicative effect experienced by each mean species response, but I was just wondering if it would be helpful to display the total SE for this multiplicative effect (if it exists).

[BertvanderVeen]

Thanks Kyle.

Yes, you can determine standard errors for combinations of parameter estimates as described above, but it requires taking the covariance of those estimators into consideration. This will require two things: a matrix (let's call it) $\textbf{C}$ that forms the appropriate linear combinations, and the asymptotic covariance matrix of the model, retrieved via vcov(model).

Then, it should be something similar as in the following example:

library(gllvm)
data(spider)

# Make example covariate with 4 categories
X = data.frame(covariate = rep(c("one","two","three","four"),7))

# Fit example model
model <- gllvm(spider$abund,X,family="poisson",num.lv=0)

# Estimated coefficients
bs <- cbind("Reference/Intercept" = model$params$beta0, model$params$Xcoef)

# Get asymptotic covariance matrix
S <- vcov(model)
S <- vcov(model)[colnames(S)=="b",colnames(S)==c("b")] # intercepts and Xcoef are stored as "b"

# Matrix that forms the combinations
C <- cbind(c(1,0,0,0),c(1,1,0,0),c(1,0,1,0),c(1,0,0,1))

# For clarity name the intercepts
colnames(S)[1:ncol(model$y)]<-row.names(S)[1:ncol(model$y)] <- "Reference/Intercept"

# Form the combinations
coefs <- bs%*%C

# Calculate SEs for the combinations
ses <- matrix(sqrt(diag(t(Matrix::bdiag(replicate(ncol(model$y),C,simplify=FALSE)))%*%S%*%Matrix::bdiag(replicate(ncol(model$y),C,simplify=FALSE)))), ncol = ncol(bs), byrow = TRUE)

# Species names on the rows, column names are case dependent..
row.names(coefs) <- row.names(ses) <- row.names(bs)

The above example is for categorical covariates only, but the principle for terms with a predictor variable is the same. For example, the SE for the intercept, and for the sum of an intercept and a term with covariate would require

$$ \textbf{C} = \left( \begin{array}{ccc} 1 & 1 \\ 0 & x_i\\ \end{array} \right) $$

or just the last column, if you only want to sum over both.

[gerverska]

Too cool--thanks for the code and explanation, Bert. I'll run your example and adapt the code to my setup once I get a better feel for what's happening.

From here, is it as straightforward as it seems to get a Z score and then a p value?
$Z = \left(\beta_1 + \beta_2 + \beta_3\right) / SE$
And then pnorm(z, lower.tail = FALSE)?

Or am I missing something here?

[BertvanderVeen]

That depends on the test you are looking to execute. For a two-tailed z-test: 2 * pnorm(-abs(Z)) with indeed Z = par/se.

[gerverska]

Hi Bert,

Apologies for not testing your code out earlier, but the above code gives the following error for the CRAN version:

# Calculate SEs for the combinations
> ses <- matrix(sqrt(diag(t(Matrix::bdiag(replicate(ncol(model$y),C,simplify=FALSE)))%*%S%*%Matrix::bdiag(replicate(ncol(model$y),C,simplify=FALSE)))), ncol = ncol(bs), byrow = TRUE)
Error in t.default(Matrix::bdiag(replicate(ncol(model$y), C, simplify = FALSE))) : 
  argument is not a matrix

The output of
Matrix::bdiag(replicate(ncol(model$y), C, simplify=FALSE) appears to be a 48 x 48 sparse Matrix of class "dtCMatrix", which seems to be where the t() error is coming from. This was fixed by converting the Matrix::bdiag(replicate(ncol(model$y), C, simplify=FALSE) output to a matrix (see the as.matrix() calls below:

# Calculate SEs for the combinations
ses <- matrix(sqrt(diag(t(Matrix::bdiag(replicate(ncol(model$y),
                                                  C,
                                                  simplify=FALSE)) |> 
                              as.matrix()) %*% S %*% 
                            Matrix::bdiag(replicate(ncol(model$y),
                                                    C,
                                                    simplify=FALSE)) |> 
                            as.matrix())),
              ncol = ncol(bs), 
              byrow = TRUE)

Here are the following standard errors I got:

> ses
               [,1]       [,2]       [,3]       [,4]
Alopacce 0.15075568 0.16012816 0.11322770 0.27735011
Alopcune 0.16903086 0.11470787 0.18569533 0.30151131
Alopfabr 0.17149858 0.31622776 0.18257419 0.20851442
Arctlute 0.99999914 0.50000007 0.26726126 0.37796449
Arctperi 0.40824829 0.40824832 0.30151134 0.21320072
Auloalbi 0.17960528 0.17407767 0.19611613 0.15811389
Pardlugu 0.11547005 0.21821789 0.21320073 0.33333331
Pardmont 0.09365858 0.08980265 0.07881104 0.14142136
Pardnigr 0.17149859 0.06950480 0.09622505 0.13245324
Pardpull 0.09853293 0.07161149 0.08137885 0.08671100
Trocterr 0.07432941 0.05735393 0.06509446 0.06324555
Zoraspin 0.17149860 0.13736056 0.20851441 0.11547005

Is this the output that I should get?

[BertvanderVeen]

Yes, that is correct.