# Output in two-part Fractional Response Model (FRM) package incorrect?

I was experimenting with the Fractional Response Model (FRM) package, and decided to replicate the results using the base GLM package to better understand the theory. I am able to replicate the results, however, I would have expected the frm$yhat2P output to be different. Consider the case where we are estimating a two-part fractional regression model with y = 1 as the relevant boundary value, using a probit binary link function and a logit fractional link function. I have generated some toy data to work with below. We create this function using the frm package as such: set.seed(1337) N <- 250 u <- rnorm(N) X <- cbind(rnorm(N), rnorm(N)) dimnames(X)[[2]] <- c("X1", "X2") ym <- exp(X[, 1]+X[, 2]+u)/(1+exp(X[, 1]+X[, 2]+u)) y <- rbeta(N, ym*20, 20*(1-ym)) y[y > 0.9] <- 1 temp <- frm(y, X, linkbin="probit", linkfrac="logit", type="2P", inf=1) df2 <- as.data.frame(round(temp$yhat2P, 8))
# store results for future reference.


The first-part of the frm model above can be replicated using glm in the following way:

X2 <- as.data.frame(cbind(X,y))
X2$y_binary <- ifelse(y == 1, 1, 0) temp2 <- glm(y_binary ~ X1+X2, family = binomial("probit"), data = X2) # Generates the likelihood of y = 1. summary(temp2) # Produces identical coefficients etc. as the FRM package # first-part.  The second-part of the frm model above can be replicated using glm in the following way: X3 <- X2[X2$y_binary!=1,]
temp3 <- glm(y~X1+X2, family=quasibinomial("logit"), data = X3)
# this generates the fractional response for y[0,1).
summary(temp3)
# Produces identical coefficients etc.
# as the FRM package second-part.


Finally, replicating the frm temp$yhat2P results (df2 produced above), i.e., the overall fitted mean values using glm can be done as such: test <- as.data.frame(round(predict(temp2,newdata=X2, type = "response"),8)) # 1st part. names(test)[1] <- "first_part" test$second_part <- round(predict(temp3,newdata=X2,
type = "response"),8) # 2nd part.
test$final_part <- round(test$first_part*testsecond_part, 8) # final step.  The formula to generate these results is essentially the following (first_part*second_part): $$E(y|x) = Pr[y = 1|x]*E[y | x, y[0,1)]$$ This output, however, is not consistent with how I interpret the frm theory, as found in Ramalho et al's 2011 paper (the author of the package) "Alternative estimating and testing empirical strategies for fractional regression models" and other supporting papers. I would expect the output should align with the following formula instead (for this particular case): $$E(y|x) = E[y | x, y[0,1)]*Pr[y[0,1) | x] + E[y | x, y = 1]*Pr[y = 1|x]$$ Which I believe is produced using glm as in the following way: testwhat_I_expect <- round(test$second_part*(1-test$first_part)
+ (1*test$first_part), 8)  Why is the frm$yhat2P output not consistent with the formula above? Am I misunderstanding the theory behind two-part FRM models? Or is it possible I am misinterpreting the frm package output? Any guidance and/or clarity would be greatly appreciated.