binomial()$linkinv(1L) now works

git-svn-id: https://svn.r-project.org/R/trunk@87546 00db46b3-68df-0310-9c12-caf00c1e9a41
r-devel · Jan 9, 2025 · aec3e51 · aec3e51
1 parent 13ede87
commit aec3e51
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 4 deletions.
diff --git a/doc/NEWS.Rd b/doc/NEWS.Rd
@@ -534,6 +534,9 @@
       \item \code{arima(.., seasonal = <wrong-vector>)} correctly errors
       now, ditto for \code{arima0()}, thanks to \I{Norbert Kuder}'s report
       on the R-devel list.
+
+      \item \code{binomial(<link>)$linkinv(eta)} now also works for
+      \code{"logit"} link when \code{is.integer(eta)}.
     }
   }
 }

diff --git a/src/library/stats/src/family.c b/src/library/stats/src/family.c
@@ -65,16 +65,17 @@ SEXP logit_link(SEXP mu)
     double *rans = REAL(ans), *rmu=REAL(mu);
 
     for (i = 0; i < n; i++)
-	rans[i] = log(x_d_omx(rmu[i]));
+	rans[i] = log(x_d_omx(rmu[i]));// log( x/(1-x) )
     UNPROTECT(1);
     return ans;
 }
 
 SEXP logit_linkinv(SEXP eta)
 {
-    int i, n = LENGTH(eta);
-    if (!n || !isReal(eta))
+    int i, n = LENGTH(eta), nprot = 1;
+    if (!n || !isNumeric(eta))
 	error(_("Argument %s must be a nonempty numeric vector"), "eta");
+    if (!isReal(eta)) {eta = PROTECT(coerceVector(eta, REALSXP)); nprot++;}
     SEXP ans = PROTECT(shallow_duplicate(eta));
     double *rans = REAL(ans), *reta = REAL(eta);
 
@@ -84,7 +85,7 @@ SEXP logit_linkinv(SEXP eta)
 	    ((etai > THRESH) ? INVEPS : exp(etai));
 	rans[i] = x_d_opx(tmp);
     }
-    UNPROTECT(1);
+    UNPROTECT(nprot);
     return ans;
 }
 

diff --git a/tests/reg-tests-1e.R b/tests/reg-tests-1e.R
@@ -1722,6 +1722,14 @@ stopifnot(exprs = {
 ## gave solve.default() error (as wrong model failed fitting)
 
 
+##  binomial()$linkinv(<int>)
+lnks <- c("logit", "probit", "cloglog", "cauchit", "log")
+binIlink <- function(eta) sapply(lnks, function(lnk) binomial(lnk)$linkinv(eta))
+stopifnot(identical(binIlink(          0:3),
+                    binIlink(as.double(0:3))))
+## integer type was not allowed for logit (only) in R <= 4.4.2
+
+
 
 ## keep at end
 rbind(last =  proc.time() - .pt,