In my previous post on GLMs, I discussed power link functions. But there are much more links that can be used :
- The square root link (for the Poisson model)
Consider some random variable Y with mean \mu and variance \sigma^2. Using Taylor’s expansion,g(Y)\sim g(\mu)+(Y-\mu)g'(\mu)+\frac{1}{2}(Y-\mu)^2g''(\mu)we can write\mathbb{E}[g(Y)]\sim g(\mu)+\frac{\sigma^2}{2}g''(\mu) \text{Var}[g(Y)]\sim [g'(\mu)]^2\sigma^2
Assume that Y\sim\mathcal{P}(\lambda), a consider a square root transformation, g(y)=\sqrt{y}, then the second equality becomes \text{Var}[\sqrt{Y}]\sim \left[\frac{1}{2\sqrt{\mathbb{E}[Y]}}\right]^2\text{Var}[Y]=\frac{1}{4}
So, somehow, with a square-root transformation, we have variance stability, which might be interpreted as some homoscedasticity.
- The complementary log-log function for the Bernoulli model
Assume that the true variable of interest is a Poisson one, N|\mathbf{X}=\mathbf{x}\sim\mathcal{P}(\lambda_{\mathbf{x}}) where \lambda_{\mathbf{x}}=\exp[\mathbf{x}^T\mathbf{\beta}]Thus,\mathbb{P}[N=0|\mathbf{X}=\mathbf{x}]=\exp[-\lambda_{\mathbf{x}}]=\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]while\mathbb{P}[N>0|\mathbf{X}=\mathbf{x}]=1-\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]=H(\mathbf{x}^T\mathbf{\beta})where H(s)=1-\exp[-\exp(s)]. Let Y=\mathbf{1}(N>0). The previous model seems like a Bernoulli regression with H as link function,\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=H(\mathbf{x}^T\mathbf{\beta})
So, assume now that instead of observing N we observe Y=\boldsymbol{1}(N>0). In that case, running a Bernoulli regression with a complementary log-log link function would be the same (?) as running first a Poisson regression on the original data, and then use it on our binary variable, zero vs. non-zero. Let us generate some data, and see what’s going on. Let us compare e^{\lambda_{\mathbf{x}}} and p_{\mathbf{x}} obtained from a standard logistic regression
n=563 set.seed(1) base=data.frame(X1=rnorm(n),X2=rnorm(n)) lambda=base$X1+base$X2 base$Y=rpois(n,exp(lambda)) regPois = glm(Y~.,data=base,family=poisson(link="log")) lambda = predict(regPois,type="response") regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit")) prob = predict(regBinom, type="response") plot(prob,exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") |
What if p_{\mathbf{x}} was obtained from a Bernoulli regression, with a cloglog link function ?
regBinom = glm((Y>0)~.,data=base,family=binomial(link="cloglog")) prob = predict(regBinom, type="response") plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") |
It looks like the fit is very good here ! Now, what if we have real data, like the dataset from A Theory of Extramarital Affairs, by Ray Fair, published in 1978 in the Journal of Political Economy (with 563 observations, and nine variables)
base = read.table("http://freakonometrics.free.fr/baseaffairs.txt",header=TRUE) str(base) x=base$SEX base$SEX="M" base$SEX[x=="0"]="F" x=base$CHILDREN base$CHILDREN="YES" base$CHILDREN[x==0]="NO" regPois = glm(Y~.,data=base,family=poisson(link="log")) lambda = predict(regPois,type="response") regBinom = glm((Y==0)~.,data=base,family=binomial(link="probit")) prob = predict(regBinom, type="response") plot(prob,exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") |
In that case the two models are very different. But actually, so is the second one
regBinom = glm((Y>0)~.,data=base,family=binomial(link="cloglog")) prob = predict(regBinom, type="response") plot(prob,1-exp(-lambda),xlim=0:1,ylim=0:1) abline(a=0,b=1,lty=2,col="red") |
How can we interpret that ? Could it be because the Poisson model is not good ? Actually, if we run a zero-inflated model here,
library(pscl) regZIP = zeroinfl(Y ~ . | ., data = base) summary(regZIP) Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) -0.002274 0.048413 -0.047 0.963 X1 1.019814 0.026186 38.945 <2e-16 *** X2 1.004814 0.024172 41.570 <2e-16 *** Zero-inflation model coefficients (binomial with logit link): Estimate Std. Error z value Pr(>|z|) (Intercept) -4.90190 2.07846 -2.358 0.0184 * X1 -2.00227 0.86897 -2.304 0.0212 * X2 -0.01545 0.96121 -0.016 0.9872 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 |
Hence, we reject here the Poisson distribution assumption, because of the inflation of zeros… It looks like the cloglog link can be used to check if the Poisson distribution is a good model, or not…
(either the company defaults, or not), so that
is said to be exchangeable if for all 
for any permutation 
of 
. A standard example is the case where
, with
and
Since 
, a necessary condition is that
i.e.
.
. 
on 
such that
where 
. Note that 
A nice proof of that result can be found in 
of 
,
that can be inverted in
The idea is then to extend the size of the vector 
, i.e. for all 
, define
so that, if we condition on 
,
but since given the sum of components of 
elements are equally likely, we can write
The first idea is to work on the blue term, and to invocate a theorem of approximation of the hypergeometric distribution 
to a binomial distribution 
, when 
becomes large. Then
Let 
and let 
denote the cumulative distribution function of 
.
The idea is then to write the sum as an integral, with respect to that distribution,
The theorem is then obtained since 
, i.e.
In the case of non-binary sequences, there is an extension of the previous result,
. 
on 
such that
where 
and
For instance, in the Gaussian case mentioned earlier, if
then
where
i.e. conditionally on 
, the 
. The proof can be found in 
where 
are i.i.d. random variables. To go further on exchangeability and related topics, see 
Since
if we assume that – given the latent factor – 
(either the company defaults, or not),
i.e.
Thus, we can derive the (unconditional) distribution of the sum
i.e.
. In the upper corner is plotted the density of the Beta distribution.
of 
random variables is exchangeable if and only if 
‘s are conditionnally independent, conditionnally on some random variable 
.
i.i.d. with distribution 
. Here we note
is a random variable. The idea is to assume that 
. Thus, we need to compute the distribution of 
which is here extremely simple (due to properties of the Gaussian distribution), i.e.
as an estimator of given our sample data (and thus, we also have a confidence interval since we know the distribution of 
. Note, first, that \theta has support 
. So we need a distribution on that support. Why not a beta distribution ? E.g.
