{"word_count": 435, "character_count": 4968}
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 [latex]Y[/latex] with mean [latex]\mu[/latex] and variance [latex]\sigma^2[/latex]. Using Taylor's expansion,[latex display="true"]g(Y)\sim g(\mu)+(Y-\mu)g'(\mu)+\frac{1}{2}(Y-\mu)^2g''(\mu)[/latex]we can write[latex display="true"]\mathbb{E}[g(Y)]\sim g(\mu)+\frac{\sigma^2}{2}g''(\mu)[/latex] [latex display="true"]\text{Var}[g(Y)]\sim [g'(\mu)]^2\sigma^2[/latex]Assume that [latex]Y\sim\mathcal{P}(\lambda)[/latex], a consider a square root transformation, [latex]g(y)=\sqrt{y}[/latex], then the second equality becomes [latex display="true"]\text{Var}[\sqrt{Y}]\sim \left[\frac{1}{2\sqrt{\mathbb{E}[Y]}}\right]^2\text{Var}[Y]=\frac{1}{4}[/latex]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 modelAssume that the true variable of interest is a Poisson one, [latex]N|\mathbf{X}=\mathbf{x}\sim\mathcal{P}(\lambda_{\mathbf{x}})[/latex] where [latex]\lambda_{\mathbf{x}}=\exp[\mathbf{x}^T\mathbf{\beta}][/latex]Thus,[latex display="true"]\mathbb{P}[N=0|\mathbf{X}=\mathbf{x}]=\exp[-\lambda_{\mathbf{x}}]=\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])][/latex]while[latex display="true"]\mathbb{P}[N&#62;0|\mathbf{X}=\mathbf{x}]=1-\exp[-(\exp[\mathbf{x}^T\mathbf{\beta}])]=H(\mathbf{x}^T\mathbf{\beta})[/latex]where [latex]H(s)=1-\exp[-\exp(s)][/latex]. Let [latex]Y=\mathbf{1}(N&#62;0)[/latex]. The previous model seems like a Bernoulli regression with [latex]H[/latex] as link function,[latex display="true"]\mathbb{P}[Y=1|\mathbf{X}=\mathbf{x}]=H(\mathbf{x}^T\mathbf{\beta})[/latex]So, assume now that instead of observing [latex]N[/latex] we observe [latex]Y=\boldsymbol{1}(N&#62;0)[/latex]. 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 [latex]e^{\lambda_{\mathbf{x}}}[/latex] and [latex]p_{\mathbf{x}}[/latex] obtained from a standard logistic regressionn=563set.seed(1)base=data.frame(X1=rnorm(n),X2=rnorm(n))lambda=base$X1+base$X2base$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 [latex]p_{\mathbf{x}}[/latex] was obtained from a Bernoulli regression, with a cloglog link function ?regBinom = glm((Y&#62;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$SEXbase$SEX="M"base$SEX[x=="0"]="F"x=base$CHILDRENbase$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 oneregBinom = glm((Y&#62;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(&#62;|z|)    (Intercept) -0.002274   0.048413  -0.047    0.963    X1           1.019814   0.026186  38.945   &#60;2e-16 ***X2           1.004814   0.024172  41.570   &#60;2e-16 *** Zero-inflation model coefficients (binomial with logit link):             Estimate Std. Error z value Pr(&#62;|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 ' ' 1Hence, 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...