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)

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:


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 

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)



1 Answer 1


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
  (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),
   (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

*see here for more info.

**see here for more info.

  • $\begingroup$ Thanks Cyrus Mohammadian. If setting the contrasts to list(target = "contr.sum") is the way to use the grand mean as reference any idea how come releveling input to exclude B instead of A gives different estimates for C-E? $\endgroup$
    – dan
    Sep 10, 2016 at 19:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.