2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You'll need to write a family and link-glm object. Check out ?family and ?make.link and the corresponding code. $\endgroup$ – Momo Aug 24 '14 at 13:30
  • $\begingroup$ @Momo Thank you. I got that far before posting today. I am stuck on how to write the inverse link and whether the search space needs to be constrained. $\endgroup$ – Mark Miller Aug 24 '14 at 13:33
  • $\begingroup$ This is good to know for a good answer. Perhaps edit the question to reflect that you need those elements. The constrained search space bit confuses me too; shouldnt the model allow all real values? $\endgroup$ – Momo Aug 24 '14 at 13:50
  • $\begingroup$ @Momo Thank you. I have added the fifth paragraph to emphasize those points. $\endgroup$ – Mark Miller Aug 24 '14 at 14:00
1
$\begingroup$

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$.

$\endgroup$
  • 1
    $\begingroup$ Do you not have any interest in the questions of why you would want to use it and why the ordinary logistic model already works? $\endgroup$ – Frank Harrell Aug 24 '14 at 13:32
  • 1
    $\begingroup$ Sounds good. Be sure to precisely describe the challenge you face, because the apparent problem may not really be a problem. Note that it is often the case that one chooses a simple interpretable model such as logistic and relaxes assumptions about functional form by expanding continuous predictors using regression splines. $\endgroup$ – Frank Harrell Aug 24 '14 at 15:34
  • 1
    $\begingroup$ I have never before seen it used for this purpose. $\endgroup$ – Frank Harrell Aug 24 '14 at 18:58
  • 1
    $\begingroup$ Mark, have you seen Section 4 of Chapter 6 in the MARK user's manual (available at phidot.org/software/mark/docs/book)? It says, again and again and again, that usually you really don't want to use a sine link. Apparently its sole purpose is to deal with some numerical problems with an "identity" design matrix (which the example in your answer is not) to obtain "initial estimates" (see p. 6-23). In fact, I cannot find any clear explanation in the manual for why one would want to use a sine link in an actual analysis--I find only extensive explanations of the problems it can cause. $\endgroup$ – whuber Aug 26 '14 at 15:11
  • 1
    $\begingroup$ @whuber Thank you. I will look at that section. My understanding was that Program MARK more-or-less using the sine link as the default if there are no covariates, i.e., if there is only one entry per row in the design matrix (regardless of whether that matrix is an identity matrix). I know lots of problems will arise if you use the sine link along with covariates. Maybe I am wrong. $\endgroup$ – Mark Miller Aug 26 '14 at 15:21
1
$\begingroup$

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
$\endgroup$
  • 2
    $\begingroup$ An example of how to implement a complicated link for a GLM in R appears in my answer at stats.stackexchange.com/a/64039. It avoids implementing the likelihood explicitly, leaving as much work as possible to code that's already incorporated in R. $\endgroup$ – whuber Aug 26 '14 at 14:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.