1
$\begingroup$

My question regards the use of the gamlss package. I am using gamlss package to fit a dataset to a logistic function. There is only one predictor, let me denote it with x and because the overall dependence is not exactly sigmoid, a better model is achieved by wrapping the predictor x in a smoother function. I chose the cubic spline smoother (cs). The response is binomial. I'll denote it with y: 

y = cbind(number of successful events, total number - number of successful events).

The R code is the following:

cs15 <- gamlss(y~cs(x, df=15), sigma.fo=~1, family=BI,  data=mydata, trace=FALSE)

I want to estimate predicted response for a set of predictors not contained in the original data set. I know this can be achieved with the function: 

cs15fit = predict(cs25, newdata=data.frame(x=xnew), type="response")

However, my problem is that I also want to estimate the standard errors, which should be done by adding se.fit=T:

cs15fit=predict(cs25, newdata=data.frame(x=xnew), type="response", se.fit=T)

But the addition of se.fit = T produces the following error:

se.fit = TRUE is not supported for new data values at the moment

Anyone know how I can still find standard errors for the new values?

Improve this question
$\endgroup$
2
  • $\begingroup$ Is your question about how to get the code to work, or about the issues underlying getting standard errors here? $\endgroup$
    – gung - Reinstate Monica
    Mar 10, 2015 at 21:09
  • $\begingroup$ The code works. The issue is that obtaining standard errors through this package (gamlss) for this case is not possible . So if there is any work around, i.e. how can I find them myself. $\endgroup$
    – user70844
    Mar 10, 2015 at 21:39

2 Answers 2

Reset to default
2
$\begingroup$

It sounds like you want prediction intervals. The se.fit argument name reflects that it will give you standard errors extracted from the original fit, and it sounds like it does this only if you do not specify newdata -- i.e. want the prediction of the original fit. (And my suspicion is that it might be giving you standard errors of your coefficients rather than prediction intervals.)

So you want to search for 'prediction interval of GAM' rather than 'standard errors'. I found a couple of papers that way, but that's just the GAM case and GAMLSS (being more flexible) is going to be even more difficult to get right.

I'd suggest falling back to a GLM -- are you sure that your "better model" isn't overfitting and that it makes a substantial difference if it's not -- where things are more straightforward and the posting @boro141 links to is applicable. (Note that you have to do more than get the se.fit results even in that case.)

http://opensiuc.lib.siu.edu/cgi/viewcontent.cgi?article=1095&context=gs_rp

Improve this answer
$\endgroup$
1
$\begingroup$

See the answer to this question: How are the standard errors computed for the fitted values from a logistic regression?

The linked answer relates to a glm. I am not sure whether using a smoothing term as you have done affects the standard error calculation, but it should give you an idea of how standard errors are computed in predict.

Improve this answer
$\endgroup$
1
  • $\begingroup$ +1. They key here is that these "standard errors" are prediction intervals, not standard errors of coefficients. The description by the OP and the title of the linked posting are both misleading, but the linked posting does give the correct answer. $\endgroup$
    – Wayne
    Feb 8, 2016 at 16:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.