recently, I want to know more details about GWmodel package.I typed gwr.binomial.wt
at console in Rstudio:
function (y, x, bw, W.mat, verbose = T)
{
tol = 1e-05
maxiter = 20
dp.n <- nrow(x)
var.n <- ncol(x)
betas <- matrix(nrow = dp.n, ncol = var.n)
betas1 <- betas
S <- matrix(nrow = dp.n, ncol = dp.n)
n = rep(1, length(y))
it.count <- 0
llik <- 0
mu <- 0.5
nu <- 0
if (verbose)
cat(" Iteration Log-Likelihood:(With bandwidth: ",
bw, ")\n=========================\n")
wt2 <- rep(1, dp.n)
repeat {
y.adj <- nu + (y - n * mu)/(n * mu * (1 - mu))
for (i in 1:dp.n) {
W.i <- W.mat[, i]
gwsi <- gw_reg(x, y.adj, W.i * wt2, FALSE, i)
betas1[i, ] <- gwsi[[1]]
}
nu <- gw_fitted(x, betas1)
mu <- exp(nu)/(1 + exp(nu))
old.llik <- llik
llik <- sum(lchoose(n, y) + (n - y) * log(1 - mu/n) +
y * log(mu/n))
if (is.na(llik))
llik <- old.llik
if (verbose)
cat(paste(" ", formatC(it.count, digits = 4, width = 4),
" ", formatC(llik, digits = 4, width = 7),
"\n"))
if (abs((old.llik - llik)/llik) < tol)
break
wt2 <- n * mu * (1 - mu)
it.count <- it.count + 1
if (it.count == maxiter)
break
}
res <- list(wt2, llik, y.adj)
res
}
Then, I don't quite understand the "gw_reg" function. I continued to type gw_reg
in the console of Rstudio:
function (x, y, w, hatmatrix, focus)
{
.Call("GWmodel_gw_reg", PACKAGE = "GWmodel", x, y, w, hatmatrix,
focus)
}
I find that this function doesn't tell me much about how it works. As I understand it, .call()
means that the function is implemented in C/C ++ code. So, how do I look at this function in more detail? Or how do I look at the C code that implements this function?
The C++ code of the function can be found in src/GWmodel.cpp
in the package's source files. The actual C++ function is also called gw_reg
but gets renamed at the point of export by Rcpp.
You can look at the function here on Github to see how it works in context, but the function itself is
List gw_reg(mat x, vec y, vec w, bool hatmatrix, int focus)
{
mat wspan(1, x.n_cols, fill::ones);
mat xtw = trans(x % (w * wspan));
mat xtwx = xtw * x;
mat xtwy = trans(x) * (w % y);
mat xtwx_inv = inv(xtwx);
vec beta = xtwx_inv * xtwy;
if (hatmatrix)
{
mat ci = xtwx_inv * xtw;
mat s_ri = x.row(focus - 1) * ci;
return List::create(
Named("beta") = beta,
Named("S_ri") = s_ri,
Named("Ci") = ci);
}
else
{
return List::create(
Named("beta") = beta);
}
}
I find the easiest way to track down such source code is to find the package on Github and use the search tool to look up the name of the C++ function.