# GAM log link does not work without starting values

I am trying to estimate a GAM regression model using the implementation of gam from the mgcv package. I have a working Gaussian model for the dispersion and a log link for the linear predictors but I receive the error

>"Error in eval(expr, envir, enclos) : cannot find valid starting values: please specify some".


Edit 1 - The exact syntax is

splineWAR <- gam(WAR ~ s(zAge, bs="cr") + s(zAdjProd, bs="cr") + s(zSOPct, bs="cr") + s(zBBPct, bs="cr"), family=gaussian(link="log"), data = mydata,  start=c(0, 0, 0, 0, 0))


I have read the relevant threads here and here but have unable to apply the steps suggested to a multiple regression. For instance, when I try and set start values for the 5 variables in my regression (1 dependent and 4 independent) by adding the start=c(n1, n2, n3, n4, n5) argument (where the n's are the mean of the relevant variable), I receive the same error even though I am seemingly copying the syntax exactly from the first link. Can anyone make a suggestion as to what I should try next? Thanks.

Edit 2 The code in the gam.fit function that runs right before the error is -

if (!(validmu(mu) && valideta(eta)))

stop("Can't find valid starting values: please specify some")

• The five variables to supply starting values for do not include the dependent variable! (It's not multiplied by any parameter in the model, of course.) They should be the constant together with one for each of the four independent variables. Perhaps that is the cause of your problem? If you can, first run an ordinary regression of log(DV) against the IVs and then use the (five) coefficients from that as your start. What happens?
– whuber
Commented Jun 11, 2013 at 19:43
• As it turns out, I was not using the glm function but the gam function. When I run the glm function using the specification described in my question, the model runs as expected. However, when I switch back the gam function, I get the error I described. Not sure what the issue is. Commented Jun 11, 2013 at 19:58
• Given that a GAM of this sort will have far more than 5 coefficients (there will be a column in the model matrix for each basis function that comprises each cubic regression spline), you aren't supplying enough starting values. How about you git the model without the family, save that fit as say foo and then refit your model with start = coef(foo). The last bit might need start = log(coef(foo))...? Commented Jun 11, 2013 at 20:34
• @whuber There's an intercept in the model, so that will be why it is working with 5 starting values in the GLM case. Commented Jun 11, 2013 at 20:35
• I've added an Answer now which I think explains things. This looks like a bug to me. I can email Simon Wood, author of mgcv and clarify for you if you'd like? Commented Jun 11, 2013 at 21:16

Update Simon Wood has fixed this bug in mgcv in version 1.7-25. The entry reads:

* bugs fixed whereby etastart etc were not passed to initial.spg and
get.null.coefs.


There are a number of places where an error message of the kind you show might be produced. The first is in the initialisation of the GAM:

> gaussian()$initialize expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$link == "inverse" && any(y == 0)) || (family$link == "log" && any(y <= 0)))) stop("cannot find valid starting values: please specify some") mustart <- y })  Note the last clause there: (family$link == "log" && any(y <= 0))

The first part of the clause is TRUE in your case, what about the second part? That will fail, but the question then is, why did it fail as is.null(start) should have been FALSE in your case. This code actually gets called in gam.fit through mgcv:::estimate.gam and thence mgcv:::initial.spg as shown by the traceback():

> traceback()
6: stop("cannot find valid starting values: please specify some")
5: eval(expr, envir, enclos)
4: eval(family$initialize) 3: initial.spg(G$X, G$y, G$w, G$family, G$S, G$off, G$L, G$lsp0) 2: estimate.gam(G, method, optimizer, control, in.out, scale, gamma, ...) 1: gam(y ~ s(x0) + s(x1) + s(x2) + s(x3), data = dat, family = gaussian(link = "log"), start = c(1, 2, 3, 4, 5))  If we look in mgcv:::initial.spg we note these lines (ignore the nobs one) start <- etastart <- mustart <- NULL nobs <- nrow(X) eval(family$initialize)


i.e. the code above evaluates the expression I showed earlier. But it does so after scrubbing out start.

I think this is a bug here as I don't see how you can fit the model you want to with the way this is coded. For example, using

 library(mgcv)
set.seed(2) ## simulate some data...
dat <- gamSim(1,n=400,dist="normal",scale=2)
b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat, start = 1:5)


we note that that fails but it tells us how many starting values to provide.

