Skip to content

Commit 918ebaa

Browse files
lukesonnetlukesonnet
authoredFeb 11, 2018
Merge pull request #21 from lukesonnet/lukesonnet/refactor
Restructure CPP files; rewrite inference
2 parents 37f2384 + 08b1073 commit 918ebaa

20 files changed

+809
-732
lines changed
 

Diff for: ‎.Rbuildignore

+1
Original file line numberDiff line numberDiff line change
@@ -5,3 +5,4 @@ src/deprecated
55
^\.travis\.yml$
66
^appveyor\.yml$
77
^tests/testthat/test_krlogit_fns\.R$
8+
^tests/testthat/test_old_krls\.R$

Diff for: ‎.travis.yml

+1
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@ matrix:
2929

3030
r_packages:
3131
- numDeriv
32+
- KRLS
3233

3334
r_github_packages:
3435
- r-lib/covr

Diff for: ‎NAMESPACE

+2-18
Original file line numberDiff line numberDiff line change
@@ -1,40 +1,23 @@
11
# Generated by roxygen2: do not edit by hand
22

33
S3method(plot,krls2)
4+
S3method(predict,krls2)
45
export(Ktrunc)
56
export(bsearch)
67
export(chunk)
7-
export(euc_dist)
8-
export(fdskrls)
98
export(generateK)
109
export(getDhat)
1110
export(inference.krls2)
12-
export(kern_gauss)
13-
export(kern_gauss_1d)
14-
export(krlogit_fn_trunc)
15-
export(krlogit_gr_trunc)
16-
export(krlogit_hess_trunc_inv)
1711
export(krls)
18-
export(krls_gr_trunc)
19-
export(krls_hess_trunc_inv)
2012
export(lambdab.fn)
2113
export(lambdabsearch)
2214
export(lambdaline)
2315
export(lambdasearch)
2416
export(logistic)
2517
export(looloss)
26-
export(mult_diag)
2718
export(newKernel)
28-
export(new_gauss_kern)
29-
export(partial_logit)
30-
export(predict.krls2)
31-
export(pwmfx)
32-
export(pwmfx_novar)
3319
export(solveForDOptim)
34-
export(solve_for_d_ls)
35-
export(solve_for_d_ls_w)
3620
export(summary.krls2)
37-
export(trace_mat)
3821
import(RSpectra)
3922
importFrom(Rcpp,sourceCpp)
4023
importFrom(grDevices,devAskNewPage)
@@ -44,6 +27,7 @@ importFrom(graphics,plot)
4427
importFrom(stats,as.formula)
4528
importFrom(stats,optim)
4629
importFrom(stats,optimize)
30+
importFrom(stats,plogis)
4731
importFrom(stats,predict)
4832
importFrom(stats,pt)
4933
importFrom(stats,quantile)

Diff for: ‎R/RcppExports.R

+12-32
Original file line numberDiff line numberDiff line change
@@ -1,83 +1,63 @@
11
# Generated by using Rcpp::compileAttributes() -> do not edit by hand
22
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
33

4-
#' @export
5-
mult_diag <- function(x, d) {
6-
.Call(`_KRLS2_mult_diag`, x, d)
7-
}
8-
9-
#' @export
10-
trace_mat <- function(x) {
11-
.Call(`_KRLS2_trace_mat`, x)
12-
}
13-
14-
#' @export
154
krls_gr_trunc <- function(U, D, y, w, fitted, dhat, lambda) {
165
.Call(`_KRLS2_krls_gr_trunc`, U, D, y, w, fitted, dhat, lambda)
176
}
187

19-
#' @export
208
krls_hess_trunc_inv <- function(U, D, w, lambda) {
219
.Call(`_KRLS2_krls_hess_trunc_inv`, U, D, w, lambda)
2210
}
2311

24-
#' @export
2512
krlogit_fn_trunc <- function(par, U, D, y, w, lambda) {
2613
.Call(`_KRLS2_krlogit_fn_trunc`, par, U, D, y, w, lambda)
2714
}
2815

29-
#' @export
3016
krlogit_gr_trunc <- function(par, U, D, y, w, lambda) {
3117
.Call(`_KRLS2_krlogit_gr_trunc`, par, U, D, y, w, lambda)
3218
}
3319

34-
#' @export
3520
partial_logit <- function(K, coef, beta0) {
3621
.Call(`_KRLS2_partial_logit`, K, coef, beta0)
3722
}
3823

39-
#' @export
4024
krlogit_hess_trunc_inv <- function(par, U, D, y, w, lambda) {
4125
.Call(`_KRLS2_krlogit_hess_trunc_inv`, par, U, D, y, w, lambda)
4226
}
4327

44-
#' @export
28+
mult_diag <- function(x, d) {
29+
.Call(`_KRLS2_mult_diag`, x, d)
30+
}
31+
32+
trace_mat <- function(x) {
33+
.Call(`_KRLS2_trace_mat`, x)
34+
}
35+
4536
euc_dist <- function(x1, x2) {
4637
.Call(`_KRLS2_euc_dist`, x1, x2)
4738
}
4839

49-
#' @export
5040
kern_gauss_1d <- function(x1, x2, b) {
5141
.Call(`_KRLS2_kern_gauss_1d`, x1, x2, b)
5242
}
5343

54-
#' @export
5544
kern_gauss <- function(x, b) {
5645
.Call(`_KRLS2_kern_gauss`, x, b)
5746
}
5847

59-
#' @export
6048
new_gauss_kern <- function(newx, oldx, b) {
6149
.Call(`_KRLS2_new_gauss_kern`, newx, oldx, b)
6250
}
6351

64-
#' @export
52+
pwmfx <- function(k, x, coefhat, vcovc_mat, p, p2, b) {
53+
.Call(`_KRLS2_pwmfx`, k, x, coefhat, vcovc_mat, p, p2, b)
54+
}
55+
6556
solve_for_d_ls <- function(y, U, D, lambda) {
6657
.Call(`_KRLS2_solve_for_d_ls`, y, U, D, lambda)
6758
}
6859

69-
#' @export
7060
solve_for_d_ls_w <- function(y, U, D, w, lambda) {
7161
.Call(`_KRLS2_solve_for_d_ls_w`, y, U, D, w, lambda)
7262
}
7363

74-
#' @export
75-
pwmfx <- function(k, x, coefhat, vcovc, p, b) {
76-
.Call(`_KRLS2_pwmfx`, k, x, coefhat, vcovc, p, b)
77-
}
78-
79-
#' @export
80-
pwmfx_novar <- function(k, x, coefhat, p, b) {
81-
.Call(`_KRLS2_pwmfx_novar`, k, x, coefhat, p, b)
82-
}
83-

Diff for: ‎R/inference.R

+315-230
Large diffs are not rendered by default.

Diff for: ‎R/kernels.R

