Sine link with binary regression

Question

I have used the SIN link to estimate probabilities, mostly with Program MARK. However, I am not sure how the SIN link works. I know the SIN link enables parameter estimation on boundaries (when the probability is 0 or 1), while the logit link seems to have problems in those situations.

Given that the SIN link works in Program MARK, it seems to me that it should also work with logistic regression in glm in Program R since both are estimating proportions of successes out of n trials. However, there does not appear to be a built-in SIN link in R for glm. I asked a question about this on Stack Overflow when I thought it would be a simple programming issue and learned R does allow user-specified links in glm.

However, nobody posted the R code for the SIN link and I have encountered problems with the inverse of the link when trying to do it myself. Namely, there may be multiple solutions given the periodic nature of the SIN function. I imagine the search space must be constrained somehow to avoid multiple solutions.

The issue now seems more like a statistical problem than a programming problem and I am posting here hoping someone will either post the R code for the SIN link in glm or provide a description of how the SIN link works sufficient for me to codify it myself. The other option I considered was to post a bounty on my Stack Overflow question and I might still do that later.

Specifically, I am stuck on how to create the inverse of the link and whether the search space must be constrained to avoid multiple solutions resulting from the periodic nature of the SIN function. I am aware of the family statement in glm and the make.link statement and that user-specified links are possible in R, for example with the logexp link (?family in R).

Here is a link to my earlier Stack Overflow question:

https://stackoverflow.com/questions/25431923/using-the-sin-link-with-logistic-regression

The SIN link is:

(sin(X * Beta) + 1) / 2

where:

X = design matrix
Beta = parameter estimates

From my earlier post, note that sin(z) is bounded between -1 and +1 inclusive. Therefore, (sin(z) + 1) is bounded between 0 and 2 inclusive. Therefore, (sin(z) + 1)/2 is bounded between 0 and 1 inclusive. This inclusivity is probably what enables the link to return estimates of proportions even when probability is 0 or 1.

My Stack Overflow post also contains R code for an 'exhaustive search' further suggesting the link probably does work for estimation of proportions if the search space is somehow constrained.

Thank you for any help. Hopefully this is an acceptable use of cross-posting given the statistical focus of the present question and the programming focus of the earlier question.

I should point out that in Program MARK the SIN link seems to be reserved for models that contain only one entry per row in the design matrix, i.e., without covariates.

EDIT

The inverse of the SIN link, as near as I can tell, is: lin.pred <- asin((2*y) - 1) This returns the correct linear predictor as long as the linear predictor, z, falls in the range -1 <= z <= 1. If covariates are not allowed, as alluded to above, then perhaps this inverse function will always be acceptable. Of course, I am still going to assume for now that dummy variables are allowed in the linear predictor (at least if there is no intercept).

The above proposed inverse link might solve the immediate problem of the inverse and constraints on the search space. The next step is to incorporate this inverse link into a user-specified link function in glm in R along with any additional required information to estimate proportions.

Answer 1

I'm not clear on why you are starting a new topic on the subject, nor am I clear on the motivation for pursuing a $\sin$ link. A binary model with a $\sin$ link is no longer a logistic model. And there are no boundary problems with ordinary logistic regression. You seem to find something wrong with infinite parameter estimates in order to achieve probabilities of zero or one but in fact this does not present a problem for the logistic model.

If you wanted to fit a variance-stabilizing model you would use the $X\beta = \arcsin(\sqrt p)$ model, yielding $p = \sin^{2}(X\beta)$. This model looks very appealing (information matrix = $X'X$ times a constant) but when you plot the log-likelihood surface you find it is too flat to lead to efficient estimates of $\beta$.

Answer 2

Instead of trying to implement the SIN link into glm in R I decided to try to write the likelihood function for the logit link and solve with optim. Then I repeated that approach simply substituting in the SIN link. Below is the R code for both. Both give the correct point estimate to three decimal places. The logit link did a better job. I was expecting better performance from the SIN link than what I observed. Perhaps I made an error somewhere.

The code for the logit link is modified from that found here:

Calculate coefficients in a logistic regression with R

I will update my answer if I learn more.

my.data <- read.table(text='
   y  x
   0  1
   0  1
   1  1
   1  1
   1  1
   1  1
   1  1
   1  1
   1  1
   1  1
', header = TRUE)

# create the design matrices
vY = as.matrix(my.data$y)
mX = as.matrix(my.data$x)

logLikelihoodLogit = function(vBeta, mX, vY) {
  return( -sum( vY * log( (1/(1+exp(-(mX %*% vBeta)))) ) + (1-vY) * log(1 - (1/(1+exp(-(mX %*% vBeta)))) )))
}

logLikelihoodSin = function(vBeta, mX, vY) {
  return( -sum( vY * log( ((sin(mX * vBeta) + 1) / 2)  ) + (1-vY) * log(1 - ((sin(mX * vBeta) + 1) / 2 ) ))) 
}

# set initial parameter values
vBeta0 = c(0.5) # arbitrary starting parameters

# minimize the negative log-likelihood for logit link
optimLogit = optim(vBeta0, logLikelihoodLogit, mX = mX, vY = vY, method = 'BFGS', hessian=TRUE)
optimLogit$par
# [1] 1.386305

(1/(1+exp(-(optimLogit$par))))
# [1] 0.8000017

# minimize the negative log-likelihood for SIN link
optimSin   = optim(vBeta0, logLikelihoodSin  , mX = mX, vY = vY, method = 'BFGS', hessian=TRUE)
optimSin$par
# [1] 0.6439062

(sin(optimSin$par) + 1) / 2
# [1] 0.800162

EDIT - August 26, 2014

I forgot that when only estimating one parameter it is better to use method='Brent'.

When I do that I obtain better point estimates. Both links now return a point estimate matching the expected value:

It might be okay to use lower = 0, upper = 1 with optimSin.b.

optimLogit.b = optim(vBeta0, logLikelihoodLogit, mX = mX, vY = vY, method='Brent', lower = -20, upper = 20)
optimLogit.b

(1/(1+exp(-(optimLogit.b$par))))
# [1] 0.8

optimSin.b = optim(vBeta0, logLikelihoodSin, mX = mX, vY = vY, method='Brent', lower = -20, upper = 20)
optimSin.b

(sin(optimSin.b$par) + 1) / 2
# [1] 0.8