> b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3), data=dat, start = 1:5)
Error in gam.fit(G, family = G$family, control = control, gamma = gamma, : Length of start should equal 37 and correspond to initial coefs.  OK. change this to a "log" link and refit with 37 starting values > b <- gam(y~s(x0)+s(x1)+s(x2)+s(x3),data=dat, + family = gaussian(link = "log"), start = runif(37)) Error in eval(expr, envir, enclos) : cannot find valid starting values: please specify some  and it still fails. Once this problem gets fixed by the author, you'll still need to specify the correct starting values for the number of terms in your spline model. • Thank you so much for the thorough answer. I just tried re-scaling my y values from so that all y > 0 and it worked. Given what you wrote, I have no idea why but I'm glad to finally have made some progress. Commented Jun 11, 2013 at 21:23 • In your case, the last clause (family$link == "log" && any(y <= 0)) would evaluate to FALSE now that you rescaled the data and that would mean the entire set of clauses evaluates to FALSE and hence the error passes. This clearly seems like a bug so I'll email Simon and refer to this Q&A. Commented Jun 11, 2013 at 21:25

GLMs require both a family (or what I call a variance family) and a link function to define the Fisher Scoring algorithm that solves for your parameter estimates. With Poisson variance and log link, this is Poisson regression or regression, attained with the argument family=poisson to glm. However, use the following argument family=binomial(link="log") and you get relative risk regression.

Most families in the R GLM function allow you to specify link="log" as an optional argument to the family object (e.g. gaussian, gamma, poisson). Irregular GLMs, especially ones for which the range of the link function is bigger than that of the fitted mean in the variance, have bizarre constrictions imposed on the parameter space which Fisher Scoring cannot accommodate.

Using the traceback() function is always useful with errors like these. You can also find the iteration where the algorithm diverges by specifying glm.control=list(maxit=1) for a 1 step estimator, glm.control=list(maxit=2) for a 2 step, and so on and so forth. Plotting your $\beta^{(i)}$ estimates for the (i)-th iteration will help you see what's happening before the Hessian becomes singular, Fisher Scoring diverges, and R explodes.

Your problem with starting values may be because you're supplying the means of the response variables when the variables are contrasts between the logs of the means. Hence ratios of logs of averages would be a better starting place. Personally, if this were such an issue, I'd fit a regular GLM to make sure the algorithm isn't universally divergent and use the parameter estimates for that model to start another.

For instance: Feeding forward logistic regression estimates (ORs) to obtain risk ratios (RRs)

## retrospective incidence of something nasty
data <- data.frame(cases=rpois(10, 10), controls=rpois(10, 1000), age=factor(seq(10), labels=c('0-10', '10-20', '20-30', '30-40', '40-50', '50-60', '60-70', '70-80', '80-90', '90-100')))

## logistic regression
fit <- glm(cbind(cases, controls) ~ age, data=data, family=binomial)

## relative risk regression
fit2 <- glm(cbind(cases, controls) ~ age, data=data, family=binomial(link='log'), start=coef(fit))

• Thanks for the detailed reply. In the comment to my question above, I mentioned that I made an error in my original question and was actually using the gam function instead of glm function. When I ran the help function on gam (127.0.0.1:29054/library/mgcv/html/gam.html), I noticed that there was no start parameter. This makes the error all the more confusing. Commented Jun 11, 2013 at 20:10
• You can ignore my comment about it not having a start parameter. The gam.fit internal module does accept the start parameter. The synax I am using unsucessfully is splineWAR <- gam(WAR ~ s(zAge, bs="cr") + s(zAdjProd, bs="cr") + s(zSOPct, bs="cr") + s(zBBPct, bs="cr"), family=gaussian(link="log"), data = mydata, start=c(0, 0, 0, 0, 0)) Commented Jun 11, 2013 at 20:14
• This is most likely because you're using splines in your formula, which changes 1 variable into several (3 or more usually) in the model matrix. Alas, if they're adaptive smoothing splines, I'm certain of the necessary degrees of freedom. You can arbitrarily increase the number of zeros in the model. I cannot replicate the s(x, bs) since s only takes arguments df and spar. Setting df would probably let you know how many zeroes are necessary to start the algorithm. Commented Jun 11, 2013 at 20:23
• @AdamO Whay would that matter. mgcv:::gam() is just a penalised GLM. Non-standard links may need good starting values, as you mention. Commented Jun 11, 2013 at 20:26
• @AdamO Given the syntax used, I'm pretty sure this is gam() from package mgcv, not the same from package gam. Commented Jun 11, 2013 at 20:27

The parameter "start" takes values for the parameters not for the variables in the regression. The model only has 4 parameters (one for each dependent variables) so you should try start=c(0,0,0,0).

• I have tried this as well but I get the same error. Commented Jun 11, 2013 at 20:35
• No, this is not correct at all. Ignoring the spline aspect there are 5 terms in the model, the 4 variables named plus the model intercept or constant term. However, as this is a GAM with cubic splines, the actual number of terms in the model is a lot more because each spline is determined by a set of $k$ basis functions. Commented Jun 11, 2013 at 20:40