Hi, I am not familiar with gamlss package nor beta inflated distribution but I think it should do the job. From your reference, p0 = nu /(1+nu+tau) & p1 = tau/(1+nu+tau).

In raw in the piece of code below, you will find the lines which should solve your problem.

Best regards,

Geoffrey

## Code for Janka V. ##

rm(list = ls())

require(gamlss)

set.seed(28122016)

rho = 0.2

N.obs = 2000

X1 = rnorm(N.obs,0,1)

X2 = rho *X1+sqrt(1-rho^2)*rnorm(N.obs,0,1)

Y = matrix(X1^2-X2^2, ncol = 1)

Y = apply(Y, MARGIN = 1, FUN = min, 1 )

Y = apply(matrix(Y, ncol = 1), MARGIN = 1, FUN = max, 0)

tab = as.data.frame(cbind(X1,X2,Y))

gam.regression

Hi, I am not familiar with gamlss package nor beta inflated distribution but I think it should do the job. From your reference, p0 = nu /(1+nu+tau) & p1 = tau/(1+nu+tau).

In raw in the piece of code below, you will find the lines which should solve your problem.

Best regards,

Geoffrey

## Code for Janka V. ##

rm(list = ls())

require(gamlss)

set.seed(28122016)

rho = 0.2

N.obs = 2000

X1 = rnorm(N.obs,0,1)

X2 = rho *X1+sqrt(1-rho^2)*rnorm(N.obs,0,1)

Y = matrix(X1^2-X2^2, ncol = 1)

Y = apply(Y, MARGIN = 1, FUN = min, 1 )

Y = apply(matrix(Y, ncol = 1), MARGIN = 1, FUN = max, 0)

tab = as.data.frame(cbind(X1,X2,Y))

gam.regression

Geoffry is right, I think.

I simplified his example to an "intercept-only model" to mak it a bit easier to understand:

require(gamlss)

set.seed(28122016)

N.beta

 Hi Geoffrey and Jochen,

Thank you very much for your clear answers. I have two small follow up questions if that is fine?

1. With regard to Geoffrey's answer: what exactly does the delta.0 and delta.1 mean that you calculated? The closer to zero, the more accurate the models estimates are?

2. With regard to Jochens answer: you say I have to compare p.beta with mean(Y). I didn't know this. I thought I had to compare it with the mean(Y) of the range (0,1) (so the dataset where 0 and 1 are excluded since they are covered by nu and tau). Is that correct? If you are not familiar with the problem and you just jumped in for the coding issues: no problem and thank you very much!!

Thank you very much!

Yes you are correct, that was a typo. p.beta estimates mean (Y.beta). Sorry.

No problem! Thanks a lot!

Janka,

By construction of the model, these two deltas should be very close to zero.

This means that : On the population which has been used to calibrate the model, the mean predicted probability of being 0 is close to the observed probability of being 0 (in this specific population). Same for "probability of being 1".

If you are more familiar with linear regression, you can do the analogy with the role of the intercept in the linear regression model. Inlcuding an intercept in the linear regression model allows you to match the mean of Y (on the observations used to calibrate the model) with the mean of predicted Y (predictions of the model on the observations used to calibrate the model).

Please, find a tool example in the code below.

Regards,

Geoffrey

############ Code for Janka V.

rm(list = ls())

set.seed(30122016)

N.obs

Perfect! Thank you very much for your help!

@ Janka Vanschoenwinkel

Perhaps my remark is too far from the problem, in particular no code is suggested. But I hope it can help in other similar problems with technicalities. Simply, I would like to remind the well know fact, that Maximum Likelihood Method (MLM) estimators  of parameters p_0 and p_1 of the inflated beta pd   are equal n_0/n and n_1/n, respectively, where n_0 and n_1 are the sizes of the separated parts of zero- and one-valued data, respectively. However this holds only in the cases, when the parameters of the core-beta pd (or any other pd concetrated outside the two poits 0 and 1) are independent of p_0 and p_1 . In particular, this will not work in any case when the estimated inflated beta pd belongs a'priori to the family with, say \alpha= A \cdot p_0 and \beta = B \cdot p_0 (in terms of formula 10.47).

Other problem is if the built-in algorithms for BEINF of gamlss are based not on the MLM, e.g. completely based on some moments method. Then the estimators of p_0 and p_1 may differ from the above implied by MLM.

The aim of my remark is to suggest using the prediction of the parameters in two steps:

1-st: to calculate estimators of p_0 and p_1 as above mentioned and to separate the data contained in the open interval (0;1) as a new_data of size n - n_1 - n_0.

2-nd: to estimate parameters \alpha and \beta from the new_data with the use of BE.

The other models for BEINF0 or BEINF1 can be estimated similarly after assumption that p_1=0 or p_0=0, respectively. The new_data should then contain the 1-valued or 0-valued original data.

Regards

An additional remark: 

Obviously, a direct application of BEINF is irreplaceable when the observations corresponding to  1  or  0  are only near these values (i.e. when the fraction of land fertilized  equals almost  100% or  almost 0%, respectively), unless rounding of  the rough data to number  1  or  0 is possible and acceptable.

Regards