Você pode obter as estimativas de cada casela de diversas formas: 1) removendo intercepto da fórmula, 2) mudando o tipo de constraste, 3) calculando as funções lineares por fora com matrizes apropriadas e 4) fazendo o item 3 com a popMeans() e/ou popMatrix(), que considero a opção mais simples. Veja o CMR.

##-----------------------------------------------------------------------------

## dados artificiais
da <- expand.grid(A=gl(4,4), B=gl(3,1))
da$y <- rpois(nrow(da), lambda=10)

## constrates usados por padrão
options()$contrasts

## removendo o intercepto, primeiro fator com estimativas por nível
m0 <- glm(y~-1+A+B, data=da, family=poisson)
summary(m0)
m0 <- glm(y~-1+B+A, data=da, family=poisson)
summary(m0)

## mudando a opção de contraste localmente
m0 <- glm(y~A+B, data=da, family=poisson,
contrasts=list(A=contr.helmert, B=contr.SAS))
summary(m0)

M <- unique(model.matrix(m0))
M

## comportamento padrão
m0 <- glm(y~A+B, data=da, family=poisson)
summary(m0)

M <- unique(model.matrix(m0))
M

##-----------------------------------------------------------------------------

## estimativas para cada combinação A*B
M%*%coef(m0)

## para soma direta de matrizes
source("http://fisher.osu.edu/~schroeder_9/AMIS900/blockdiag.R")

mA <- contrasts(da$A)
mB <- contrasts(da$B)
M <- list(mA, mB)
M <- cbind(1, blockdiag(M))
M <- unique(M)
M

M%*%coef(m0)
coef(m0)

require(doBy)

## efeito marginal de A (ou seja, na média dos níveis de B)
popMeans(m0, effect="A")

## efeito marginal de B (ou seja, na média dos níveis de A)
popMeans(m0, effect="B")

## efeito de A condicional à B=1
popMeans(m0, effect="A", at=list(B="1"))

## aqui é para ver as matrizes envolvidas
popMatrix(m0, effect="A")
popMatrix(m0, effect="A", at=list(B="1"))

##-----------------------------------------------------------------------------

À disposição.
Walmes.