1
#' Class affinity maximum likelihood text scaling model
2
#'
3
#' `textmodel_affinity` implements the maximum likelihood supervised text
4
#' scaling method described in Perry and Benoit (2017).
5
#' @param x the [dfm] or [bootstrap_dfm] object on which the model
6
#'   will be fit.  Does not need to contain only the training documents, since
7
#'   the index of these will be matched automatically.
8
#' @param y vector of training classes/scores associated with each document
9
#'   identified in `data`
10
#' @param exclude a set of words to exclude from the model
11
#' @param smooth a smoothing parameter for class affinities; defaults to 0.5
12
#'   (Jeffreys prior). A plausible alternative would be 1.0 (Laplace prior).
13
#' @param ref_smooth a smoothing parameter for token distributions;
14
#'   defaults to 0.5
15
#' @param verbose logical; if `TRUE` print diagnostic information during
16
#'   fitting.
17
#' @author Patrick Perry and Kenneth Benoit
18
#' @references Perry, P.O. & Benoit, K.R. (2017). Scaling Text with
19
#'   the Class Affinity Model.
20
#'   [arXiv:1710.08963 \[stat.ML\]](http://arxiv.org/abs/1710.08963).
21
#' @examples
22
#' (af <- textmodel_affinity(data_dfm_lbgexample, y = c("L", NA, NA, NA, "R", NA)))
23
#' predict(af)
24
#' predict(af, newdata = data_dfm_lbgexample[6, ])
25
#'
26
#' \dontrun{
27
#' # compute bootstrapped SEs
28
#' dfmat <- quanteda::bootstrap_dfm(data_corpus_dailnoconf1991, n = 10, remove_punct = TRUE)
29
#' textmodel_affinity(dfmat, y = c("Govt", "Opp", "Opp", rep(NA, 55)))
30
#' }
31
#' @export
32
#' @keywords textmodel experimental
33
#' @importFrom methods as
34
#' @importFrom stats sd predict
35
#' @importFrom quanteda dfm_group as.dfm
36
#' @seealso [predict.textmodel_affinity()] for methods of applying a
37
#'   fitted [textmodel_affinity] model object to predict quantities from
38
#'   (other) documents.
39
textmodel_affinity <- function(x, y, exclude = NULL,
40
                               smooth = 0.5, ref_smooth = 0.5,
41
                               verbose = quanteda_options("verbose")) {
42 1
    UseMethod("textmodel_affinity")
43
}
44

45
#' @export
46
textmodel_affinity.default <- function(x, y, exclude = NULL,
47
                                       smooth = 0.5, ref_smooth = 0.5,
48
                                       verbose = quanteda_options("verbose")) {
49 1
    stop(friendly_class_undefined_message(class(x), "textmodel_affinity"))
50
}
51

52

53
#' @export
54
#' @importFrom quanteda colSums rowSums t
55
textmodel_affinity.dfm <- function(x, y, exclude = NULL,
56
                                   smooth = 0.5, ref_smooth = 0.5,
57
                                   verbose = quanteda_options("verbose")) {
58 1
    x <- as.dfm(x)
59 1
    if (!sum(x)) stop(message_error("dfm_empty"))
60

61
    # estimate reference distributions
62 1
    counts <- dfm_group(x, groups = y, force = TRUE)
63

64
    # determine support set
65 1
    support <- colSums(counts) > 0
66

67
    # remove 'exclude' words from support
68 1
    if (!is.null(exclude)) {
69 0
        if (is.character(exclude)) {
70 0
            words <- colnames(x)
71 0
            support[words %in% exclude] <- FALSE
72
        } else {
73 0
            support[exclude] <- FALSE
74
        }
75 0
        counts[, !support] <- 0
76 0
        tots <- rowSums(counts)
77
    }
78

79
    # add Laplace-style smoothing, but only to words in the support
80 1
    counts[, support] <- counts[, support] + ref_smooth
81 1
    tots <- rowSums(counts)
82 1
    probs <- counts / tots
83

84 1
    p <- t(probs) # affinity expects transposed probabilities
85

86 1
    fitted <- affinity(p, t(x), smooth = smooth)
87 1
    result <- list(
88 1
        smooth = fitted$smooth,
89 1
        x = x,
90 1
        y = y,
91 1
        p = p,
92 1
        support = fitted$support,
93 1
        method = "affinity",
94 1
        call = match.call()
95
    )
96 1
    class(result) <- "textmodel_affinity"
97 1
    return(result)
98
}
99

