使用 glm 拟合逻辑回归的默认起始值
Default starting values fitting logistic regression with glm
我想知道 glm 中的默认起始值是如何指定的。
这个post suggests that default values are set as zeros. This one说它背后有一个算法，但是相关的link被破坏了。
我尝试用算法跟踪拟合简单逻辑回归模型：
set.seed(123)
x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)
# to see parameter estimates in each step
trace(glm.fit, quote(print(coefold)), at = list(c(22, 4, 8, 4, 19, 3)))
首先，没有指定初始值：
glm(y ~ x, family = "binomial")
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
NULL
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995188 1.1669508
第一步，初始值为NULL。
其次，我将起始值设置为零：
glm(y ~ x, family = "binomial", start = c(0, 0))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0 0
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3177530 0.9097521
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3909975 1.1397163
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3994147 1.1666173
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995191 1.1669518
我们可以看到第一种方法和第二种方法之间的迭代不同。
要查看 glm 指定的初始值，我尝试仅用一次迭代来拟合模型：
glm(y ~ x, family = "binomial", control = list(maxit = 1))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
NULL
Call: glm(formula = y ~ x, family = "binomial", control = list(maxit = 1))
Coefficients:
(Intercept) x
0.3864 1.1062
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 134.6
Residual Deviance: 115 AIC: 119
参数估计（毫不奇怪）对应于第二次迭代中第一种方法的估计，即 [1] 0.386379 1.106234
将这些值设置为初始值会导致与第一种方法相同的迭代序列：
glm(y ~ x, family = "binomial", start = c(0.386379, 1.106234))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995188 1.1669508
那么问题来了，这些值是怎么计算出来的？
TL;DR
start=c(b0,b1) 将 eta 初始化为 b0+x*b1 (mu to 1/(1+exp(-eta)))
start=c(0,0) 将 eta 初始化为 0（mu 为 0.5），无论 y 或 x 值如何。
start=NULL 如果 y=1，无论 x 值如何，都会初始化 eta=1.098612 (mu=0.75)。
start=NULL 如果 y=0，无论 x 值如何，都会初始化 eta=-1.098612 (mu=0.25)。
一旦计算出 eta（以及随后的 mu 和 var(mu)），就会计算出 w 和 z 并将其发送到 QR 求解器，本着qr.solve(cbind(1,x) * w, z*w).
长格式
基于 Roland 的评论：我做了一个 glm.fit.truncated()，我把 glm.fit 带到了 C_Cdqrls 调用，然后将其注释掉。 glm.fit.truncated 输出 z 和 w 值（以及用于计算 z 和 w 的数量值），然后将其传递给C_Cdqrls 通话：
## call Fortran code via C wrapper
fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * w, z * w,
min(1e-7, control$epsilon/1000), check=FALSE)
可以阅读更多有关 C_Cdqrls here 的内容。幸运的是，基础 R 中的函数 qr.solve 直接接入了在 glm.fit().
中调用的 LINPACK 版本
所以我们 运行 glm.fit.truncated 为不同的起始值规范，然后用 w 和 z 值调用 qr.solve，我们看到“起始值”如何"（或第一个显示的迭代值）被计算出来。正如 Roland 指出的那样，在 glm() 中指定 start=NULL 或 start=c(0,0) 会影响 w 和 z 的计算，not for start.
对于 start=NULL：z 是一个向量，其中元素的值为 2.431946 或 -2.431946，w 是一个向量，其中所有元素均为 0.4330127：
start.is.null <- glm.fit.truncated(x,y,family=binomial(), start=NULL)
start.is.null
w <- start.is.null$w
z <- start.is.null$z
## if start is NULL, the first displayed values are:
qr.solve(cbind(1,x) * w, z*w)
# > qr.solve(cbind(1,x) * w, z*w)
# x
# 0.386379 1.106234
对于 start=c(0,0)：z 是一个向量，其中元素的值为 2 或 -2，w 是一个向量，其中所有元素均为 0.5：
## if start is c(0,0)
start.is.00 <- glm.fit.truncated(x,y,family=binomial(), start=0)
start.is.00
w <- start.is.00$w
z <- start.is.00$z
## if start is c(0,0), the first displayed values are:
qr.solve(cbind(1,x) * w, z*w)
# > qr.solve(cbind(1,x) * w, z*w)
# x
# 0.3177530 0.9097521
这很好，但是我们如何计算 w 和 z？在 glm.fit.truncated() 的底部附近，我们看到
z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])
查看以下用于计算 z 和 w 的量的输出值之间的比较：
cbind(y, start.is.null$mu, start.is.00$mu)
cbind(y, start.is.null$eta, start.is.00$eta)
cbind(start.is.null$var_mu, start.is.00$var_mu)
cbind(start.is.null$mu.eta.val, start.is.00$mu.eta.val)
请注意，start.is.00 的矢量 mu 只有值 0.5，因为 eta 设置为 0 且 mu(eta) = 1/(1+exp(-0))= 0.5 . start.is.null 将 y=1 的设置为 mu=0.75（对应于 eta=1.098612），将 y=0 的设置为 mu=0.25（对应于 eta=-1.098612），因此 var_mu = 0.75*0.25 = 0.1875.
然而，值得注意的是，我更改了种子并重新运行所有内容，y=1 的 mu=0.75 和 y=0 的 mu=0.25（因此其他数量保持不变）。也就是说，无论 y 和 x 是什么，start=NULL 都会产生相同的 w 和 z，因为它们初始化 eta=1.098612 (mu=0.75 ) 如果 y=1 且 eta=-1.098612 (mu=0.25) 如果 y=0.
因此，截距系数和 X 系数的起始值似乎未设置为 start=NULL，而是根据 y 值并独立于 x- 将初始值赋予 eta价值。从那里计算 w 和 z，然后与 x 一起发送到 qr.solver。
代码到 运行 在 上面的块之前：
set.seed(123)
x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)
glm.fit.truncated <- function(x, y, weights = rep.int(1, nobs),
start = 0,etastart = NULL, mustart = NULL,
offset = rep.int(0, nobs),
family = binomial(),
control = list(),
intercept = TRUE,
singular.ok = TRUE
){
control <- do.call("glm.control", control)
x <- as.matrix(x)
xnames <- dimnames(x)[[2L]]
ynames <- if(is.matrix(y)) rownames(y) else names(y)
conv <- FALSE
nobs <- NROW(y)
nvars <- ncol(x)
EMPTY <- nvars == 0
## define weights and offset if needed
if (is.null(weights))
weights <- rep.int(1, nobs)
if (is.null(offset))
offset <- rep.int(0, nobs)
## get family functions:
variance <- family$variance
linkinv <- family$linkinv
if (!is.function(variance) || !is.function(linkinv) )
stop("'family' argument seems not to be a valid family object", call. = FALSE)
dev.resids <- family$dev.resids
aic <- family$aic
mu.eta <- family$mu.eta
unless.null <- function(x, if.null) if(is.null(x)) if.null else x
valideta <- unless.null(family$valideta, function(eta) TRUE)
validmu <- unless.null(family$validmu, function(mu) TRUE)
if(is.null(mustart)) {
## calculates mustart and may change y and weights and set n (!)
eval(family$initialize)
} else {
mukeep <- mustart
eval(family$initialize)
mustart <- mukeep
}
if(EMPTY) {
eta <- rep.int(0, nobs) + offset
if (!valideta(eta))
stop("invalid linear predictor values in empty model", call. = FALSE)
mu <- linkinv(eta)
## calculate initial deviance and coefficient
if (!validmu(mu))
stop("invalid fitted means in empty model", call. = FALSE)
dev <- sum(dev.resids(y, mu, weights))
w <- sqrt((weights * mu.eta(eta)^2)/variance(mu))
residuals <- (y - mu)/mu.eta(eta)
good <- rep_len(TRUE, length(residuals))
boundary <- conv <- TRUE
coef <- numeric()
iter <- 0L
} else {
coefold <- NULL
eta <-
if(!is.null(etastart)) etastart
else if(!is.null(start))
if (length(start) != nvars)
stop(gettextf("length of 'start' should equal %d and correspond to initial coefs for %s", nvars, paste(deparse(xnames), collapse=", ")),
domain = NA)
else {
coefold <- start
offset + as.vector(if (NCOL(x) == 1L) x * start else x %*% start)
}
else family$linkfun(mustart)
mu <- linkinv(eta)
if (!(validmu(mu) && valideta(eta)))
stop("cannot find valid starting values: please specify some", call. = FALSE)
## calculate initial deviance and coefficient
devold <- sum(dev.resids(y, mu, weights))
boundary <- conv <- FALSE
##------------- THE Iteratively Reweighting L.S. iteration -----------
for (iter in 1L:control$maxit) {
good <- weights > 0
varmu <- variance(mu)[good]
if (anyNA(varmu))
stop("NAs in V(mu)")
if (any(varmu == 0))
stop("0s in V(mu)")
mu.eta.val <- mu.eta(eta)
if (any(is.na(mu.eta.val[good])))
stop("NAs in d(mu)/d(eta)")
## drop observations for which w will be zero
good <- (weights > 0) & (mu.eta.val != 0)
if (all(!good)) {
conv <- FALSE
warning(gettextf("no observations informative at iteration %d",
iter), domain = NA)
break
}
z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])
# ## call Fortran code via C wrapper
# fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * w, z * w,
# min(1e-7, control$epsilon/1000), check=FALSE)
#
#print(iter)
#print(z)
#print(w)
}
}
return(list(z=z, w=w, mustart=mustart, etastart=etastart, eta=eta, offset=offset, mu=mu, mu.eta.val=mu.eta.val,
weight=weights, var_mu=variance(mu)))
}
我想知道 glm 中的默认起始值是如何指定的。
这个post suggests that default values are set as zeros. This one说它背后有一个算法，但是相关的link被破坏了。
我尝试用算法跟踪拟合简单逻辑回归模型：
set.seed(123)
x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)
# to see parameter estimates in each step
trace(glm.fit, quote(print(coefold)), at = list(c(22, 4, 8, 4, 19, 3)))
首先，没有指定初始值：
glm(y ~ x, family = "binomial")
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
NULL
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995188 1.1669508
第一步，初始值为NULL。
其次，我将起始值设置为零：
glm(y ~ x, family = "binomial", start = c(0, 0))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0 0
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3177530 0.9097521
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3909975 1.1397163
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3994147 1.1666173
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995191 1.1669518
我们可以看到第一种方法和第二种方法之间的迭代不同。
要查看 glm 指定的初始值，我尝试仅用一次迭代来拟合模型：
glm(y ~ x, family = "binomial", control = list(maxit = 1))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
NULL
Call: glm(formula = y ~ x, family = "binomial", control = list(maxit = 1))
Coefficients:
(Intercept) x
0.3864 1.1062
Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
Null Deviance: 134.6
Residual Deviance: 115 AIC: 119
参数估计（毫不奇怪）对应于第二次迭代中第一种方法的估计，即 [1] 0.386379 1.106234
将这些值设置为初始值会导致与第一种方法相同的迭代序列：
glm(y ~ x, family = "binomial", start = c(0.386379, 1.106234))
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.386379 1.106234
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3991135 1.1653971
Tracing glm.fit(x = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, .... step 22,4,8,4,19,3
[1] 0.3995188 1.1669508
那么问题来了，这些值是怎么计算出来的？
TL;DR
start=c(b0,b1)将 eta 初始化为b0+x*b1(mu to 1/(1+exp(-eta)))start=c(0,0)将 eta 初始化为 0（mu 为 0.5），无论 y 或 x 值如何。start=NULL如果 y=1，无论 x 值如何，都会初始化 eta=1.098612 (mu=0.75)。start=NULL如果 y=0，无论 x 值如何，都会初始化 eta=-1.098612 (mu=0.25)。一旦计算出 eta（以及随后的 mu 和 var(mu)），就会计算出
w和z并将其发送到 QR 求解器，本着qr.solve(cbind(1,x) * w, z*w).
长格式
基于 Roland 的评论：我做了一个 glm.fit.truncated()，我把 glm.fit 带到了 C_Cdqrls 调用，然后将其注释掉。 glm.fit.truncated 输出 z 和 w 值（以及用于计算 z 和 w 的数量值），然后将其传递给C_Cdqrls 通话：
## call Fortran code via C wrapper
fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * w, z * w,
min(1e-7, control$epsilon/1000), check=FALSE)
可以阅读更多有关 C_Cdqrls here 的内容。幸运的是，基础 R 中的函数 qr.solve 直接接入了在 glm.fit().
所以我们 运行 glm.fit.truncated 为不同的起始值规范，然后用 w 和 z 值调用 qr.solve，我们看到“起始值”如何"（或第一个显示的迭代值）被计算出来。正如 Roland 指出的那样，在 glm() 中指定 start=NULL 或 start=c(0,0) 会影响 w 和 z 的计算，not for start.
对于 start=NULL：z 是一个向量，其中元素的值为 2.431946 或 -2.431946，w 是一个向量，其中所有元素均为 0.4330127：
start.is.null <- glm.fit.truncated(x,y,family=binomial(), start=NULL)
start.is.null
w <- start.is.null$w
z <- start.is.null$z
## if start is NULL, the first displayed values are:
qr.solve(cbind(1,x) * w, z*w)
# > qr.solve(cbind(1,x) * w, z*w)
# x
# 0.386379 1.106234
对于 start=c(0,0)：z 是一个向量，其中元素的值为 2 或 -2，w 是一个向量，其中所有元素均为 0.5：
## if start is c(0,0)
start.is.00 <- glm.fit.truncated(x,y,family=binomial(), start=0)
start.is.00
w <- start.is.00$w
z <- start.is.00$z
## if start is c(0,0), the first displayed values are:
qr.solve(cbind(1,x) * w, z*w)
# > qr.solve(cbind(1,x) * w, z*w)
# x
# 0.3177530 0.9097521
这很好，但是我们如何计算 w 和 z？在 glm.fit.truncated() 的底部附近，我们看到
z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])
查看以下用于计算 z 和 w 的量的输出值之间的比较：
cbind(y, start.is.null$mu, start.is.00$mu)
cbind(y, start.is.null$eta, start.is.00$eta)
cbind(start.is.null$var_mu, start.is.00$var_mu)
cbind(start.is.null$mu.eta.val, start.is.00$mu.eta.val)
请注意，start.is.00 的矢量 mu 只有值 0.5，因为 eta 设置为 0 且 mu(eta) = 1/(1+exp(-0))= 0.5 . start.is.null 将 y=1 的设置为 mu=0.75（对应于 eta=1.098612），将 y=0 的设置为 mu=0.25（对应于 eta=-1.098612），因此 var_mu = 0.75*0.25 = 0.1875.
然而，值得注意的是，我更改了种子并重新运行所有内容，y=1 的 mu=0.75 和 y=0 的 mu=0.25（因此其他数量保持不变）。也就是说，无论 y 和 x 是什么，start=NULL 都会产生相同的 w 和 z，因为它们初始化 eta=1.098612 (mu=0.75 ) 如果 y=1 且 eta=-1.098612 (mu=0.25) 如果 y=0.
因此，截距系数和 X 系数的起始值似乎未设置为 start=NULL，而是根据 y 值并独立于 x- 将初始值赋予 eta价值。从那里计算 w 和 z，然后与 x 一起发送到 qr.solver。
代码到 运行 在 上面的块之前：
set.seed(123)
x <- rnorm(100)
p <- 1/(1 + exp(-x))
y <- rbinom(100, size = 1, prob = p)
glm.fit.truncated <- function(x, y, weights = rep.int(1, nobs),
start = 0,etastart = NULL, mustart = NULL,
offset = rep.int(0, nobs),
family = binomial(),
control = list(),
intercept = TRUE,
singular.ok = TRUE
){
control <- do.call("glm.control", control)
x <- as.matrix(x)
xnames <- dimnames(x)[[2L]]
ynames <- if(is.matrix(y)) rownames(y) else names(y)
conv <- FALSE
nobs <- NROW(y)
nvars <- ncol(x)
EMPTY <- nvars == 0
## define weights and offset if needed
if (is.null(weights))
weights <- rep.int(1, nobs)
if (is.null(offset))
offset <- rep.int(0, nobs)
## get family functions:
variance <- family$variance
linkinv <- family$linkinv
if (!is.function(variance) || !is.function(linkinv) )
stop("'family' argument seems not to be a valid family object", call. = FALSE)
dev.resids <- family$dev.resids
aic <- family$aic
mu.eta <- family$mu.eta
unless.null <- function(x, if.null) if(is.null(x)) if.null else x
valideta <- unless.null(family$valideta, function(eta) TRUE)
validmu <- unless.null(family$validmu, function(mu) TRUE)
if(is.null(mustart)) {
## calculates mustart and may change y and weights and set n (!)
eval(family$initialize)
} else {
mukeep <- mustart
eval(family$initialize)
mustart <- mukeep
}
if(EMPTY) {
eta <- rep.int(0, nobs) + offset
if (!valideta(eta))
stop("invalid linear predictor values in empty model", call. = FALSE)
mu <- linkinv(eta)
## calculate initial deviance and coefficient
if (!validmu(mu))
stop("invalid fitted means in empty model", call. = FALSE)
dev <- sum(dev.resids(y, mu, weights))
w <- sqrt((weights * mu.eta(eta)^2)/variance(mu))
residuals <- (y - mu)/mu.eta(eta)
good <- rep_len(TRUE, length(residuals))
boundary <- conv <- TRUE
coef <- numeric()
iter <- 0L
} else {
coefold <- NULL
eta <-
if(!is.null(etastart)) etastart
else if(!is.null(start))
if (length(start) != nvars)
stop(gettextf("length of 'start' should equal %d and correspond to initial coefs for %s", nvars, paste(deparse(xnames), collapse=", ")),
domain = NA)
else {
coefold <- start
offset + as.vector(if (NCOL(x) == 1L) x * start else x %*% start)
}
else family$linkfun(mustart)
mu <- linkinv(eta)
if (!(validmu(mu) && valideta(eta)))
stop("cannot find valid starting values: please specify some", call. = FALSE)
## calculate initial deviance and coefficient
devold <- sum(dev.resids(y, mu, weights))
boundary <- conv <- FALSE
##------------- THE Iteratively Reweighting L.S. iteration -----------
for (iter in 1L:control$maxit) {
good <- weights > 0
varmu <- variance(mu)[good]
if (anyNA(varmu))
stop("NAs in V(mu)")
if (any(varmu == 0))
stop("0s in V(mu)")
mu.eta.val <- mu.eta(eta)
if (any(is.na(mu.eta.val[good])))
stop("NAs in d(mu)/d(eta)")
## drop observations for which w will be zero
good <- (weights > 0) & (mu.eta.val != 0)
if (all(!good)) {
conv <- FALSE
warning(gettextf("no observations informative at iteration %d",
iter), domain = NA)
break
}
z <- (eta - offset)[good] + (y - mu)[good]/mu.eta.val[good]
w <- sqrt((weights[good] * mu.eta.val[good]^2)/variance(mu)[good])
# ## call Fortran code via C wrapper
# fit <- .Call(C_Cdqrls, x[good, , drop = FALSE] * w, z * w,
# min(1e-7, control$epsilon/1000), check=FALSE)
#
#print(iter)
#print(z)
#print(w)
}
}
return(list(z=z, w=w, mustart=mustart, etastart=etastart, eta=eta, offset=offset, mu=mu, mu.eta.val=mu.eta.val,
weight=weights, var_mu=variance(mu)))
}