# Betareg error when running bootstrap-like simulations (Error in chol.default(K) : the leading minor of order 16 is not positive definite

A previous post discussed a similar case on non-positive definite covariance matrices resulting when producing half-normal residual plots using the package betareg. However, I would like additional comments on the implications of this error in the context of beta-regression analysis.

I start with the assumption that the data set at hand is a realization of the random process we aim to represent with the beta regression model. If this is the case, would not be expected that, if the model is properly specified, any realization of such a data set (like bootstrap samples) should allow for model fitting without major issues? In other words, if we found that certain random realizations of our data set result in non-positive definite covariance matrices, is this an indication that the model might not be properly specified?

I've been implementing a beta regression analysis to predict the chemical composition of suspended sediments. For data visualization purposes I've decided to resample my data set 5000 times, fit the model, and predict the chemical composition to create high-resolution density plots. I've found that if I specify more complex model structures for the precision submodel (e.g., by including interaction terms), the error in question occurs after a few iterations (5-10). If I exclude some of those interaction terms, I can run the 5000 iterations without any problems. However, the more complex model, when fitted a single time, has lower AIC values and can be significantly different from the simpler model (when compared with lrtest).

Thus, I'm puzzled by these results and wonder if I should stick to the simpler models that do not run into non-positive definite covariance matrices, or should I use the complex models even though they seem to run into non-positive covariance matrices.

The code below provides access to the actual data set and to model fitting:

#Dataset:Available in google drive
install.packages("RCurl")
require(RCurl)
con <- textConnection( url )

#Fitting models
model_0<-betareg(y~treatment*(x1+x2)+season*(x2)+x3+I(x3^2)+x4|treatment*(x1+x2),
data=data,
summary(model_0)

model_1<-betareg(y~treatment*(x1+x2)+season*(x2)+x3+I(x3^2)+x4|treatment*(x1+x2)+season*(x2),
data=data,
summary(model_1)

model_2<-betareg(y~treatment*(x1+x2)+season*(x2)+x3+I(x3^2)+x4|treatment*(x1+x2)+x3+I(x3^2),
data=data,
summary(model_2)

#Half-normal plot of rediduals
plot(model_0, which = 5, type = "deviance", sub.caption = NULL)
plot(model_1, which = 5, type = "deviance", sub.caption = NULL)
plot(model_2, which = 5, type = "deviance", sub.caption = NULL)#It runs into non-positive covariance matrix

#Model selection
AIC(model_0,model_1,model_2)#Model 2 has the lowest AIC
lrtest(model_0,model_2)#It also suggests preference for model_2

#Bootstrap-like predictions
model<-model_2#It runs into non-positive covariance matrix and also warning about optimization failing to converge

list_data_b<-list()

for (i in 1:100){
if (i==101){
break
}
data_b=data[sample(nrow(data),size=nrow(data),replace=TRUE),]
m<-betareg(as.formula(model),data = data_b)
data_b\$y_p<-predict(m,data_b,na.action = na.pass)
list_data_b[[i]]<-data_b
}



About the data: predictor variables were log transformed since processes leading to mobilization of particulate material in streams are typically multiplicative.