100
#' @export
101
textmodel_affinity.dfm_bootstrap <- function(x, y, exclude = NULL,
102
                                             smooth = 0.5, ref_smooth = 0.5,
103
                                             verbose = quanteda_options("verbose")) {
104 0
    if (verbose)
105 0
        message("Bootstrapping textmodel_affinity for ", ndoc(x[[1]]), " documents:")
106

107
    # compute the model for the original corpus
108 0
    if (verbose)
109 0
        message("   ...computing model for dfm from original texts")
110 0
    fitted <- textmodel_affinity(x[[1]], y, exclude, smooth, ref_smooth)
111 0
    predicted <- predict(fitted)
112

113
    # compute the replicates, save coefficients to an array
114 0
    coeff_replicates <- array(NA, dim = c(dim(coef(predicted)), length(x)))
115 0
    if (verbose)
116 0
        message("   ...computing bootstrapped models and saving coefficients: ", appendLF = FALSE)
117 0
    for (i in seq_len(length(x))) {
118 0
        message(i, " ", appendLF = FALSE)
119 0
        temp <- textmodel_affinity(x[[i]], y, exclude, smooth, ref_smooth)
120 0
        coeff_replicates[, , i] <- coef(predict(temp))
121
    }
122 0
    message("")
123 0
    if (verbose)
124 0
        message("   ...replacing original SEs with bootstrapped SEs")
125
    # replace analytical coefficients with sd of replicates of coefficients
126 0
    predicted$se <- apply(coeff_replicates, c(1, 2), sd)
127 0
    if (verbose)
128 0
        message("   ...finished.")
129 0
    predicted
130
}
131

132
## ============= Internal computation methods ========
133

