Yes, a deviance test is still valid. Some more details: Since the general theory is not specific for binomial models, I will start out with some general theory, but use binomial examples (and R.)
GLM's is based on the exponential dispersion model
$$
f(y_i;\theta_i,\phi)= \exp\left\{ w_i [y_i \theta_i -\gamma(\theta_i)]/\phi +\tau(y_i,\phi/w_i)\right\}
$$
where $y_i$ is the observation, $\theta_i$ parameter which depends on a linear predictor $\eta_i=x^T\beta$, $\phi$ a scale parameter and $w_i$ a prior weight. To understand the notation, think about a normal theory model, which glm's generalize. There $\phi$ is the variance $\sigma^2$, and if $y_i$ is the mean of a group of $w_i$ independent observations with the same covariables, then the variance is $\phi/w_i$. The last term $\tau(y_i,\phi/w_i)$ is often of little interest since it does not depend on the interest parameters $\theta_i$ (or $\beta$,) so we will treat it cavalierly.
So the binomial case. If we have an observation $y_i^* \sim \mathcal{Binom}(w_i,p_i)$ then we will treat $y_i=y_i^*/w_i$ as the observation, so that the expectation of $y_i$ is $p_i$ and its variance $\frac{p_i(1-p_i)}{w_i}$. The binomial pmf can then be written as
$$
f(y_i;\theta_i,\phi)=\exp\left\{ w_i[y_i\theta_i-\log(1+e^{\theta_i})]/\phi + \log\binom{w_i/\phi}{y_i w_i/\phi} \right\}
$$ where $\phi=1$ and $\theta_i=\log\frac{p_i}{1-p_i}$. We can identify $\gamma(\theta_i)=\log(1+e^{\theta_i})$ and $\tau(y_i,\phi/w_i) = \log\binom{w_i/\phi}{y_i w_i/\phi} $. This form is chosen such that we can get the quasi-model just by allowing $\phi>0$ to vary freely.
The quasi-likelihood we then get from this model, is constructed *to function as a likelihood for the $\theta_i$ (or $\beta$) parameters, it will not work as a likelihood for $\phi$. This means that the quasi-likelihood function shares enough of the properties of a true likelihood function that the usual likelihood asymptotic theory still goes through, see also Idea and intuition behind quasi maximum likelihood estimation (QMLE). Since it does not have this properties as a function of $\phi$, inference about $\phi$ must be treated outside that framework. Specifically, there is no reason to hope that maximizing the qlikelihood in $\phi$ to give good results.
Now, the analysis of deviance. We define the saturated model S by giving each observation its own parameter, so setting $\hat{\mu}_i=\gamma'(\hat{\theta}_i)=y_i$. Then by assuming for the moment that $\phi=1$ we get
$$
D_M=2\sum_i \left\{ w_i[( y_i \theta(y_i)-\gamma(\theta(y_i)))-( y_i\hat{\theta}_i-\gamma(\hat{\theta_i
}) ) ]\right\}
$$
which is twice the loglikelihood-ratio for testing the reduced model M within the saturated model S. Note that this does not depend on the function $\tau$ at all. For the case of normal-theory models, this is the residual sum of squares (RSS), which is not a function of the scale parameter $\phi=\sigma^2$ either.
$D_M/\phi$ is the scaled deviance while $D_M$ often is called the residual deviance, since in normal models it corresponds to the RSS. In normal models we have $D_M/\phi \sim \chi^2_{n-p}$ so an unbiased estimator of the variance parameter $\phi$ in this case is $\hat{\phi}=D_M/(n-p)$ and this might hold as an approximation also in other cases, but often better is
$$
\tilde{\phi}=\frac1{n-p}\sum_i \frac{(y_i-\hat{\mu}_i)^2}{V(\hat{\mu_i})/w_i}
$$ where $V$ is the variance function, in the binomial case $V(\mu)=\mu(1-\mu)$. In the binomial case, this is considered to be better, and is the scale estimate used by R.
If we are interested in a submodel $M_0 \subset M$, with $q < p$ regression parameters, then the likelihood ratio test is
$$
\frac{D_{M_0}-D_M}{\phi} \stackrel{\text{approx}}{\sim} \chi^2_{p-q}
$$
and with estimated scale we might use
$$
\frac{D_{M_0}-D_M}{\hat{\phi}(p-q)} \stackrel{\text{approx}}{\sim} \mathcal{F}_{p-q,n-p}
$$
in analogy with the normal theory.
So, let us look at a simulated example.
set.seed(7*11*13) # My public seed
n <- 200
k <- 5
N <- n*k
intercept <- rnorm(n, 0, 1)
x <- rnorm(n, 1, 1.5)
beta <- 0.1
expit <- function(x) 1/(1+exp(-x))
eta <- intercept + beta*x
p <- expit(eta)
Y <- rbinom(n, k, p)
This creates overdispersion by simulating a random intercept for each of the $n=200$ groups of size $k=5$. Then we will estimate a simple model two ways, by using a binomial likelihood, and then a quasibinomial likelihood:
mod0 <- glm( cbind(Y, k-Y) ~ x, family=binomial)
modq <- glm( cbind(Y, k-Y) ~ x, family=quasibinomial)
Then the model summaries:
summary(mod0)
Call:
glm(formula = cbind(Y, k - Y) ~ x, family = binomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.053 -1.180 -0.103 1.180 2.836
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.02787 0.07632 -0.365 0.71496
x 0.12941 0.04170 3.103 0.00192 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 441.41 on 199 degrees of freedom
Residual deviance: 431.62 on 198 degrees of freedom
AIC: 749.1
Number of Fisher Scoring iterations: 3
> summary(modq)
Call:
glm(formula = cbind(Y, k - Y) ~ x, family = quasibinomial)
Deviance Residuals:
Min 1Q Median 3Q Max
-3.053 -1.180 -0.103 1.180 2.836
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.02787 0.10117 -0.275 0.7832
x 0.12941 0.05529 2.341 0.0202 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for quasibinomial family taken to be 1.757479)
Null deviance: 441.41 on 199 degrees of freedom
Residual deviance: 431.62 on 198 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 3
Compare the two summaries. They are very similar, the differences are in coefficient standard errors, and the printed scale parameter estimate, and lacking AIC of the modq summary. Check that you can calculate, "by hand", the standard errors for the quasimodel modq from the standard errors for mod0 and the estimated scale.
The printed deviances, and deviance residuals, are identical. This is because the residual deviance is defined by taking $\phi=1$ in both cases. The null deviance is the residual deviance for the null model, the model with only an intercept. The scaled deviance is not printed, but can be calculated from the output.
The analysis of deviance is calculated by the anova() function. Here we will see differences. First the model based on a binomial likelihood:
anova(mod0, test="Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: cbind(Y, k - Y)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 199 441.41
x 1 9.7883 198 431.62 0.001756 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
so here the regression seems significant. Then for the quasi-model:
anova(modq, test="F")
Analysis of Deviance Table
Model: quasibinomial, link: logit
Response: cbind(Y, k - Y)
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev F Pr(>F)
NULL 199 441.41
x 1 9.7883 198 431.62 5.5695 0.01925 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
What is printed as F here is (in this case) the scaled deviance (since $p-q=1$.)
(I will come back to your second question)
