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
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
### iid2.R ---
#----------------------------------------------------------------------
## author: Brice Ozenne
## created: okt 12 2017 (13:16)
## Version:
## last-updated: jan 18 2022 (09:38)
## By: Brice Ozenne
## Update #: 625
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
## * Documentation - iid2
#' @title Influence Function With Small Sample Correction.
#' @description Extract the influence function from a latent variable model.
#' It is similar to \code{lava::iid} but with small sample correction.
#' @name iid2
#'
#' @param object,x a \code{lvmfit} or \code{lvmfit2} object (i.e. output of \code{lava::estimate} or \code{lavaSearch2::estimate2}).
#' @param cluster [integer vector] the grouping variable relative to which the observations are iid.
#' @param robust [logical] if \code{FALSE}, the influence function is rescaled such its the squared sum equals the model-based standard error (instead of the robust standard error).
#' Do not match the model-based correlation though.
#' @param as.lava [logical] if \code{TRUE}, uses the same names as when using \code{stats::coef}.
#' @param ssc [character] method used to correct the small sample bias of the variance coefficients (\code{"none"}, \code{"residual"}, \code{"cox"}). Only relevant when using a \code{lvmfit} object.
#' @param ... additional argument passed to \code{estimate2} when using a \code{lvmfit} object.
#'
#' @details When argument object is a \code{lvmfit} object, the method first calls \code{estimate2} and then extract the variance-covariance matrix.
#'
#' @seealso \code{\link{estimate2}} to obtain \code{lvmfit2} objects.
#'
#' @return A matrix containing the 1st order influence function relative to each sample (in rows)
#' and each model coefficient (in columns).
#'
#' @examples
#' #### simulate data ####
#' n <- 5e1
#' p <- 3
#' X.name <- paste0("X",1:p)
#' link.lvm <- paste0("Y~",X.name)
#' formula.lvm <- as.formula(paste0("Y~",paste0(X.name,collapse="+")))
#'
#' m <- lvm(formula.lvm)
#' distribution(m,~Id) <- Sequence.lvm(0)
#' set.seed(10)
#' d <- sim(m,n)
#'
#' #### latent variable model ####
#' e.lvm <- estimate(lvm(formula.lvm),data=d)
#' iid.tempo <- iid2(e.lvm)
#'
#'
#' @concept extractor
#' @keywords smallSampleCorrection
#' @export
`iid2` <-
function(object, ...) UseMethod("iid2")
## * iid2.lvmfit
#' @rdname iid2
#' @export
iid2.lvmfit <- function(object, robust = TRUE, cluster = NULL, as.lava = TRUE, ssc = lava.options()$ssc,...){
return(lava::iid(estimate2(object, ssc = ssc, ...), robust = robust, cluster = cluster, as.lava = as.lava))
}
## * iid2.lvmfit2
#' @rdname iid2
#' @export
iid2.lvmfit2 <- function(object, robust = TRUE, cluster = NULL, as.lava = TRUE, ...){
dots <- list(...)
if(length(dots)>0){
warning("Argument(s) \'",paste(names(dots),collapse="\' \'"),"\' not used by ",match.call()[1],". \n")
}
## ** compute iid
object.score <- score2(object, indiv = TRUE, cluster = cluster, as.lava = FALSE)
object.vcov <- vcov2(object, as.lava = FALSE)
object.iid <- object.score %*% object.vcov
if(robust == FALSE){
object.chol.vcov <- matrixPower(object.vcov, power = 1/2, symmetric = TRUE)
object.chol.rvcov <- matrixPower(crossprod(object.iid[rowSums(!is.na(object.iid))>0,,drop=FALSE]), power = -1/2, symmetric = TRUE)
object.iid <- object.iid %*% object.chol.rvcov %*% object.chol.vcov
}
## restaure names
colnames(object.iid) <- colnames(object.score)
if(!is.null(cluster)){
attr(object.iid,"cluster") <- attr(object.score,"cluster")
}
## ** export
if(as.lava){
object.iid <- object.iid[,names(object$sCorrect$skeleton$originalLink2param),drop=FALSE]
colnames(object.iid) <- as.character(object$sCorrect$skeleton$originalLink2param)
}
return(object.iid)
}
## * iid.lvmfit2
#' @rdname iid2
#' @export
iid.lvmfit2 <- function(x, robust = TRUE, cluster = NULL, as.lava = TRUE, ...){
return(iid2(x, robust = robust, cluster = cluster, as.lava = as.lava, ...)) ## necessary as first argument of iid must be x
}
## * Documentation - iid2plot
#' @title Display the i.i.d. Decomposition
#' @description Extract the i.i.d. decomposition and display it along with the corresponding coefficient.
#'
#' @param object a \code{lvmfit} or \code{lvmfit2} object (i.e. output of \code{lava::estimate} or \code{lavaSearch2::estimate2}).
#' @param param [character] name of one of the model parameters.
#'
#' @export
iid2plot <- function(object, param){
all.param <- coef2(object)
if(length(param) != 1){
stop("Incorrect argument \'param\'. \n",
"Should have length 1. \n")
}
if(param %in% names(all.param) == FALSE){
stop("Incorrect argument \'param\'. \n",
"Should be one of the model parameters. \n")
}
object.iid <- iid2(object)[,param]
h <- graphics::hist(object.iid, main = paste0(param,"=",all.param[param]), xlab = "iid")
return(invisible(h))
}
##----------------------------------------------------------------------
### iid2.R ends here
|