134
#' Internal function to fit the likelihood scaling mixture model.
135
#'
136
#' Ken recommends you use [textmodel_affinity()] instead.
137
#' @param p word likelihoods within classes, estimated from training data
138
#' @param x term-document matrix for document(s) to be scaled
139
#' @param smooth a misnamed smoothing parameter, either a scalar or a vector
140
#'   equal in length to the number of documents
141
#' @author Patrick Perry
142
#' @return a list with stuff
143
#' @examples
144
#' p <- matrix(c(c(5/6, 0, 1/6), c(0, 4/5, 1/5)), nrow = 3,
145
#'             dimnames = list(c("A", "B", "C"), NULL))
146
#' theta <- c(.2, .8)
147
#' q <- drop(p %*% theta)
148
#' x <- 2 * q
149
#' (fit <- affinity(p, x))
150
#' @keywords internal textmodel
151
#' @export
152
affinity <- function(p, x, smooth = 0.5, verbose = FALSE) {
153 1
    p <- as.matrix(p)
154 1
    ncat <- nrow(p)   # number of features
155 1
    ndist <- ncol(p)  # k or number of classes
156 1
    nfit <- NCOL(x)   # number of documents to be predicted
157

158 1
    distnames <- colnames(p)  # class labels in the training set
159

160 1
    xdim <- dim(x)  # dim of the test set (nfeat x ndoc)
161 1
    stopifnot(NROW(x) == ncat)
162

163 1
    if (length(smooth) == 1) {
164 1
        smooth <- rep(smooth, ndist)
165
    }
166 1
    stopifnot(length(smooth) == ndist)
167 1
    stopifnot(all(smooth >= 0))
168

169 1
    if (!inherits(x, "CsparseMatrix")) {
170 0
        x <- as.matrix(x)
171
    }
172 1
    fitnames <- colnames(x)  # document names to be fitted
173

174 1
    support <- !apply(p == 0, 1, all) # ignore categories with no support
175 1
    p <- p[support, , drop = FALSE]
176 1
    x <- x[support, , drop = FALSE]
177

178
    # fit to each column of x
179 1
    fit <- vector("list", nfit)
180 1
    if (is.matrix(x)) {
181 0
        for (j in seq_len(nfit)) {
182 0
            if (verbose) {
183 0
                cat("Fitting column ", j, "\n", sep = "")
184
            }
185 0
            fit[[j]] <- affinity1(p, x[, j], smooth, verbose)
186
        }
187
    } else {
188 1
        x <- as(x, "dgCMatrix")
189 1
        val <- x@x
190 1
        row_ind <- x@i + 1 # switch from 0- to 1-based indexing
191 1
        col_ptr <- x@p + 1 #
192

193 1
        for (j in seq_len(nfit)) {
194 1
            if (verbose) {
195 0
                cat("Fitting column ", j, "\n", sep = "")
196
            }
197

198 1
            start <- col_ptr[[j]]
199 1
            end <- col_ptr[[j + 1]]
200 1
            len <- end - start
201 1
            ix <- seq.int(start, length.out = len)
202

203 1
            xj <- val[ix]
204 1
            ij <- row_ind[ix]
205 1
            pj <- p[ij, , drop = FALSE]
206

207 1
            fit[[j]] <- affinity1(pj, xj, smooth, verbose)
208
        }
209
    }
210

211
    # simplify results
212 1
    coefficients <- matrix(NA, nfit, ndist, dimnames = list(fitnames, distnames))
213 1
    se <- matrix(NA, nfit, ndist, dimnames = list(fitnames, distnames))
214 1
    cov <- array(NA, c(ndist, ndist, nfit),
215 1
                 dimnames = list(distnames, distnames, fitnames))
216

217 1
    for (j in seq_len(nfit)) {
218 1
        fitj <- fit[[j]]
219 1
        coefficients[j, ] <- fitj$theta
220 1
        se[j, ] <- fitj$theta_se
221 1
        cov[, , j] <- fitj$cov
222
    }
223

224
    # drop dimension if input x was a vector - NOT ANY MORE -kb
225 1
    if (is.null(xdim) && nfit == 1) {
226 0
        coefficients <- coefficients[1, drop = FALSE]
227 0
        se <- se[1, drop = FALSE]
228 0
        cov <- cov[, , 1, drop = FALSE]
229
    }
230

231 1
    list(coefficients = coefficients,
232 1
         se = se,
233 1
         cov = cov,
234 1
         smooth = smooth,
235 1
         support = support)
236
}
237

238