+1-1
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ generateK <- function(X,
2929
U <- truncDat$Utrunc
3030
D <- truncDat$eigvals
3131
} else {
32-
eigobj <- eigen(K, symmetric = T)
32+
eigobj <- eigen(K)
3333
eigvaliszero <- eigobj$values == 0
3434
if(any(eigvaliszero)) {
3535

Diff for: ‎R/krls2.R

+5-5
Original file line numberDiff line numberDiff line change
@@ -279,7 +279,7 @@ krls <- function(# Data arguments
279279
## todo: build distance matrix E first so that we can more quickly compute
280280
## K below
281281

282-
if(!is.null(lambda)) {
282+
if (!is.null(lambda)) {
283283
# Allow range in lambda argument
284284
if (length(lambda) > 1) {
285285
lambdarange <- lambda
@@ -293,7 +293,7 @@ krls <- function(# Data arguments
293293
}
294294

295295
## set chunks for any ranges
296-
if(is.null(lambda) | is.null(b)) {
296+
if (is.null(lambda) || is.null(b)) {
297297
chunks <- chunk(sample(n), hyperfolds)
298298

299299
control <- c(control,
@@ -316,14 +316,14 @@ krls <- function(# Data arguments
316316
###----------------------------------------
317317

318318
## If b ix fixed, compute all of the kernels
319-
if(!is.null(b)) {
319+
if (!is.null(b)) {
320320

321321
## Compute kernel matrix and truncated objects
322322
Kdat <- generateK(X=X,
323323
b=b,
324324
control=control)
325325

326-
if(is.null(lambda)) {
326+
if (is.null(lambda)) {
327327
lambda <- lambdasearch(y=y,
328328
X=X,
329329
Kdat=Kdat,
@@ -332,7 +332,7 @@ krls <- function(# Data arguments
332332
}
333333
} else { # if(is.null(b)
334334

335-
if(is.null(lambda)) {
335+
if (is.null(lambda)) {
336336

337337
hyperOut <- lambdabsearch(y=y,
338338
X=X,

Diff for: ‎R/lambdabsearch.R

+3-3
Original file line numberDiff line numberDiff line change
@@ -264,9 +264,9 @@ lambdaline <-
264264
#Lbound <- 0
265265
Lbound = .Machine$double.eps #CJH: to avoid Inf in next statement
266266

267-
#while(sum(D / (D + Lbound)) > q){
268-
# Lbound <- Lbound+.05
269-
#}
267+
while(sum(D / (D + Lbound)) > q){
268+
Lbound <- Lbound+.05
269+
}
270270
} else {
271271
stopifnot(is.vector(Lbound),
272272
length(Lbound)==1,

Diff for: ‎R/predict.krls2.R

+101-74
Original file line numberDiff line numberDiff line change
@@ -1,74 +1,101 @@
1-
## This file contains functions for predicting outcome variables
2-
## Functions:
3-
## predict.krls2
4-
## logistic
5-
6-
7-
## a predict function for class 'krls2'
8-
#' @export
9-
predict.krls2 <- function(object, newdata, se.fit = FALSE, ...) {
10-
if (class(object) != "krls2") {
11-
warning("Object not of class 'krls2'")
12-
UseMethod("predict")
13-
return(invisible(NULL))
14-
}
15-
16-
if(is.null(object$vcov.c) & se.fit)
17-
stop("se.fit requires the vcov.c. Use summary object with predict if you want stand error for the fits")
18-
# todo: we should probably call inference here if they want it?
19-
20-
newdata <- as.matrix(newdata)
21-
if (ncol(object$X) != ncol(newdata)) {
22-
stop("ncol(newdata) differs from ncol(X) from fitted krls object")
23-
}
24-
25-
newdataK <- newKernel(X = object$X, newData = newdata, whichkernel = object$kernel, b = object$b)
26-
27-
if(object$loss == "logistic") {
28-
yfitted <- logistic(K = newdataK, coeff = object$coeffs, beta0 = object$beta0hat)
29-
30-
if(se.fit){
31-
32-
newU <- newdataK %*% mult_diag(object$U, 1/object$D)
33-
# todo: could move to cpp if slow
34-
35-
partiallogit <- partial_logit(newU, object$dhat, object$beta0hat)
36-
# each row is dy/dd with dy/db in the last column
37-
deriv.logit <- cbind(t(mult_diag(t(newU), partiallogit)),
38-
partiallogit)
39-
40-
vcov.fit <- tcrossprod(deriv.logit %*% object$vcov.db0, deriv.logit)
41-
se.fit <- matrix(sqrt(diag(vcov.fit)),ncol=1)
42-
} else {
43-
vcov.fit <- se.fit <- deriv.logit <- NULL
44-
}
45-
46-
} else if (object$loss == "leastsquares") {
47-
48-
yfitted <- newdataK %*% object$coeffs
49-
50-
# ses for fitted
51-
if(se.fit){
52-
# transform to variance of c's on standarized scale
53-
vcov.c.raw <- object$vcov.c * as.vector((1/var(object$y)))
54-
vcov.fitted <- tcrossprod(newdataK%*%vcov.c.raw,newdataK)
55-
vcov.fit <- (apply(object$y,2,sd)^2)*vcov.fitted
56-
se.fit <- matrix(sqrt(diag(vcov.fit)),ncol=1)
57-
} else {
58-
vcov.fit <- se.fit <- NULL
59-
}
60-
61-
# bring back to original scale
62-
yfitted <- (yfitted * apply(object$y,2,sd))+mean(object$y)
63-
64-
}
65-
fit <- list(fit = yfitted,
66-
se.fit = se.fit, vcov.fit = vcov.fit,
67-
newdataK = newdataK)
68-
69-
if(object$loss == "logistic") {
70-
fit$deriv.logit <- deriv.logit
71-
}
72-
73-
return(fit)
74-
}
1+
## This file contains functions for predicting outcome variables
2+
## Functions:
3+
## predict.krls2
4+
## logistic
5+
6+
7+
## a predict function for class 'krls2'
8+
#' @export
9+
predict.krls2 <- function(object, newdata, se.fit = FALSE, ...) {
10+
11+
if(is.null(object$vcov.c) & se.fit)
12+
stop("se.fit requires the vcov.c. Use summary object with predict if you want stand error for the fits")
13+
# todo: we should probably call inference here if they want it?
14+
15+
newdata <- as.matrix(newdata)
16+
if (ncol(object$X) != ncol(newdata)) {
17+
stop("ncol(newdata) differs from ncol(X) from fitted krls object")
18+
}
19+
20+
newdataK <- newKernel(X = object$X, newData = newdata, whichkernel = object$kernel, b = object$b)
21+
22+
if (object$loss == "logistic") {
23+
fit <- predict_logistic(
24+
newdataK = newdataK,
25+
dhat = object$dhat,
26+
coeffs = object$coeffs,
27+
beta0hat = object$beta0hat,
28+
U = object$U,
29+
D = object$D,
30+
vcov.d = object$vcov.d,
31+
se.fit = se.fit
32+
)
33+
} else if (object$loss == "leastsquares") {
34+
fit <- predict_leastsquares(
35+
newdataK = newdataK,
36+
y = object$y,
37+
coeffs = object$coeffs,
38+
vcov.c = object$vcov.c,
39+
se.fit = se.fit
40+
)
41+
}
42+
43+
return(fit)
44+
}
45+
46+
predict_logistic <- function(newdataK, dhat, coeffs, beta0hat, U, D, vcov.d, se.fit) {
47+
48+
yfitted <- logistic(K = newdataK, coeff = coeffs, beta0 = beta0hat)
49+
50+
if(se.fit){
51+
52+
newU <- newdataK %*% mult_diag(U, 1 / D)
53+
# todo: could move to cpp if slow
54+
55+
partiallogit <- partial_logit(newU, dhat, beta0hat)
56+
# each row is dy/dd with dy/db in the last column
57+
deriv.logit <- cbind(t(mult_diag(t(newU), partiallogit)),
58+
partiallogit)
59+
60+
vcov.fit <- tcrossprod(deriv.logit %*% vcov.d, deriv.logit)
61+
se.fit <- matrix(sqrt(diag(vcov.fit)),ncol=1)
62+
} else {
63+
vcov.fit <- se.fit <- deriv.logit <- NULL
64+
}
65+
66+
fit <- list(
67+
fit = yfitted,
68+
se.fit = se.fit,
69+
vcov.fit = vcov.fit,
70+
newdataK = newdataK,
71+
deriv.logit = deriv.logit
72+
)
73+
}
74+
75+
predict_leastsquares <- function(newdataK, coeffs, vcov.c, y, se.fit) {
76+
77+
yfitted <- newdataK %*% coeffs
78+
79+
# ses for fitted
80+
if (se.fit) {
81+
# transform to variance of c's on standarized scale
82+
vcov.c.raw <- vcov.c * as.vector((1 / var(y)))
83+
vcov.fitted <-
84+
tcrossprod(newdataK %*% vcov.c.raw, newdataK)
85+
vcov.fit <- (apply(y, 2, sd) ^ 2) * vcov.fitted
86+
se.fit <- matrix(sqrt(diag(vcov.fit)), ncol = 1)
87+
} else {
88+
vcov.fit <- se.fit <- NULL
89+
}
90+
91+
# bring back to original scale
92+
yfitted <- (yfitted * apply(y, 2, sd)) + mean(y)
93+
94+
fit <- list(
95+
fit = yfitted,
96+
se.fit = se.fit,
97+
vcov.fit = vcov.fit,
98+
newdataK = newdataK
99+
)
100+
101+
}

Diff for: ‎R/summary.krls2.R

+2-2
Original file line numberDiff line numberDiff line change
@@ -25,8 +25,8 @@ summary.krls2 <- function(object,
2525
object <- inference.krls2(object, ...)
2626
}
2727
# average marginal effects
28-
est <- t(object$avgderivatives)
29-
se <- sqrt((object$var.avgderivatives))
28+
est <- object$avgderivatives
29+
se <- sqrt(object$var.avgderivatives)
3030
tval <- est/(se)
3131
avgcoefficients <- cbind(est, se, tval, 2 * pt(abs(tval),n-d, lower.tail = FALSE))
3232
colnames(avgcoefficients) <- c("Est", "Std. Error", "t value", "Pr(>|t|)")

Diff for: ‎src/RcppExports.cpp

+43-58
Original file line numberDiff line numberDiff line change
@@ -6,29 +6,6 @@
66

77
using namespace Rcpp;
88

9-
// mult_diag
10-
arma::mat mult_diag(const arma::mat& x, const arma::vec& d);
11-
RcppExport SEXP _KRLS2_mult_diag(SEXP xSEXP, SEXP dSEXP) {
12-
BEGIN_RCPP
13-
Rcpp::RObject rcpp_result_gen;
14-
Rcpp::RNGScope rcpp_rngScope_gen;
15-
Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
16-
Rcpp::traits::input_parameter< const arma::vec& >::type d(dSEXP);
17-
rcpp_result_gen = Rcpp::wrap(mult_diag(x, d));
18-
return rcpp_result_gen;
19-
END_RCPP
20-
}
21-
// trace_mat
22-
double trace_mat(const arma::mat& x);
23-
RcppExport SEXP _KRLS2_trace_mat(SEXP xSEXP) {
24-
BEGIN_RCPP
25-
Rcpp::RObject rcpp_result_gen;
26-
Rcpp::RNGScope rcpp_rngScope_gen;
27-
Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
28-
rcpp_result_gen = Rcpp::wrap(trace_mat(x));
29-
return rcpp_result_gen;
30-
END_RCPP
31-
}
329
// krls_gr_trunc
3310
arma::vec krls_gr_trunc(const arma::mat& U, const arma::vec& D, const arma::vec& y, const arma::vec& w, const arma::vec& fitted, const arma::vec& dhat, const double& lambda);
3411
RcppExport SEXP _KRLS2_krls_gr_trunc(SEXP USEXP, SEXP DSEXP, SEXP ySEXP, SEXP wSEXP, SEXP fittedSEXP, SEXP dhatSEXP, SEXP lambdaSEXP) {
@@ -121,6 +98,29 @@ BEGIN_RCPP
12198
return rcpp_result_gen;
12299
END_RCPP
123100
}
101+
// mult_diag
102+
arma::mat mult_diag(const arma::mat& x, const arma::vec& d);
103+
RcppExport SEXP _KRLS2_mult_diag(SEXP xSEXP, SEXP dSEXP) {
104+
BEGIN_RCPP
105+
Rcpp::RObject rcpp_result_gen;
106+
Rcpp::RNGScope rcpp_rngScope_gen;
107+
Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
108+
Rcpp::traits::input_parameter< const arma::vec& >::type d(dSEXP);
109+
rcpp_result_gen = Rcpp::wrap(mult_diag(x, d));
110+
return rcpp_result_gen;
111+
END_RCPP
112+
}
113+
// trace_mat
114+
double trace_mat(const arma::mat& x);
115+
RcppExport SEXP _KRLS2_trace_mat(SEXP xSEXP) {
116+
BEGIN_RCPP
117+
Rcpp::RObject rcpp_result_gen;
118+
Rcpp::RNGScope rcpp_rngScope_gen;
119+
Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
120+
rcpp_result_gen = Rcpp::wrap(trace_mat(x));
121+
return rcpp_result_gen;
122+
END_RCPP
123+
}
124124
// euc_dist
125125
double euc_dist(const arma::rowvec& x1, const arma::rowvec& x2);
126126
RcppExport SEXP _KRLS2_euc_dist(SEXP x1SEXP, SEXP x2SEXP) {
@@ -171,6 +171,23 @@ BEGIN_RCPP
171171
return rcpp_result_gen;
172172
END_RCPP
173173
}
174+
// pwmfx
175+
arma::vec pwmfx(const arma::mat& k, const arma::vec& x, const arma::vec& coefhat, const Rcpp::Nullable<Rcpp::NumericMatrix>& vcovc_mat, const arma::vec& p, const arma::vec& p2, const double& b);
176+
RcppExport SEXP _KRLS2_pwmfx(SEXP kSEXP, SEXP xSEXP, SEXP coefhatSEXP, SEXP vcovc_matSEXP, SEXP pSEXP, SEXP p2SEXP, SEXP bSEXP) {
177+
BEGIN_RCPP
178+
Rcpp::RObject rcpp_result_gen;
179+
Rcpp::RNGScope rcpp_rngScope_gen;
180+
Rcpp::traits::input_parameter< const arma::mat& >::type k(kSEXP);
181+
Rcpp::traits::input_parameter< const arma::vec& >::type x(xSEXP);
182+
Rcpp::traits::input_parameter< const arma::vec& >::type coefhat(coefhatSEXP);
183+
Rcpp::traits::input_parameter< const Rcpp::Nullable<Rcpp::NumericMatrix>& >::type vcovc_mat(vcovc_matSEXP);
184+
Rcpp::traits::input_parameter< const arma::vec& >::type p(pSEXP);
185+
Rcpp::traits::input_parameter< const arma::vec& >::type p2(p2SEXP);
186+
Rcpp::traits::input_parameter< const double& >::type b(bSEXP);
187+
rcpp_result_gen = Rcpp::wrap(pwmfx(k, x, coefhat, vcovc_mat, p, p2, b));
188+
return rcpp_result_gen;
189+
END_RCPP
190+
}
174191
// solve_for_d_ls
175192
Rcpp::List solve_for_d_ls(const arma::vec& y, const arma::mat& U, const arma::vec& D, const double& lambda);
176193
RcppExport SEXP _KRLS2_solve_for_d_ls(SEXP ySEXP, SEXP USEXP, SEXP DSEXP, SEXP lambdaSEXP) {
@@ -200,55 +217,23 @@ BEGIN_RCPP
200217
return rcpp_result_gen;
201218
END_RCPP
202219
}
203-
// pwmfx
204-
arma::mat pwmfx(const arma::mat& k, const arma::mat& x, const arma::vec& coefhat, const arma::mat& vcovc, const arma::vec& p, const double& b);
205-
RcppExport SEXP _KRLS2_pwmfx(SEXP kSEXP, SEXP xSEXP, SEXP coefhatSEXP, SEXP vcovcSEXP, SEXP pSEXP, SEXP bSEXP) {
206-
BEGIN_RCPP
207-
Rcpp::RObject rcpp_result_gen;
208-
Rcpp::RNGScope rcpp_rngScope_gen;
209-
Rcpp::traits::input_parameter< const arma::mat& >::type k(kSEXP);
210-
Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
211-
Rcpp::traits::input_parameter< const arma::vec& >::type coefhat(coefhatSEXP);
212-
Rcpp::traits::input_parameter< const arma::mat& >::type vcovc(vcovcSEXP);
213-
Rcpp::traits::input_parameter< const arma::vec& >::type p(pSEXP);
214-
Rcpp::traits::input_parameter< const double& >::type b(bSEXP);
215-
rcpp_result_gen = Rcpp::wrap(pwmfx(k, x, coefhat, vcovc, p, b));
216-
return rcpp_result_gen;
217-
END_RCPP
218-
}
219-
// pwmfx_novar
220-
arma::mat pwmfx_novar(const arma::mat& k, const arma::mat& x, const arma::vec& coefhat, const arma::vec& p, const double& b);
221-
RcppExport SEXP _KRLS2_pwmfx_novar(SEXP kSEXP, SEXP xSEXP, SEXP coefhatSEXP, SEXP pSEXP, SEXP bSEXP) {
222-
BEGIN_RCPP
223-
Rcpp::RObject rcpp_result_gen;
224-
Rcpp::RNGScope rcpp_rngScope_gen;
225-
Rcpp::traits::input_parameter< const arma::mat& >::type k(kSEXP);
226-
Rcpp::traits::input_parameter< const arma::mat& >::type x(xSEXP);
227-
Rcpp::traits::input_parameter< const arma::vec& >::type coefhat(coefhatSEXP);
228-
Rcpp::traits::input_parameter< const arma::vec& >::type p(pSEXP);
229-
Rcpp::traits::input_parameter< const double& >::type b(bSEXP);
230-
rcpp_result_gen = Rcpp::wrap(pwmfx_novar(k, x, coefhat, p, b));
231-
return rcpp_result_gen;
232-
END_RCPP
233-
}
234220

235221
static const R_CallMethodDef CallEntries[] = {
236-
{"_KRLS2_mult_diag", (DL_FUNC) &_KRLS2_mult_diag, 2},
237-
{"_KRLS2_trace_mat", (DL_FUNC) &_KRLS2_trace_mat, 1},
238222
{"_KRLS2_krls_gr_trunc", (DL_FUNC) &_KRLS2_krls_gr_trunc, 7},
239223
{"_KRLS2_krls_hess_trunc_inv", (DL_FUNC) &_KRLS2_krls_hess_trunc_inv, 4},
240224
{"_KRLS2_krlogit_fn_trunc", (DL_FUNC) &_KRLS2_krlogit_fn_trunc, 6},
241225
{"_KRLS2_krlogit_gr_trunc", (DL_FUNC) &_KRLS2_krlogit_gr_trunc, 6},
242226
{"_KRLS2_partial_logit", (DL_FUNC) &_KRLS2_partial_logit, 3},
243227
{"_KRLS2_krlogit_hess_trunc_inv", (DL_FUNC) &_KRLS2_krlogit_hess_trunc_inv, 6},
228+
{"_KRLS2_mult_diag", (DL_FUNC) &_KRLS2_mult_diag, 2},
229+
{"_KRLS2_trace_mat", (DL_FUNC) &_KRLS2_trace_mat, 1},
244230
{"_KRLS2_euc_dist", (DL_FUNC) &_KRLS2_euc_dist, 2},
245231
{"_KRLS2_kern_gauss_1d", (DL_FUNC) &_KRLS2_kern_gauss_1d, 3},
246232
{"_KRLS2_kern_gauss", (DL_FUNC) &_KRLS2_kern_gauss, 2},
247233
{"_KRLS2_new_gauss_kern", (DL_FUNC) &_KRLS2_new_gauss_kern, 3},
234+
{"_KRLS2_pwmfx", (DL_FUNC) &_KRLS2_pwmfx, 7},
248235
{"_KRLS2_solve_for_d_ls", (DL_FUNC) &_KRLS2_solve_for_d_ls, 4},
249236
{"_KRLS2_solve_for_d_ls_w", (DL_FUNC) &_KRLS2_solve_for_d_ls_w, 5},
250-
{"_KRLS2_pwmfx", (DL_FUNC) &_KRLS2_pwmfx, 6},
251-
{"_KRLS2_pwmfx_novar", (DL_FUNC) &_KRLS2_pwmfx_novar, 5},
252237
{NULL, NULL, 0}
253238
};
254239

Diff for: ‎src/krls.cpp

+4
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
// [[Rcpp::depends(RcppArmadillo)]]
2+
// [[Rcpp::plugins(cpp11)]]
3+
4+
#define ARMA_64BIT_WORD 1

Diff for: ‎src/krls_grad_hess.cpp

+111
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
#include <RcppArmadillo.h>
2+
#include "krls_helpers.h"
3+
using namespace Rcpp;
4+
5+
// [[Rcpp::export]]
6+
arma::vec krls_gr_trunc(
7+
const arma::mat& U,
8+
const arma::vec& D,
9+
const arma::vec& y,
10+
const arma::vec& w,
11+
const arma::vec& fitted,
12+
const arma::vec& dhat,
13+
const double& lambda) {
14+
15+
arma::vec score = -2 * (mult_diag(U.t(), w) * (y - fitted) + lambda * (accu(w) / y.n_elem) * (dhat / D));
16+
17+
return score;
18+
}
19+
20+
// Hessian used for sandwich estimator for KRLS
21+
// [[Rcpp::export]]
22+
arma::mat krls_hess_trunc_inv(const arma::mat& U,
23+
const arma::vec& D,
24+
const arma::vec& w,
25+
const double& lambda) {
26+
27+
arma::mat hess = 2 * (mult_diag(U.t(), w) * U + arma::diagmat(lambda * (1 / D)));
28+
29+
return arma::inv_sympd(hess);
30+
31+
}
32+
33+
// [[Rcpp::export]]
34+
double krlogit_fn_trunc(const arma::vec& par,
35+
const arma::mat& U,
36+
const arma::vec& D,
37+
const arma::vec& y,
38+
const arma::vec& w,
39+
const double& lambda) {
40+
41+
arma::vec coef = par.subvec(0, par.n_elem - 2);
42+
double beta0 = par(par.n_elem-1);
43+
arma::mat Ud = U * coef;
44+
45+
double ret = accu(w % (y % log(1 + exp(-(beta0 + Ud))) +
46+
(1 - y) % log(1 + exp(beta0 + Ud))) +
47+
lambda * arma::as_scalar(coef.t() * ((1.0/D) % coef)));
48+
49+
return ret;
50+
51+
}
52+
53+
// [[Rcpp::export]]
54+
arma::vec krlogit_gr_trunc(const arma::vec& par,
55+
const arma::mat& U,
56+
const arma::vec& D,
57+
const arma::vec& y,
58+
const arma::vec& w,
59+
const double& lambda) {
60+
61+
arma::vec coef = par.subvec(0, par.n_elem - 2);
62+
double beta0 = par(par.n_elem-1);
63+
arma::vec resid = y - (1 / (1 + exp(-U * coef - beta0)));
64+
65+
arma::vec ret(par.n_elem);
66+
67+
ret.subvec(0, par.n_elem - 2) = -U.t() * (w % resid) + 2 * y.n_elem * (lambda / D) % coef;
68+
ret(par.n_elem - 1) = -accu(w % resid);
69+
70+
return ret;
71+
}
72+
73+
// Generates (exp(-Kc - b)) / (1+exp(-Kc-b))^2 used for the krlogit hessian
74+
// and for for predicted value SEs and first difference SEs
75+
// Should work with either truncated data or full data
76+
// [[Rcpp::export]]
77+
arma::vec partial_logit(const arma::mat& K,
78+
const arma::vec& coef,
79+
const double& beta0) {
80+
arma::vec ret = exp(-K * coef - beta0) / pow((1 + exp(-K * coef - beta0)), 2);
81+
82+
return ret;
83+
}
84+
85+
// [[Rcpp::export]]
86+
arma::mat krlogit_hess_trunc_inv(const arma::vec& par,
87+
const arma::mat& U,
88+
const arma::vec& D,
89+
const arma::vec& y,
90+
const arma::vec& w,
91+
const double& lambda) {
92+
93+
arma::vec coef = par.subvec(0, par.n_elem - 2);
94+
double beta0 = par(par.n_elem-1);
95+
arma::vec meat = w % partial_logit(U, coef, beta0);
96+
97+
arma::mat ret(par.n_elem, par.n_elem);
98+
99+
arma::mat dcdc = mult_diag(U.t(), meat) * U + diagmat(2 * y.n_elem * (lambda / D));
100+
101+
arma::vec dcdb = U.t() * meat;
102+
103+
double dbdb = accu(meat);
104+
105+
ret.submat(0, 0, coef.n_elem - 1, coef.n_elem - 1) = dcdc;
106+
ret.submat(coef.n_elem, 0, coef.n_elem, coef.n_elem - 1) = dcdb.t();
107+
ret.submat(0, coef.n_elem, coef.n_elem - 1, coef.n_elem) = dcdb;
108+
ret(coef.n_elem, coef.n_elem) = dbdb;
109+
110+
return arma::inv_sympd(ret);
111+
}

Diff for: ‎src/krls_helpers.cpp

+4-305
Original file line numberDiff line numberDiff line change
@@ -1,324 +1,23 @@
1-
// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
2-
3-
// [[Rcpp::depends(RcppArmadillo)]]
4-
// [[Rcpp::plugins(cpp11)]]
5-
6-
#define ARMA_64BIT_WORD 1
71
#include <RcppArmadillo.h>
2+
#include "krls_helpers.h"
83
using namespace Rcpp;
94

105
//----------
116
// Helper functions
127
//----------
138

14-
//' @export
159
// [[Rcpp::export]]
1610
arma::mat mult_diag(const arma::mat& x, const arma::vec& d) {
17-
11+
1812
arma::mat out(x.n_rows, x.n_cols);
1913
for (unsigned j = 0; j < x.n_cols; ++j) {
2014
out.col(j) = x.col(j) * d(j);
2115
}
22-
16+
2317
return out;
2418
}
2519

26-
//' @export
2720
// [[Rcpp::export]]
2821
double trace_mat(const arma::mat& x) {
2922
return arma::trace(x);
30-
}
31-
32-
//' @export
33-
// [[Rcpp::export]]
34-
arma::vec krls_gr_trunc(
35-
const arma::mat& U,
36-
const arma::vec& D,
37-
const arma::vec& y,
38-
const arma::vec& w,
39-
const arma::vec& fitted,
40-
const arma::vec& dhat,
41-
const double& lambda) {
42-
43-
arma::vec score = -2 * (mult_diag(U.t(), w) * (y - fitted) + lambda * (accu(w) / y.n_elem) * (dhat / D));
44-
45-
return score;
46-
}
47-
48-
// Hessian used for sandwich estimator for KRLS
49-
//' @export
50-
// [[Rcpp::export]]
51-
arma::mat krls_hess_trunc_inv(const arma::mat& U,
52-
const arma::vec& D,
53-
const arma::vec& w,
54-
const double& lambda) {
55-
56-
arma::mat hess = 2 * (mult_diag(U.t(), w) * U + arma::diagmat(lambda * (1 / D)));
57-
58-
return arma::inv_sympd(hess);
59-
60-
}
61-
62-
63-
//' @export
64-
// [[Rcpp::export]]
65-
double krlogit_fn_trunc(const arma::vec& par,
66-
const arma::mat& U,
67-
const arma::vec& D,
68-
const arma::vec& y,
69-
const arma::vec& w,
70-
const double& lambda) {
71-
72-
arma::vec coef = par.subvec(0, par.n_elem - 2);
73-
double beta0 = par(par.n_elem-1);
74-
arma::mat Ud = U * coef;
75-
76-
double ret = accu(w % (y % log(1 + exp(-(beta0 + Ud))) +
77-
(1 - y) % log(1 + exp(beta0 + Ud))) +
78-
lambda * arma::as_scalar(coef.t() * ((1.0/D) % coef)));
79-
80-
return ret;
81-
82-
}
83-
84-
//' @export
85-
// [[Rcpp::export]]
86-
arma::vec krlogit_gr_trunc(const arma::vec& par,
87-
const arma::mat& U,
88-
const arma::vec& D,
89-
const arma::vec& y,
90-
const arma::vec& w,
91-
const double& lambda) {
92-
93-
arma::vec coef = par.subvec(0, par.n_elem - 2);
94-
double beta0 = par(par.n_elem-1);
95-
arma::vec resid = y - (1 / (1 + exp(-U * coef - beta0)));
96-
97-
arma::vec ret(par.n_elem);
98-
99-
ret.subvec(0, par.n_elem - 2) = -U.t() * (w % resid) + 2 * y.n_elem * (lambda / D) % coef;
100-
ret(par.n_elem - 1) = -accu(w % resid);
101-
102-
return ret;
103-
}
104-
105-
// Generates (exp(-Kc - b)) / (1+exp(-Kc-b))^2 used for the krlogit hessian
106-
// and for for predicted value SEs and first difference SEs
107-
// Should work with either truncated data or full data
108-
//' @export
109-
// [[Rcpp::export]]
110-
arma::vec partial_logit(const arma::mat& K,
111-
const arma::vec& coef,
112-
const double& beta0) {
113-
arma::vec ret = exp(-K * coef - beta0) / pow((1 + exp(-K * coef - beta0)), 2);
114-
115-
return ret;
116-
}
117-
118-
//' @export
119-
// [[Rcpp::export]]
120-
arma::mat krlogit_hess_trunc_inv(const arma::vec& par,
121-
const arma::mat& U,
122-
const arma::vec& D,
123-
const arma::vec& y,
124-
const arma::vec& w,
125-
const double& lambda) {
126-
127-
arma::vec coef = par.subvec(0, par.n_elem - 2);
128-
double beta0 = par(par.n_elem-1);
129-
arma::vec meat = w % partial_logit(U, coef, beta0);
130-
131-
arma::mat ret(par.n_elem, par.n_elem);
132-
133-
arma::mat dcdc = mult_diag(U.t(), meat) * U + diagmat(2 * y.n_elem * (lambda / D));
134-
135-
arma::vec dcdb = U.t() * meat;
136-
137-
double dbdb = accu(meat);
138-
139-
ret.submat(0, 0, coef.n_elem - 1, coef.n_elem - 1) = dcdc;
140-
ret.submat(coef.n_elem, 0, coef.n_elem, coef.n_elem - 1) = dcdb.t();
141-
ret.submat(0, coef.n_elem, coef.n_elem - 1, coef.n_elem) = dcdb;
142-
ret(coef.n_elem, coef.n_elem) = dbdb;
143-
144-
return arma::inv_sympd(ret);
145-
}
146-
147-
// ----------
148-
// GAUSSIAN KERNELS
149-
// ----------
150-
151-
// Euclidean distance function
152-
153-
//' @export
154-
// [[Rcpp::export]]
155-
double euc_dist(const arma::rowvec& x1, const arma::rowvec& x2) {
156-
double out = 0.0;
157-
unsigned n = x1.n_elem;
158-
159-
for (unsigned i = 0; i < n; ++i) {
160-
out += pow(x1(i) - x2(i), 2);
161-
}
162-
return sqrt(out);
163-
}
164-
165-
// Gaussian kernel between two row vectors
166-
167-
//' @export
168-
// [[Rcpp::export]]
169-
double kern_gauss_1d(const arma::rowvec& x1, const arma::rowvec& x2, const double& b)
170-
{
171-
return exp(-pow(euc_dist(x1, x2), 2) / (b));
172-
}
173-
174-
// Gaussian kernel matrix for a matrix with itself
175-
// b is sigma, or 2P by default
176-
177-
//' @export
178-
// [[Rcpp::export]]
179-
arma::mat kern_gauss(const arma::mat& x, const double& b)
180-
{
181-
unsigned n = x.n_rows;
182-
double val;
183-
// Diagonal will remain ones
184-
arma::mat out(n, n, arma::fill::ones);
185-
186-
for (unsigned i = 0; i < n; ++i) {
187-
188-
for (unsigned j = i + 1; j < n; ++j) {
189-
val = kern_gauss_1d(x.row(i), x.row(j), b);
190-
out(i, j) = val;
191-
out(j, i) = val;
192-
}
193-
194-
}
195-
return out;
196-
}
197-
198-
// Kernel matrix for distance between two matrices
199-
200-
//' @export
201-
// [[Rcpp::export]]
202-
arma::mat new_gauss_kern(const arma::mat& newx, const arma::mat& oldx, const double& b) {
203-
unsigned n1 = newx.n_rows;
204-
unsigned n2 = oldx.n_rows;
205-
arma::mat out(n1, n2);
206-
207-
for (unsigned i = 0; i < n1; ++i) {
208-
for (unsigned j = 0; j < n2; ++j) {
209-
out(i, j) = kern_gauss_1d(newx.row(i), oldx.row(j), b);
210-
}
211-
212-
}
213-
return out;
214-
}
215-
216-
217-
//' @export
218-
// [[Rcpp::export]]
219-
Rcpp::List solve_for_d_ls(const arma::vec& y,
220-
const arma::mat& U,
221-
const arma::vec& D,
222-
const double& lambda) {
223-
224-
arma::vec Ginv = 1 / (1 + lambda / D);
225-
226-
arma::vec dhat = Ginv % (U.t() * y);
227-
// This is the same as above
228-
//arma::vec tempLoss = (y - U * coeffs) / diagvec(arma::eye(y.n_elem, y.n_elem) - mult_diag(U, Ginv) * U.t());
229-
230-
arma::vec tempLoss = (y - U * mult_diag(U, Ginv).t() * y) / diagvec(arma::eye(y.n_elem, y.n_elem) - U * mult_diag(U, Ginv).t());
231-
double Le = as_scalar(tempLoss.t() * tempLoss);
232-
233-
return Rcpp::List::create(Rcpp::Named("dhat") = dhat,
234-
Rcpp::Named("Le") = Le);
235-
}
236-
237-
//' @export
238-
// [[Rcpp::export]]
239-
Rcpp::List solve_for_d_ls_w(const arma::vec& y,
240-
const arma::mat& U,
241-
const arma::vec& D,
242-
const arma::vec& w,
243-
const double& lambda) {
244-
245-
arma::mat Uw = mult_diag(U.t(), w);
246-
arma::mat Ginv = arma::inv_sympd(Uw * U + diagmat(lambda / D));
247-
248-
arma::vec dhat = Ginv * (Uw * y);
249-
// This is the same as above
250-
//arma::vec tempLoss = (y - U * coeffs) / diagvec(arma::eye(y.n_elem, y.n_elem) - mult_diag(U, Ginv) * U.t());
251-
arma::mat UGinvUw = U * Ginv * Uw;
252-
253-
arma::vec tempLoss = (y - UGinvUw * y) / diagvec(arma::eye(y.n_elem, y.n_elem) - UGinvUw);
254-
double Le = as_scalar(tempLoss.t() * tempLoss);
255-
256-
return Rcpp::List::create(Rcpp::Named("dhat") = dhat,
257-
Rcpp::Named("Le") = Le);
258-
}
259-
260-
// Compute pointwise marginal effects and the var avg pwmfx
261-
//' @export
262-
// [[Rcpp::export]]
263-
arma::mat pwmfx(const arma::mat& k,
264-
const arma::mat& x,
265-
const arma::vec& coefhat,
266-
const arma::mat& vcovc,
267-
const arma::vec& p,
268-
const double& b)
269-
{
270-
271-
double n = x.n_rows;
272-
arma::mat out(n + 1, x.n_cols);
273-
arma::mat distmat(n, n);
274-
arma::mat distk(n, n);
275-
arma::vec p2 = pow(p, 2);
276-
double val;
277-
278-
for (unsigned j = 0; j < x.n_cols; ++j) {
279-
for (unsigned i = 0; i < n; ++i) {
280-
val = 0;
281-
for (unsigned i2 = 0; i2 < n; ++i2) {
282-
distmat(i, i2) = x(i, j) - x(i2, j);
283-
284-
val += coefhat(i2) * k(i, i2) * distmat(i, i2);
285-
}
286-
287-
out(i, j) = - (p(i) / b) * val;
288-
}
289-
290-
distk = k % distmat;
291-
out(n, j) = 1 / pow(b * n, 2) * accu(p2.t() * distk.t() * vcovc * distk);
292-
}
293-
294-
return out;
295-
}
296-
297-
298-
// Compute pointwise marginal effects but not var avg pwmfx
299-
//' @export
300-
// [[Rcpp::export]]
301-
arma::mat pwmfx_novar(const arma::mat& k,
302-
const arma::mat& x,
303-
const arma::vec& coefhat,
304-
const arma::vec& p,
305-
const double& b)
306-
{
307-
308-
double n = x.n_rows;
309-
arma::mat out(n, x.n_cols);
310-
double val;
311-
312-
for (unsigned j = 0; j < x.n_cols; ++j) {
313-
for (unsigned i = 0; i < n; ++i) {
314-
val = 0;
315-
for (unsigned i2 = 0; i2 < n; ++i2) {
316-
val += coefhat(i2) * k(i, i2) * (x(i, j) - x(i2, j));
317-
}
318-
319-
out(i, j) = - (p(i) / b) * val;
320-
}
321-
}
322-
323-
return out;
324-
}
23+
}

Diff for: ‎src/krls_helpers.h

+8
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
#ifndef KRLS_HELPER_H
2+
#define KRLS_HELPER_H
3+
4+
#include <RcppArmadillo.h>
5+
arma::mat mult_diag(const arma::mat& x, const arma::vec& d);
6+
double trace_mat(const arma::mat& x);
7+
8+
#endif

Diff for: ‎src/krls_kernels.cpp

+65
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
#include <RcppArmadillo.h>
2+
using namespace Rcpp;
3+
4+
// ----------
5+
// GAUSSIAN KERNELS
6+
// ----------
7+
8+
// Euclidean distance function
9+
10+
// [[Rcpp::export]]
11+
double euc_dist(const arma::rowvec& x1, const arma::rowvec& x2) {
12+
double out = 0.0;
13+
unsigned n = x1.n_elem;
14+
15+
for (unsigned i = 0; i < n; ++i) {
16+
out += pow(x1(i) - x2(i), 2);
17+
}
18+
return sqrt(out);
19+
}
20+
21+
// Gaussian kernel between two row vectors
22+
23+
// [[Rcpp::export]]
24+
double kern_gauss_1d(const arma::rowvec& x1, const arma::rowvec& x2, const double& b)
25+
{
26+
return exp(-pow(euc_dist(x1, x2), 2) / (b));
27+
}
28+
29+
// Gaussian kernel matrix for a matrix with itself
30+
// b is sigma, or 2P by default
31+
// [[Rcpp::export]]
32+
arma::mat kern_gauss(const arma::mat& x, const double& b)
33+
{
34+
unsigned n = x.n_rows;
35+
double val;
36+
// Diagonal will remain ones
37+
arma::mat out(n, n, arma::fill::ones);
38+
39+
for (unsigned i = 0; i < n; ++i) {
40+
41+
for (unsigned j = i + 1; j < n; ++j) {
42+
val = kern_gauss_1d(x.row(i), x.row(j), b);
43+
out(i, j) = val;
44+
out(j, i) = val;
45+
}
46+
47+
}
48+
return out;
49+
}
50+
51+
// Kernel matrix for distance between two matrices
52+
// [[Rcpp::export]]
53+
arma::mat new_gauss_kern(const arma::mat& newx, const arma::mat& oldx, const double& b) {
54+
unsigned n1 = newx.n_rows;
55+
unsigned n2 = oldx.n_rows;
56+
arma::mat out(n1, n2);
57+
58+
for (unsigned i = 0; i < n1; ++i) {
59+
for (unsigned j = 0; j < n2; ++j) {
60+
out(i, j) = kern_gauss_1d(newx.row(i), oldx.row(j), b);
61+
}
62+
63+
}
64+
return out;
65+
}

Diff for: ‎src/krls_pwmfx.cpp

+41
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
#include <RcppArmadillo.h>
2+
using namespace Rcpp;
3+
4+
// Compute pointwise marginal effects and the var avg pwmfx for one variable
5+
// [[Rcpp::export]]
6+
arma::vec pwmfx(const arma::mat& k,
7+
const arma::vec& x,
8+
const arma::vec& coefhat,
9+
const Rcpp::Nullable<Rcpp::NumericMatrix>& vcovc_mat,
10+
const arma::vec& p,
11+
const arma::vec& p2,
12+
const double& b)
13+
{
14+
15+
arma::mat vcovc;
16+
if (vcovc_mat.isNotNull()) {
17+
vcovc = Rcpp::as<arma::mat>(vcovc_mat);
18+
}
19+
double n = x.n_elem;
20+
arma::vec out(n + 1);
21+
arma::mat distmat(n, n);
22+
double val = 0.0;
23+
24+
for (unsigned i = 0; i < n; ++i) {
25+
val = 0.0;
26+
for (unsigned i2 = 0; i2 < n; ++i2) {
27+
distmat(i, i2) = x(i) - x(i2);
28+
29+
val += coefhat(i2) * k(i, i2) * distmat(i, i2);
30+
}
31+
32+
out(i) = - (p(i) / b) * val;
33+
}
34+
35+
if (vcovc_mat.isNotNull()) {
36+
arma::mat distk = k % distmat;
37+
out(n) = 1 / pow(b * n, 2) * accu(p2.t() * distk.t() * vcovc * distk);
38+
}
39+
40+
return out;
41+
}

Diff for: ‎src/krls_solve.cpp

+44
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
#include <RcppArmadillo.h>
2+
#include "krls_helpers.h"
3+
using namespace Rcpp;
4+
5+
// [[Rcpp::export]]
6+
Rcpp::List solve_for_d_ls(const arma::vec& y,
7+
const arma::mat& U,
8+
const arma::vec& D,
9+
const double& lambda) {
10+
11+
arma::vec Ginv = 1 / (1 + lambda / D);
12+
13+
arma::vec dhat = Ginv % (U.t() * y);
14+
// This is the same as above
15+
//arma::vec tempLoss = (y - U * coeffs) / diagvec(arma::eye(y.n_elem, y.n_elem) - mult_diag(U, Ginv) * U.t());
16+
17+
arma::vec tempLoss = (y - U * mult_diag(U, Ginv).t() * y) / diagvec(arma::eye(y.n_elem, y.n_elem) - U * mult_diag(U, Ginv).t());
18+
double Le = as_scalar(tempLoss.t() * tempLoss);
19+
20+
return Rcpp::List::create(Rcpp::Named("dhat") = dhat,
21+
Rcpp::Named("Le") = Le);
22+
}
23+
24+
// [[Rcpp::export]]
25+
Rcpp::List solve_for_d_ls_w(const arma::vec& y,
26+
const arma::mat& U,
27+
const arma::vec& D,
28+
const arma::vec& w,
29+
const double& lambda) {
30+
31+
arma::mat Uw = mult_diag(U.t(), w);
32+
arma::mat Ginv = arma::inv_sympd(Uw * U + diagmat(lambda / D));
33+
34+
arma::vec dhat = Ginv * (Uw * y);
35+
// This is the same as above
36+
//arma::vec tempLoss = (y - U * coeffs) / diagvec(arma::eye(y.n_elem, y.n_elem) - mult_diag(U, Ginv) * U.t());
37+
arma::mat UGinvUw = U * Ginv * Uw;
38+
39+
arma::vec tempLoss = (y - UGinvUw * y) / diagvec(arma::eye(y.n_elem, y.n_elem) - UGinvUw);
40+
double Le = as_scalar(tempLoss.t() * tempLoss);
41+
42+
return Rcpp::List::create(Rcpp::Named("dhat") = dhat,
43+
Rcpp::Named("Le") = Le);
44+
}

Diff for: ‎tests/testthat/test_old_krls.R

+42
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
# This file is .Rbuildignore'd as it needs the old KRLS package
2+
3+
context("Match KRLS 1.0.0")
4+
5+
n <- 20
6+
X <- cbind(rnorm(n), rbinom(n, 1, 0.5))
7+
y <- rnorm(20)
8+
9+
test_that("", {
10+
11+
# Default sigma (b) is different
12+
oldout <- KRLS::krls(X = X, y = y, sigma = 2*ncol(X), print.level = 3)
13+
newout <- KRLS2::inference.krls2(KRLS2::krls(X = X, y = y, printlevel = 3))
14+
15+
expect_equivalent(
16+
oldout$lambda,
17+
newout$lambda
18+
)
19+
20+
# Lambda and coeffs are different because we use eigtrunc
21+
expect_equal(
22+
oldout$coeffs,
23+
newout$coeffs
24+
)
25+
26+
expect_equal(
27+
oldout$derivatives,
28+
newout$derivatives
29+
)
30+
31+
expect_equivalent(
32+
oldout$avgderivatives,
33+
newout$avgderivatives
34+
)
35+
36+
expect_equivalent(
37+
oldout$var.avgderivatives,
38+
newout$var.avgderivatives * c(1, 2) # old variance of FDs had been multiplied by 2
39+
)
40+
41+
})
42+

Diff for: ‎tests/testthat/test_truncation.R

+4-4
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,9 @@ test_that("Truncation returns roughly same results as no truncation", {
44
# This largely tests if anything goes really wrong in truncation
55
# not that it is a good approximation
66

7-
krlso <- krls(y = mtcars$am, X = mtcars$mpg, truncate = FALSE)
7+
krlso <- krls(y = mtcars$am, X = mtcars[, c("mpg", "wt")], truncate = FALSE)
88

9-
krlso_trunc <- krls(y = mtcars$am, X = mtcars$mpg, epsilon = 0.0000001)
9+
krlso_trunc <- krls(y = mtcars$am, X = mtcars[, c("mpg", "wt")], epsilon = 0.0000001)
1010

1111
expect_true(
1212
ncol(krlso_trunc$U) < ncol(krlso$U)
@@ -16,14 +16,14 @@ test_that("Truncation returns roughly same results as no truncation", {
1616
max(
1717
inference.krls2(krlso)$derivatives -
1818
inference.krls2(krlso_trunc)$derivatives
19-
) < 1e-6
19+
) < 1e-4
2020
)
2121

2222
expect_true(
2323
max(
2424
krlso$fitted -
2525
krlso_trunc$fitted
26-
) < 1e-7
26+
) < 1e-6
2727
)
2828

2929
# Overfitting without truncation with logistic seems to be a serious problem

0 commit comments

Comments
 (0)
Please sign in to comment.