You are incorrect in stating that "The only difference between the two models is that the Parm3 was removed for model 2". You are overlooking the fact that
glm.nb() is also estimating the $\theta$ parameter of the Negative Binomial model.
Here is an example from
quine.nb1 <- glm.nb(Days ~ Sex/(Age + Eth*Lrn), data = quine)
quine.nb2 <- update(quine.nb1, . ~ . + Sex:Age:Lrn)
In this setting,
quine.nb2 is your model 1 and
quine.nb1 is your model 2, and the null deviances do indeed differ. However, if we fit the model of
quine.nb1 using the estimated value of $\theta$ from
quine.nb2 we should see the same null deviances. Here I refit using
glm() and the
negative.binomial() family function so I am certain the same code is used for fitting and I can fix $\theta$ at a "known" value.
theta <- quine.nb2$theta
f1 <- formula(quine.nb1)
f2 <- formula(quine.nb2)
m1 <- glm(f1, data = quine, family = negative.binomial(theta = theta))
m2 <- glm(f2, data = quine, family = negative.binomial(theta = theta))
Now we extract the Null deviance from the two models
> m1$null.deviance # $ ignore: the code & mathjax is messing up again...
and they are the same.
The key point is that with
glm.nb() the reported null deviance is conditional upon the estimated value of $\theta$. You have to be careful that you are comparing like with like.