239
#
240
# internal function to maximize the affinity likelihood
241
# author: Patrick Perry
242
#
243
affinity1 <- function(p, x, smooth, verbose = FALSE) {
244 1
    ncat <- nrow(p)
245 1
    ndist <- ncol(p)
246

247 1
    objective <- function(theta, bweight = 0, gradient = TRUE, hessian = TRUE) {
248 1
        q <- drop(p %*% theta)
249

250 1
        loglik <- sum(x * ifelse(x == 0, 0, log(q)))
251 1
        penalty <- -sum(smooth * log(theta))
252 1
        barrier <- -sum(log(theta))
253 1
        value <- (-loglik) + penalty + bweight * barrier
254

255 1
        res <- list(loglik = loglik,
256 1
                    penalty = penalty,
257 1
                    barrier = barrier,
258 1
                    bweight = bweight,
259 1
                    value = value)
260

261 1
        if (gradient || hessian) {
262 1
            score <- drop(t(p) %*% ifelse(x == 0, 0, x / q))
263 1
            grad <- (-score) - ((smooth + bweight) / theta)
264

265 1
            res[["score"]] <- score
266 1
            res[["grad"]] <- grad
267
        }
268

269 1
        if (hessian) {
270 1
            imat_sqrt <- p * ifelse(x == 0, 0, sqrt(x) / q)
271 1
            imat <- t(imat_sqrt) %*% imat_sqrt
272 1
            hess <- imat + diag((smooth + bweight) / theta ^ 2, length(theta))
273

274 1
            res[["imat"]] <- imat
275 1
            res[["hess"]] <- hess
276
        }
277

278 1
        res
279
    }
280

281

282 1
    residuals <- function(theta, nu, bweight = 0) {
283 1
        obj <- objective(theta, bweight, hessian = FALSE)
284 1
        a <- rep(1, length(theta))
285

286 1
        g <- (obj$grad) + a * nu
287 1
        h <- sum(theta) - 1
288 1
        norm <- sqrt(sum(g^2) + h^2)
289

290 1
        list(dual = g, primal = h, norm = norm)
291
    }
292

293

294 1
    newton_step <- function(theta, nu, bweight = 0) {
295 1
        obj <- objective(theta, bweight)
296 1
        a <- rep(1, length(theta))
297

298 1
        g <- (obj$grad) + a * nu
299 1
        h <- sum(theta) - 1
300

301 1
        H <- (obj$hess)
302

303 1
        Hi_a <- solve(H, a)
304 1
        Hi_g <- solve(H, g)
305 1
        s <- drop(-(t(a) %*% Hi_a))
306 1
        w <- drop(t(a) %*% Hi_g - h) / s
307 1
        v <- -(Hi_a * w) - Hi_g
308

309 1
        list(primal = v, dual = w)
310
    }
311

312

313 1
    optimize <- function(theta, nu, bweight = 0, verbose = FALSE) {
314 1
        tol <- 1e-8
315 1
        alpha <- 0.01
316 1
        beta <- 0.5
317

318 1
        resid <- residuals(theta, nu, bweight)
319 1
        iter <- 0
320

321 1
        while (resid$norm > tol || max(abs(resid$primal)) > tol) {
322 1
            if (verbose) {
323 0
                cat("iteration: ", iter, "; residual norm: ", resid$norm, "\n",
324 0
                    sep = "")
325
            }
326 1
            step <- newton_step(theta, nu, bweight)
327

328 1
            tmax <- min(ifelse(step$primal >= 0, Inf, -theta / step$primal))
329 1
            t <- min(1, tmax)
330

331 1
            repeat {
332 1
                theta1 <- theta + t * step$primal
333 1
                nu1 <- nu + t * step$dual
334

335 1
                if (all(theta1 > 0)) {
336 1
                    resid1 <- residuals(theta1, nu1, bweight)
337 1
                    if (resid1$norm <= (1 - alpha * t) * resid$norm) {
338 1
                        break
339
                    }
340
                }
341 1
                t <- beta * t
342
            }
343 1
            theta <- theta1
344 1
            nu <- nu1
345 1
            resid <- resid1
346 1
            iter <- iter + 1
347
        }
348

349 1
        obj <- objective(theta, bweight)
350 1
        list(theta = theta,
351 1
             nu = nu,
352 1
             resid = resid,
353 1
             objective = obj,
354 1
             iter = iter)
355
    }
356

357

358 1
    shrink_barrier <- function(theta, nu, bweight = 0, verbose = FALSE) {
359 0
        shrink <- 0.1
360

361 0
        if (verbose) {
362 0
            cat("bweight: ", bweight, "\n", sep = "")
363
        }
364 0
        opt <- optimize(theta, nu, bweight, verbose = verbose)
365

366 0
        while (bweight >= 1e-8) {
367 0
            theta <- opt$theta
368 0
            nu <- opt$nu
369 0
            bweight <- shrink * bweight
370 0
            if (verbose) {
371 0
                cat("bweight: ", bweight, "\n", sep = "")
372
            }
373 0
            opt <- optimize(theta, nu, bweight, verbose = verbose)
374
        }
375 0
        opt[["bweight"]] <- bweight
376 0
        opt
377
    }
378

379

380 1
    theta <- rep(1 / ndist, ndist)
381 1
    nu <- 0
382

383 1
    if (any(smooth == 0)) {
384 0
        bweight <- 100
385 0
        res <- shrink_barrier(theta, nu, bweight, verbose = verbose)
386
    } else {
387 1
        res <- optimize(theta, nu, 0, verbose = verbose)
388
    }
389

390 1
    obj <- res$objective
391

392
    # The covariance matrix C solves the block system
393
    #
394
    # [ H   a ] [ C   w ] = [ I 0 ]
395
    # [ a^T 0 ] [ w^T s ] = [ 0 1 ]
396
    #
397
    # where a = (1, 1, ..., 1)^T.
398
    #
399 1
    H <- obj$hess
400 1
    a <- rep(1, ndist)
401

402
    # First solve for w^t:
403
    #   w^T = (a^T H^{-1} a)^{-1} a^T H^{-1}
404
    #
405 1
    Hi <- solve(H)
406 1
    Hi_a <- solve(H, a)
407 1
    a_Hi_a <- sum(a * Hi_a)
408

409
    # Then substitute for C:
410
    #   C = H^{-1} - H^{-1} a (a^T H^{-1} a)^{-1} a^T H^{-1}
411
    #
412 1
    cov <- Hi - (1 / a_Hi_a) * (Hi_a %*% t(Hi_a))
413

414 1
    res[["cov"]] <- cov
415 1
    res[["theta_se"]] <- sqrt(pmax(diag(cov), 0))
416 1
    res
417
}
418

