Description Usage Arguments Author(s) Examples
Use the adaptive LASSO to fit a GLM in the GWR setting.
1 | gwglmnet.adaptive.fit(x, y, coords, weight.matrix, s, verbose, family, prior.weights)
|
x |
|
y |
|
coords |
|
weight.matrix |
|
s |
|
verbose |
|
family |
|
prior.weights |
Wesley Brooks
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (x, y, coords, weight.matrix, s, verbose, family, prior.weights)
{
gwglmnet.object = list()
coords.unique = unique(coords)
model = list()
s.optimal = vector()
adapt.normx = list()
adapt.scale = list()
cv.error = list()
coef.scale = list()
glm.step = list()
for (i in 1:dim(coords.unique)[1]) {
colocated = which(coords[, 1] == coords.unique[i, 1] &
coords[, 2] == coords.unique[i, 2])
loow = weight.matrix[i, -colocated]
if (sum(loow) == 0) {
return(list(cv.error = Inf))
}
prior.loow = prior.weights[-colocated]
reps = length(colocated)
w <- prior.loow * loow
xx = as.matrix(x[-colocated, ])
yy = as.matrix(y[-colocated])
if (family == "binomial" && (abs(sum(yy * w) - sum(w)) <
1e-04 || sum(yy * w) < 1e-04)) {
cat(paste("Abort. i=", i, ", weighted sum=", sum(yy *
w), ", sum of weights=", sum(w), "\n", sep = ""))
model[[i]] = NULL
cv.error[[i]] = 0
s.optimal = c(s.optimal, max(s))
}
else {
m <- ncol(xx)
n <- nrow(xx)
one <- rep(1, n)
meanx <- drop(one %*% xx)/n
x.centered <- scale(xx, meanx, FALSE)
normx <- sqrt(drop(one %*% (x.centered^2)))
adapt.normx[[i]] = normx
names(normx) <- NULL
xs = x.centered
for (k in 1:dim(x.centered)[2]) {
if (normx[k] != 0) {
xs[, k] = xs[, k]/normx[k]
}
else {
xs[, k] = rep(0, dim(xs)[1])
normx[k] = Inf
}
}
out.glm = try(glm(yy ~ xs, family = family, weights = w))
if (class(out.glm) == "try-error") {
cat(paste("Had to use the last glm for location ",
i, "\n", sep = ""))
glm.step[[i]] = out.glm = glm.step[[i - 1]]
}
else {
glm.step[[i]] = out.glm
}
beta.glm = out.glm$coeff[2:(m + 1)]
adapt.weight = abs(beta.glm)
adapt.scale[[i]] = adapt.weight
for (k in 1:dim(x.centered)[2]) {
if (!is.na(adapt.weight[k])) {
xs[, k] = xs[, k] * adapt.weight[k]
}
else {
xs[, k] = rep(0, dim(xs)[1])
adapt.weight[k] = 0
}
}
coef.scale[[i]] = adapt.weight/normx
names(coef.scale[[i]]) = sapply(strsplit(names(coef.scale[[i]]),
"xs"), function(x) {
x[2]
})
if (sum(coef.scale[[i]]) < 1e-10) {
if (verbose) {
cat(paste("opted for the intercept-only model at location: ",
i, "\n", sep = ""))
}
model[[i]] = NULL
predictions = rep(coef(out.glm)[["(Intercept)"]],
length(colocated))
cv.error[[i]] = abs(matrix(predictions - matrix(y[colocated],
nrow = reps, ncol = length(s))))
s.optimal = c(s.optimal, max(s))
}
else {
if (family == "binomial") {
model[[i]] = glmnet(x = xs, y = cbind(1 - yy,
yy), lambda = s, family = family, weights = w)
}
else {
model[[i]] = glmnet(x = xs, y = yy, lambda = s,
family = family, weights = w)
}
predictions = predict(model[[i]], newx = scale(matrix(x[colocated,
], nrow = reps, ncol = dim(xx)[2]), center = meanx,
scale = normx/adapt.weight), type = "response",
s = s)
cv.error[[i]] = colSums(abs(matrix(predictions -
matrix(y[colocated], nrow = reps, ncol = length(s)),
nrow = reps, ncol = length(s))))
s.optimal = c(s.optimal, s[which.min(cv.error[[i]])])
}
}
if (verbose) {
cat(paste(i, "\n", sep = ""))
}
}
gwglmnet.object[["coef.scale"]] = coef.scale
gwglmnet.object[["model"]] = model
gwglmnet.object[["s"]] = s.optimal
gwglmnet.object[["mode"]] = mode
gwglmnet.object[["coords"]] = coords.unique
gwglmnet.object[["cv.error"]] = cv.error
gwglmnet.object[["s.range"]] = s
class(gwglmnet.object) = "gwglmnet.object"
return(gwglmnet.object)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.