The reason behind the trick of assuming $p$ and $q$ are independent, differentiating, then applying the relation. We want to simplify this expression:
$$\bar{n}=\sum_{n=0}^N W(n)n=\sum_{n=0}^N\frac{N!}{n!(N-n)!}p^{n}q^{N-n}n$$
where $q=1-p$.
The trick used is to write:
$$np^{n}=p\frac{\partial}{\partial p}(p^{n})$$
plugging in,
$$\bar{n}=\sum_{n=0}^N\frac{N!}{n!(N-n)!}\left[p\frac{\partial}{\partial p}(p^{n})\right]q^{N-n}$$
But the book does something strange here.
$$\implies \bar{n}=p\frac{\partial}{\partial p}\left[\sum_{n=0}^N\frac{N!}{n!(N-n)!}p^{n}q^{N-n}\right]$$
the, according to binomial theorem,
$$\bar{n}=p\frac{\partial}{\partial p}(p+q)^N$$
$$\bar{n}=pN(p+q)^{N-1}$$
Now, if we use $q=1-p$,
$$\bar{n}=pN$$

The question: Why did we consider $p$ and $q$ independent variables and put the partial derivative out of the sum, then differentiate and at last apply the relation between $p$ and $q$?
Why are we allowed to do so? (I don't want to know why we decided to do so (then one would simply say "because otherwise we couldn't solve it!"), but rather what makes us enable to do so.)
The book: Fundamentals of Statistical and Thermal Physics, Frederick Reif, First chapter
 A: This is a case where knowing less is better. The trick that is implied is that rather than the object of interest, one calculates something more general.
So if you would ignore the fact that it is a problem about probabilities, we have a generic expression for $\bar{n}$ which is a function in two independent variables $p$ and $q$
$$
\bar{n}(p,q) = \sum_{n=0}^N \binom{N}{n} n p^n q^{N-n}
$$
One always has $p \frac{\partial}{\partial p} p^n = n p^n$ and since the binomial does not depend on $p$ and neither does $q$ in this generic expression, we can bring the "operator" $p \frac{\partial}{\partial p}$ outside the summation because it is finite sum and hence no problems would arise from lack of convergency. With the various other steps, one shows then that
$$
\bar{n}(p,q) = p N (p+q)^{N-1}
$$
You know it is correct because you can also follow the various steps in the opposite direction. Alternatively, you could expand the last expression in to a sum and rearrange terms a bit and change the summation index by one to obtain the original expression as well.
This means that this is a correct relation for independent variables $p$ and $q$ which can have any value, for instance $(p,q)=(2,5)$ or $(\sqrt{2},\pi)$ or any other combination of numbers. 
Now it happens to be the case that the general expression we started with has a particular extra meaning in that it corresponds to the average value of $n$.  But the relation we derived is valid for for any $p$ and $q$, in particular the special case $q=1-p$.
It is crucial in the derivation that the partial derivative is introduced when $p$ and $q$ are independent. Hence it should also disappear again before one inserts the dependency $q=1-p$.
This approach is used often to remove factors that depend on the summation index and reduce summation series to simple binomials and most people find it easier to use than the alternative but equivalent method of rearranging terms and shifting summation indices.
A: I understand that the joggling with variables that at first
are dependent, then
treated as independent and finally again becomes
dependent is somewhat confusing. Although mathematically
correct, there are cleaner ways of presenting it.
The presented calculations involve in reality the following 4 expressions:


*

*$E_1(p) =  \sum_{n=0}^N W(n) n = 
   \sum_{n=0}^N 
   \frac{N!}{n!(N-n)!}
    p^n (1-p)^{N-n} n$

*$E_2(p,q) =
   \sum_{n=0}^N 
   \frac{N!}{n!(N-n)!}
    p^n q^{N-n} n$

*$E_3(p,q)= N p (p+q)^{N-1}$

*$E_4(p) = N p$
Here, $E_2(p,q)$ and $E_3(p,q)$ are functions of two independent
variables $p$ and $q$, whereas $E_1(p)=E_2(p,1-p)$ and 
$E_4(p)=E_3(p,1-p)$
are functions of the variable $p$, only.
 What we want to show is that
$E_1(p)=E_4(p)$ or equivalently that $E_2(p,1-p)=E_3(p,1-p)$.
The trick using partial derivatives (which you seem to be ok with) shows
that for all real $p$ and $q$ you have
 $$ E_2(p,q)=E_3(p,q). $$
Now, this is an identity between two expressions in the two 
variables $p$ and $q$. 
It remains valid if you substitute whatever real-valued 
expression for $p$ and $q$
that you may think of.
 For example, expanding $E_2(pt, p+t)=E_3(pt,p+t)$ would yield
a (not very useful) identity in variables $p$ and $t$. 
Of more interest we have $E_3(p,1-p)=E_4(p,1-p)$.
which immediately gives
the wanted identity between $E_1(p)=E_4(p)$.
A possibly cleaner way to do the above calculation would be to
introduce nice independent variables $x$ and $y$ and establish 
as before that $E_2(x,y)=E_3(x,y)$ or
$$ 
\sum_{n=0}^N 
   \frac{N!}{n!(N-n)!}
   x^n y^{N-n} n = x N (x+y)^N
   $$
is valid for all real $x$ and $y$. In particular, if $p$ is any
constant, the formula is valid when you substitute
$x=p$ and $y=1-p$, leading to
the desired result and avoiding the
joggling between independent and dependent variables.
A: Very tricky for me to answer clearly, I will try as an exercise. As the late R.Feynmann once said, “a bird does not need to know aerodynamics to fly”, the general truth of which applies also to toddling toddlers such as myself. 
So, if the purpose it to calculate a partial derivative of a function $f(x,y)$, where $x$ and $y$ are constrained by an expression such as $g(x,y) = 0$, one cannot of course assume one variable is free and then substitute.
This is easily verified by an example $$ z = f(x,y) = x^2 \cdot y $$ subject to the constraint $x + y = 1$.
In this case, consideration of the constraint lead to elimination of $y$, and $$z = x^2 \cdot(1-x)$$ and $$\frac{\mathrm{d}z} {\mathrm{d}x } = 2x – 3x^2$$
Unconstrainted partial derivation yields $$\frac{\partial z} {\partial x } = 2x y$$ a different result. 
This is not surprising: in the first case, the variable $x$ varies along the line $x + y = 1$: in the second, not so. The operation of partial derivation critically depends on whether two variables are constrained, on which one is held fixed, etc.
But the purpose in the example you shared is different than calculating the partial derivative of a constrained function. 
Consider the expression $$\frac{\partial}  {\partial p }  (p^n) q$$
All is wanted is to find an equivalent expression, for further manipulation. 
One then considers the case that $q$ is a constant. 
The next step, i.e. writing $$\frac{\partial}  {\partial p }  (p^n) y  = \frac{\partial}  {\partial p }  [(p^n) y]$$ is certainly valid IF it is assumed $q$ is a constant, as well as applying the binomial theorem and subsequent derivation. 
Then you get to a step when $q$ being a constant or not does not matter anymore, i.e. the validity of the pointwise evaluation of the expression $$ pN (p+q)^{N-1} $$ does not depend on whether two points are constrained or not. 
This is why the procedure works. 
In brief, one is given a function to evaluate (i.e. a pointwise-operation, completely independent of any constraints among variables), where two variables are related by a constraint. One can view this as an unconstrained function, defined on the whole of $\mathcal{R}^2$.
One assumes a variable is a constant, and performs legitimate operations under such assumption.
One then derives a final result which does not depend on the assumption taken: one safely concludes the final result applies to the case where the two variables are constrained. 
Imagine the contrary were true, i.e. the final result did not apply to a certain $p = \hat{p}$ and $\hat{q} = 1-\hat{p}$.  That would be absurd, as I could have chosen $ \hat{q}$ in the first place as constant.