419
# supporting methods for textmodel_affinity ------------
420

421
#' @method print textmodel_affinity
422
#' @export
423
print.textmodel_affinity <- function(x, ...) {
424 1
    cat("Call:\n")
425 1
    print(x$call)
426

427 1
    ref <- table(x$y)
428 1
    namez <- names(ref)
429 1
    namez[2:length(namez)] <- paste(",", namez[2:length(namez)])
430

431 1
    cat("\n",
432 1
        "Training documents per class:",
433 1
        paste0(interleave(paste0(namez, ": "), as.integer(ref)), collapse = ""), "; ",
434 1
        "total training features: ", nrow(x$p), "\n",
435 1
        sep = "")
436
}
437

438
## ============== Prediction Methods =======================================
439

440
#' Prediction for a fitted affinity textmodel
441
#'
442
#' @description
443
#' Estimate \eqn{\theta_i} for each document, from a fitted
444
#' [textmodel_affinity] object.
445
#'
446
#' Other methods below provide standard ways to extract or compute quantities
447
#' from predicted [textmodel_affinity] objects.
448
#' @param object a fitted affinity textmodel
449
#' @param level probability level for confidence interval width
450
#' @param newdata dfm on which prediction should be made
451
#' @param ... unused
452
#' @return `predict()` returns a list of predicted affinity textmodel
453
#'   quantities.
454
#' @importFrom methods new
455
#' @importFrom stats predict
456
#' @method predict textmodel_affinity
457
#' @keywords textmodel internal
458
#' @seealso [influence.predict.textmodel_affinity()] for methods of
459
#'   computing the influence of particular features from a predicted
460
#'   [textmodel_affinity] model.
461
#' @export
462
predict.textmodel_affinity <- function(object, newdata = NULL,
463
                                       level = 0.95, ...) {
464 1
    if (length(list(...)) > 0)
465 0
        stop("Arguments:", names(list(...)), "not supported.\n")
466

467 1
    if (!is.null(newdata)) {
468 0
        data <- newdata
469 0
        train <- rep(FALSE, nrow(data))
470
    } else {
471 1
        data <- object$x
472 1
        newdata <- data
473 1
        train <- !is.na(object$y)
474
    }
475

476 1
    predicted <- affinity(object$p, t(newdata), smooth = object$smooth)
477

478 1
    result <- list(
479 1
        coefficients = predicted$coefficients,
480 1
        se = predicted$se,
481 1
        cov = predicted$cov,
482 1
        smooth = object$smooth,
483 1
        newdata = newdata,
484 1
        train = train,
485 1
        level = level,
486 1
        p = object$p,
487 1
        support = object$support)
488 1
    class(result) <- c("predict.textmodel_affinity", "list")
489 1
    return(result)
490
}
491

