# Obtaining all effects in a multinomial regression where the grand mean is set as baseline

I have data from these set of experiments:

In each experiment I infect a neuron with a rabies virus. The virus climbs backwards across the dendrites of the infected neuron and jumps across the synapse to the input axons of that neuron. In the input neurons the rabies will then express a marker gene thereby labeling them. This allows me to see which neurons are inputs to the target neuron I infected and thus create a connectivity map of a certain region in the brain.

In each experiment I obtain counts of all the infected input neurons of the target neuron I infected.

Here's a simulation of the data: (3 targets and 5 inputs)

set.seed(1)
probs <- list(c(0.4,0.1,0.1,0.2,0.2),c(0.1,0.3,0.4,0.1,0.1),c(0.1,0.1,0.4,0.2,0.2))
mat <- matrix(unlist(lapply(probs,function(p) rmultinom(1, as.integer(runif(1,50,150)), p))),ncol=3)
inputs <- LETTERS[1:5]
targets <- letters[1:3]
df <- data.frame(input = c(unlist(apply(mat,2,function(x) rep(inputs ,x)))),target = rep(targets ,apply(mat,2,sum)))


What I'd like to estimate is the effect of each target neuron on these counts, relative to the grand mean. I was thinking that a multinomial regression model is appropriate in this case, where the contrasts are set to the contr.sum option:

library(foreign)
library(nnet)
library(reshape2)

df$input <- factor(df$input,levels=inputs)
df$target <- factor(df$target,levels=targets)
fit <- multinom(input ~ target, data = df,contrasts = list(target = "contr.sum"))
# weights:  20 (12 variable)
initial  value 505.363505
iter  10 value 445.057386
final  value 441.645283
converged


Which gives me:

> summary(fit)$coefficients (Intercept) target1 target2 B 0.08556288 -1.743854 1.6062660 C 0.55375003 -2.094266 1.2616939 D -0.17624590 -1.364270 0.6284231 E -0.04091248 -1.617374 0.6601274  So the effects for input A are not reported and I would like to obtain both the effects of all targets on all inputs. I'm wondering if adding a mean across targets and a mean across inputs, and setting them as baseline dummy variables is a good solution: #add target mean mat <- cbind(mat,round(apply(mat,1,mean))) colnames(mat)[ncol(mat)] <- "x" targets <- c(targets,"x") #add input mean mat <- rbind(mat,round(apply(mat,2,mean))) rownames(mat)[nrow(mat)] <- "X" inputs <- c(inputs,"X")  So x and X represent the means of targets and inputs, respectively, and are rounded so that they are counts. df <- data.frame(input = c(unlist(apply(mat,2,function(x) rep(inputs ,x)))),target = rep(targets ,apply(mat,2,sum))) df$input <- factor(df$input,levels=rev(inputs)) df$target <- factor(df$target,levels=rev(targets))  And then fit the multinom regression using dummy coding: fit <- multinom(input ~ target, data = df)  Thanks ## 1 Answer You have a couple different questions in your post -I'll address them separately. 1) You can obtain A by using relevel() to add a new reference level: df <- within(df, input <- relevel(input, ref = 2)) #Reference is now B summary(fit)$coefficients
(Intercept)    target1    target2
A -0.08551218  1.7437515 -1.6061592 #Notice A and not B
C  0.46819646 -0.3503939 -0.3445671
D -0.26180114  0.3796378 -0.9779083
E -0.12649098  0.1265176 -0.9461571


2) To obtain p-values using the Wald test (z-statistic) try the following*:

z.stat <- summary(fit)$coefficients/summary(fit)$standard.errors
(1 - pnorm(abs(z.stat), 0, 1)) * 2
(Intercept)    target1    target2
A   0.8245699 0.001986126 0.00330052
C   0.3352785 0.512106514 0.54471123
D   0.5900357 0.523038723 0.09852789
E   0.8002826 0.831427999 0.10975537


or using a t-statistic**:

fit.sum <- summary(fit)
pt(abs(fit.sum$coefficients / fit.sum$standard.errors),
df=nrow(df)-1,lower=FALSE)
(Intercept)  target1   target2
A   0.5876428 0.9989175 0.9982271
C   0.8319889 0.7437060 0.7274249
D   0.7047907 0.7382471 0.9502347
E   0.5997757 0.5842167 0.9446179