492
#' @method print predict.textmodel_affinity
493
#' @export
494
print.predict.textmodel_affinity <- function(x, ...) {
495 0
    print(unclass(x))
496
}
497

498
#' @rdname predict.textmodel_affinity
499
#' @method coef predict.textmodel_affinity
500
#' @return `coef()` returns a document \eqn{\times} class matrix of class
501
#'   affinities for each document.
502
#' @export
503
coef.predict.textmodel_affinity <- function(object, ...) {
504 1
    object$coefficients
505
}
506

507
#' @method coefficients predict.textmodel_affinity
508
#' @export
509
coefficients.predict.textmodel_affinity <- function(object, ...) {
510 0
    UseMethod("coef")
511
}
512

513
#' @rdname predict.textmodel_affinity
514
#' @return
515
#' `residuals()` returns a document-by-feature matrix of residuals.
516
#' `resid()` is an alias.
517
#' @method residuals predict.textmodel_affinity
518
#' @param type see [residuals.lm]
519
#' @importFrom stats residuals resid
520
#' @export
521
residuals.predict.textmodel_affinity <- function(object, type = c("response", "pearson"), ...) {
522

523 0
    type <- match.arg(type)
524 0
    expected <- coef(object) %*% t(object$p) * rowSums(object$newdata)
525 0
    r <- object$newdata - expected
526 0
    if (type == "response") {
527 0
        res <- r
528 0
    } else if (type == "pearson") {
529 0
        res <- r / sqrt(expected)
530
    }
531 0
    res[, !object$support] <- NA
532 0
    as.matrix(res)
533
}
534

535
#' @export
536
#' @method resid predict.textmodel_affinity
537
resid.predict.textmodel_affinity <- function(object, ...) {
538 0
    UseMethod("residuals", ...)
539
}
540

541
#' @rdname predict.textmodel_affinity
542
#' @method rstandard predict.textmodel_affinity
543
#' @return `rstandard()` is a shortcut to return the pearson residuals.
544
#' @importFrom stats rstandard sd
545
#' @export
546
rstandard.predict.textmodel_affinity <- function(model, ...) {
547 0
    residuals(model, type = "pearson")
548
}
549

550

551
# ============== Influence methods =============
552

553
#' Compute feature influence from a predicted textmodel_affinity object
554
#'
555
#' Computes the influence of features on scaled [textmodel_affinity]
556
#' applications.
557
#' @param model a predicted
558
#'   [textmodel_affinity][predict.textmodel_affinity] object
559
#' @param subset whether to use all data or a subset (for instance, exclude the
560
#'   training set)
561
#' @param ... unused
562
#' @seealso [influence.lm()]
563
#' @keywords textmodel internal
564
#' @importFrom stats influence
565
#' @method influence predict.textmodel_affinity
566
#' @examples
567
#' tmot <- textmodel_affinity(data_dfm_lbgexample, y = c("L", NA, NA, NA, "R", NA))
568
#' pred <- predict(tmot)
569
#' influence(pred)
570
#' @export
571
influence.predict.textmodel_affinity <- function(model, subset = !train, ...) {
572
    # subset/training set
573 1
    train <- model$train
574

575
    # reference distributions
576 1
    p <- model$p
577 1
    levels <- colnames(p)
578 1
    support <- model$support
579

580
    # class affinity estimates
581 1
    theta <- model$coefficients[subset, , drop = FALSE]
582 1
    cov <- model$cov[, , subset, drop = FALSE]
583

584
    # data
585 1
    x <- model$newdata[subset, ]
586 1
    x[, !support] <- 0
587 1
    x <- as(t(x), "dgCMatrix")
588 1
    nword <- nrow(x)
589 1
    words <- rownames(x)
590 1
    ntext <- ncol(x)
591 1
    texts <- colnames(x)
592

593 1
    val <- x@x
594 1
    row_ind <- x@i + 1 # switch from 0- to 1-based indexing
595 1
    col_ptr <- x@p + 1 #
596

597 1
    infl_norm <- numeric(length(val))
598 1
    infl_mode <- integer(length(val))
599 1
    rate <- numeric(length(val))
600

601 1
    for (j in seq_len(ntext)) {
602 1
        start <- col_ptr[[j]]
603 1
        end <- col_ptr[[j + 1]]
604 1
        len <- end - start
605 1
        ix <- seq.int(start, length.out = len)
606

607 1
        xj <- val[ix]
608 1
        ij <- row_ind[ix]
609

610 1
        pj <- p[ij,,drop = FALSE]
611 1
        mu <- as.vector(pj %*% theta[j, ])
612 1
        q <- pj / ifelse(mu == 0, 1, mu)
613

614 1
        infl_dir <- as.matrix(q %*% cov[, , j])
615 1
        h2 <- rowSums(q * infl_dir)
616

617
        # crude approximation for Hessian:
618
        # infl_norm[ix] <- 0.5 * xj * rowSums(abs(infl_dir))
619

620
        # more accurate approximation:
621 1
        infl_norm[ix] <- 0.5 * abs(xj / (1 - xj * h2)) * rowSums(abs(infl_dir))
622

623 1
        infl_mode[ix] <- apply(infl_dir, 1, which.max)
624 1
        rate[ix] <- xj / sum(xj)
625
    }
626

627
    # Note: why not just use Matrix::t()? KW
628 1
    transpose <- function(values, as_matrix = FALSE) {
629 1
        mat <- t(Matrix::sparseMatrix(i = x@i, p = x@p, x = values,
630 1
                                      dims = c(nword, ntext),
631 1
                                      dimnames = list(words, texts),
632 1
                                      index1 = FALSE))
633 1
        if (as_matrix) {
634 1
            ans <- mat
635
        } else {
636 1
            ans <- mat@x
637
        }
638 1
        ans
639
    }
640

641 1
    norm <- transpose(infl_norm, as_matrix = TRUE)
642 1
    rate <- transpose(rate)
643 1
    count <- transpose(val)
644 1
    mode <- transpose(infl_mode)
645

646 1
    res <- list(norm = norm, count = count, rate = rate,
647 1
                mode = mode, levels = levels, subset = subset,
648 1
                support = support)
649 1
    class(res) <- "influence.predict.textmodel_affinity"
650 1
    res
651
}
652

653
#' @title Internal methods for textmodel_affinity
654
#' @description Internal print and summary methods for derivative
655
#'   [textmodel_affinity] objects.
656
#' @name textmodel_affinity-internal
657
#' @keywords textmodel internal
658
#' @method print influence.predict.textmodel_affinity
659
#' @param n how many coefficients to print before truncating
660
#' @export
661
print.influence.predict.textmodel_affinity <- function(x, n = 30, ...) {
662 0
    ans <- summary(x, ...)
663 0
    print(ans, n)
664
}
665

666
#' @rdname textmodel_affinity-internal
667
#' @method summary influence.predict.textmodel_affinity
668
#' @importFrom stats median
669
#' @export
670
summary.influence.predict.textmodel_affinity <- function(object, ...) {
671 1
    norm <- object$norm
672 1
    ntext <- nrow(norm)
673 1
    nword <- ncol(norm)
674 1
    words <- colnames(norm)
675

676 1
    val <- norm@x
677 1
    row_ind <- norm@i + 1
678 1
    col_ptr <- norm@p + 1
679 1
    mode <- object$mode
680 1
    count <- object$count
681 1
    rate <- object$rate
682

683 1
    count_val <- numeric(nword)
684 1
    mean_val <- numeric(nword)
685 1
    med_val <- numeric(nword)
686 1
    sd_val <- numeric(nword)
687 1
    max_val <- numeric(nword)
688

689 1
    max_count <- numeric(nword)
690 1
    max_rate <- numeric(nword)
691 1
    med_rate <- numeric(nword)
692 1
    max_dir <- numeric(nword)
693

694 1
    for (j in seq_len(nword)) {
695 1
        start <- col_ptr[[j]]
696 1
        end <- col_ptr[[j + 1]]
697 1
        len <- end - start
698

699 1
        if (len > 0) {
700 1
            ix <- seq.int(start, length.out = len)
701 1
            xj <- val[ix]
702 1
            ij <- row_ind[ix]
703

704 1
            m <- which.max(xj)
705

706 1
            count_val[j] <- len
707 1
            mean_val[j] <- mean(xj)
708 1
            med_val[j] <- median(xj)
709 1
            sd_val[j] <- ifelse(len > 1, sd(xj), 0)
710 1
            max_val[j] <- xj[m] # == val[ix[m]]
711

712 1
            max_count[j] <- count[ix[m]]
713 1
            max_rate[j] <- rate[ix[m]]
714 1
            med_rate[j] <- median(rate[ix])
715 1
            max_dir[j] <- mode[ix[m]]
716
        } else {
717 1
            count_val[j] <- 0
718 1
            mean_val[j] <- 0
719 1
            med_val[j] <- 0
720 1
            sd_val[j] <- 0
721 1
            max_val[j] <- 0
722

723 1
            max_count[j] <- 0
724 1
            max_rate[j] <- 0
725 1
            med_rate[j] <- 0
726 1
            max_dir[j] <- NA
727
        }
728
    }
729

730 1
    labels <- object$levels
731 1
    levels <- seq_along(labels)
732 1
    max_dir <- factor(max_dir, levels, labels)
733

734 1
    result <- list(word = words, count = count_val,
735 1
                   mean = mean_val, median = med_val,
736 1
                   sd = sd_val, max = max_val,
737 1
                   direction = max_dir,
738 1
                   rate = med_rate,
739 1
                   support = object$support)
740 1
    class(result) <- "summary.influence.predict.textmodel_affinity"
741 1
    result
742
}
743

744
#' @rdname textmodel_affinity-internal
745
#' @method print summary.influence.predict.textmodel_affinity
746
#' @param n how many coefficients to print before truncating
747
#' @export
748
print.summary.influence.predict.textmodel_affinity <- function(x, n = 30, ...) {
749 0
    ix <- sort(x$median, decreasing = TRUE, index.return = TRUE)$ix
750 0
    influence <- x$median[ix]
751

752 0
    d <- data.frame(word = x$word,
753 0
                    count = x$count,
754 0
                    median = x$median,
755 0
                    max = x$max,
756 0
                    direction = x$direction)
757 0
    d <- d[ix,]
758 0
    rownames(d) <- d$word
759 0
    d$word <- NULL
760

761 0
    if (!is.null(n) && !is.na(n)) {
762 0
        n <- min(n, nrow(d))
763 0
        d <- d[seq_len(n),,drop=FALSE]
764

765 0
        cat("Top ", n, " influential words:\n\n", sep = "")
766
    }
767 0
    print(d)
768

769 0
    invisible(x)
770
}
771

772

773
# ========= Internal functions =========
774

775
# compute chi^2 goodness of fit
776
gof_chi2 <- function(x) {
777 0
    UseMethod("gof_chi2")
778
}
779
gof_chi2.predict.textmodel_affinity <- function(x) {
780 0
    rowSums(rstandard(x)[,x$support,drop=FALSE]^2)
781
}
782

783
# function to interleave two vector objects
784
# Example:
785
# interleave(letters[1:3], 1:3)
786
# ## [1] "a" "1" "b" "2" "c" "3"
787
interleave <- function(v1, v2) {
788 1
    ord1 <- 2 * (seq_along(v1)) - 1
789 1
    ord2 <- 2 * (seq_along(v2))
790 1
    c(v1, v2)[order(c(ord1, ord2))]
791
}
792

793

794

Read our documentation on viewing source code .

Loading