functions.R 10.7 KB
Newer Older
Rauschenberger's avatar
Rauschenberger committed
1

Rauschenberger's avatar
Rauschenberger committed
2
#' @export
Rauschenberger's avatar
Rauschenberger committed
3
#' @aliases colasso-package
Rauschenberger's avatar
Rauschenberger committed
4
#' @title
Rauschenberger's avatar
Rauschenberger committed
5
#' colasso
Rauschenberger's avatar
Rauschenberger committed
6
7
#' 
#' @description
Armin Rauschenberger's avatar
Armin Rauschenberger committed
8
#' Implements penalised regression with response moderation.
Rauschenberger's avatar
Rauschenberger committed
9
#'  
Rauschenberger's avatar
Rauschenberger committed
10
11
12
#' @param y
#' response\strong{:}
#' vector of length \eqn{n}
Rauschenberger's avatar
Rauschenberger committed
13
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
14
15
16
17
18
#' @param Y
#' response\strong{:}
#' matrix with \eqn{n} rows and \eqn{p} columns,
#' or vector of length \eqn{n} (see details)
#' 
Rauschenberger's avatar
Rauschenberger committed
19
20
21
#' @param X
#' covariates\strong{:}
#' matrix with \eqn{n} rows (samples) and \eqn{p} columns (variables)
Rauschenberger's avatar
Rauschenberger committed
22
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
23
24
25
26
27
28
#' @param alpha
#' elastic net parameter\strong{:}
#' numeric between \eqn{0} and \eqn{1};
#' \eqn{alpha=1} for lasso,
#' \eqn{alpha=0} for ridge
#' 
Rauschenberger's avatar
Rauschenberger committed
29
#' @param nfolds
Rauschenberger's avatar
Rauschenberger committed
30
#' number of folds
Rauschenberger's avatar
Rauschenberger committed
31
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
32
33
34
35
36
#' @param family
#' see glmnet
#' 
#' @param type.measure
#' see glmnet
Rauschenberger's avatar
Rauschenberger committed
37
#' 
Rauschenberger's avatar
Rauschenberger committed
38
#' @examples
Armin Rauschenberger's avatar
Armin Rauschenberger committed
39
#' n <- 100; p <- 20; q <- 10
Rauschenberger's avatar
Rauschenberger committed
40
#' X <- matrix(rnorm(n*p),nrow=n,ncol=p)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
41
#' Y <- matrix(rnorm(n*q),nrow=n,ncol=q)
Rauschenberger's avatar
Rauschenberger committed
42
43
#' #y <- rbinom(n=n,size=1,prob=0.2)
#' y <- rnorm(n=n)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
44
#' test <- colasso(y=y,Y=Y,X=X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
45
46
47
48
49
#' 
colasso <- function(y,Y,X,alpha=1,nfolds=10,family="gaussian",type.measure="deviance"){
  
  # properties
  n <- nrow(X); p <- ncol(X)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
50
51
52
  if(!family %in% c("gaussian","poisson","binomial")){
    stop("Family not implemented.")
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
53
  if(length(y)!=n){stop("sample size")}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
54
55
  #foldid <- sample(x=rep(x=seq_len(nfolds),length.out=n))
  foldid <- palasso:::.folds(y=y,nfolds=nfolds)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
56
57
58
  
  # weights
  pi <- seq(from=0,to=1,by=0.2) # adapt this
Rauschenberger's avatar
Rauschenberger committed
59

Armin Rauschenberger's avatar
Armin Rauschenberger committed
60
61
  # model fitting
  fit <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
62
  ym <- colasso::colasso_moderate(Y=Y) # trial
Armin Rauschenberger's avatar
Armin Rauschenberger committed
63
64
65
66
  for(i in seq_along(pi)){
    weights <- rep(c(1-pi[[i]],pi[[i]]),each=n)
    fit[[i]] <- glmnet::glmnet(y=c(y,ym),x=rbind(X,X),weights=weights,family=family,alpha=alpha)
  }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
67
  names(fit) <- paste0("pi",pi)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
68
69
70
  
  # inner cross-validation
  pred <- lapply(pi,function(x) matrix(data=NA,nrow=length(y),ncol=100))
Armin Rauschenberger's avatar
Armin Rauschenberger committed
71
  for(k in unique(foldid)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
72
    y0 <- y[foldid!=k]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
73
    y1 <- y[foldid==k]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
74
    Y0 <- Y[foldid!=k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
75
    Y1 <- Y[foldid==k,,drop=FALSE]
Armin Rauschenberger's avatar
Armin Rauschenberger committed
76
77
    X0 <- X[foldid!=k,,drop=FALSE]
    X1 <- X[foldid==k,,drop=FALSE]
Rauschenberger's avatar
Rauschenberger committed
78
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
79
    y0m <- colasso_moderate(Y=Y0)
Rauschenberger's avatar
Rauschenberger committed
80
    for(i in seq_along(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
81
82
83
84
      weights <- rep(c(1-pi[[i]],pi[[i]]),each=sum(foldid!=k)) # trial
      glmnet <- glmnet::glmnet(y=c(y0,y0m),x=rbind(X0,X0),weights=weights,family=family,alpha=alpha)
      temp <- stats::predict(object=glmnet,newx=X1,type="response",s=fit[[i]]$lambda)
      pred[[i]][foldid==k,seq_len(ncol(temp))] <- temp
Rauschenberger's avatar
Rauschenberger committed
85
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
86
87
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
88
  # loss sequence 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
89
  for(i in seq_along(pi)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
90
91
92
93
94
95
96
97
98
    # WATCH OUT: adapt to all loss fuctions
    #fit[[i]]$cvm <- apply(X=pred[[i]],MARGIN=2,FUN=function(x) mean((x-y)^2))
    fit[[i]]$cvm <- palasso:::.loss(y=y,fit=pred[[i]],family=family,type.measure=type.measure,foldid=foldid)[[1]]
    # WATCH OUT: minimise or maximise
    if(type.measure=="AUC"){
      fit[[i]]$lambda.min <- fit[[i]]$lambda[which.max(fit[[i]]$cvm)]
    } else {
      fit[[i]]$lambda.min <- fit[[i]]$lambda[which.min(fit[[i]]$cvm)]
    }
Armin Rauschenberger's avatar
Armin Rauschenberger committed
99
100
  }
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
101
102
103
104
105
106
  # loss sequence
  #cvm <- palasso:::.loss(y=y,fit=pred,family=family,type.measure=type.measure,foldid=foldid)
  
  # optimisation
  #model <- .extract(fit=fit.full,lambda=lambda,cvm=cvm,type.measure=args$type.measure)
  
Armin Rauschenberger's avatar
Armin Rauschenberger committed
107
108
109
110
111
112
113
114
115
116
  # selection
  cvm <- sapply(fit,function(x) x$cvm[which(x$lambda==x$lambda.min)])
  if(type.measure=="AUC"){
    sel <- which.max(cvm)
  } else {
    sel <- which.min(cvm)
  }
  fit[[length(pi)+1]] <- fit[[sel]]
  
  #graphics::plot(cvm); graphics::abline(v=sel,lty=2)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
117
  names(fit) <- c("glmnet",paste0("pi",pi[-1]),"conet")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
118
  return(fit)
Rauschenberger's avatar
Rauschenberger committed
119
120
121
}


Rauschenberger's avatar
Rauschenberger committed
122
123
#' @export
#' @title
Rauschenberger's avatar
Rauschenberger committed
124
#' moderated response
Armin Rauschenberger's avatar
Armin Rauschenberger committed
125
#'
Rauschenberger's avatar
Rauschenberger committed
126
#' @description
Rauschenberger's avatar
Rauschenberger committed
127
#' This function ...
Armin Rauschenberger's avatar
Armin Rauschenberger committed
128
#'
Rauschenberger's avatar
Rauschenberger committed
129
#' @inheritParams colasso
Armin Rauschenberger's avatar
Armin Rauschenberger committed
130
131
132
#' 
#' @param ...
#' further arguments (currently not implemented)
Rauschenberger's avatar
Rauschenberger committed
133
#' vector with entries between \eqn{0} and \eqn{1} (rename argument)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
134
135
136
#'
#'
#'
Rauschenberger's avatar
Rauschenberger committed
137
138
#' @examples
#' NA
Armin Rauschenberger's avatar
Armin Rauschenberger committed
139
colasso_moderate <- function(Y,...){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
140
141
  # (most basic version possible)
  y <- rowMeans(Y)
Armin Rauschenberger's avatar
Armin Rauschenberger committed
142
143
  y <- apply(Y,1,stats::median)
  if(all(y %in% c(0,0.5,1))){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
144
145
146
147
    y[y==0.5] <- 1
    warning("Invalid unless binomial family.")
  }
  return(y)
Rauschenberger's avatar
Rauschenberger committed
148
149
}

Rauschenberger's avatar
Rauschenberger committed
150
151
#' @export
#' @title
Rauschenberger's avatar
Rauschenberger committed
152
#' simulate data
Rauschenberger's avatar
Rauschenberger committed
153
154
155
#' 
#' @description
#' This function ...
Rauschenberger's avatar
Rauschenberger committed
156
157
158
#'  
#' @param n
#' sample size
Rauschenberger's avatar
Rauschenberger committed
159
#' 
Rauschenberger's avatar
Rauschenberger committed
160
161
#' @param p
#' number of covariates
Rauschenberger's avatar
Rauschenberger committed
162
#' 
Rauschenberger's avatar
Rauschenberger committed
163
164
#' @param cor
#' correlation structure
Rauschenberger's avatar
Rauschenberger committed
165
#' 
Rauschenberger's avatar
Rauschenberger committed
166
167
#' @param plot
#' logical
Rauschenberger's avatar
Rauschenberger committed
168
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
169
170
171
#' @param family
#' character
#' 
Rauschenberger's avatar
Rauschenberger committed
172
#' @examples
Rauschenberger's avatar
Rauschenberger committed
173
174
#' # CONTINUE HERE
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
175
colasso_simulate <- function(n=100,p=500,cor="constant",family="gaussian",plot=TRUE){
Rauschenberger's avatar
Rauschenberger committed
176
177
178
179
180
181
182
    # correlation matrix
    if(cor=="none"){
        Sigma <- matrix(data=0,nrow=p,ncol=p)
    } else if(cor=="constant"){
        Sigma <- matrix(data=0.05,nrow=p,ncol=p)
    } else if(cor=="autoregressive"){
        # adjust 0.9 to p, such that mean(Sigma)=0.05
Armin Rauschenberger's avatar
Armin Rauschenberger committed
183
        # sum(2*(p-seq_len(p)+1)*0.9^seq_len(p))/(p*p)
Rauschenberger's avatar
Rauschenberger committed
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
        Sigma <- 0.9^abs(col(diag(p))-row(diag(p)))
    } else if(cor=="unstructured"){
        Sigma <- matrix(data=stats::rbeta(n=p,shape1=0.05,shape2=1),nrow=p,ncol=p)
    }
    diag(Sigma) <- 1
    
    X <- MASS::mvrnorm(n=n,mu=rep(0,p),Sigma=Sigma)
    stats::median(abs(as.numeric(cor(X))))

    # non-sparse effects
    #beta <- stats::rnorm(n=p,mean=0,sd=1)
    #mu <- X %*% beta
    #y <- stats::rnorm(n=n,mean=mu)

    # sparse effects
Armin Rauschenberger's avatar
Armin Rauschenberger committed
199
200
    beta <- rep(x=1,times=p) 
    beta[stats::rbinom(n=p,size=1,prob=0.95)==1] <- 0
Rauschenberger's avatar
Rauschenberger committed
201
    mu <- X %*% beta
Armin Rauschenberger's avatar
Armin Rauschenberger committed
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
    
    if(family=="gaussian"){
      Y <- replicate(n=10,expr=stats::rnorm(n=n,mean=mu,sd=10))
    } else if(family=="binomial"){
      prob <- exp(mu)/(1+exp(mu))
      Y <- replicate(n=10,expr=stats::rbinom(n=n,size=1,prob=prob))
    } else if(family=="poisson"){
      lambda <- exp(mu)
      Y <- replicate(n=10,expr=stats::rpois(n=n,lambda=lambda))
    } else if(family=="cox"){
      warning("Cox regression not implemented!")
    }
    
    #y <- stats::rnorm(n=n,mean=mu,sd=10)
    y <- Y[,1]
Rauschenberger's avatar
Rauschenberger committed
217
218
219

    # predictivity -----------------------------------------------------------------
    if(plot){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
220
        graphics::par(mar=c(3,3,1,1))
Rauschenberger's avatar
Rauschenberger committed
221
222
223
        test <- glmnet::cv.glmnet(x=X,y=y)
        graphics::plot(x=log(test$lambda),y=test$cvm)
        graphics::abline(h=test$cvm[test$lambda==max(test$lambda)],lty=2)
Rauschenberger's avatar
Rauschenberger committed
224
    }
Rauschenberger's avatar
Rauschenberger committed
225
        
Armin Rauschenberger's avatar
Armin Rauschenberger committed
226
    return(list(y=y,Y=Y,X=X))
Rauschenberger's avatar
Rauschenberger committed
227
228
}

Rauschenberger's avatar
Rauschenberger committed
229
230
#' @export
#' @title
Rauschenberger's avatar
Rauschenberger committed
231
#' External cross-validation
Rauschenberger's avatar
Rauschenberger committed
232
233
#' 
#' @description
Rauschenberger's avatar
Rauschenberger committed
234
235
#' This function ...
#'  
Rauschenberger's avatar
Rauschenberger committed
236
#' @param y
Armin Rauschenberger's avatar
Armin Rauschenberger committed
237
238
239
240
#' response vector
#' 
#' @param Y
#' response moderation matrix
Rauschenberger's avatar
Rauschenberger committed
241
#' 
Rauschenberger's avatar
Rauschenberger committed
242
243
#' @param X
#' covariates
Rauschenberger's avatar
Rauschenberger committed
244
#' 
Rauschenberger's avatar
Rauschenberger committed
245
246
#' @param plot
#' logical
Rauschenberger's avatar
Rauschenberger committed
247
#' 
Rauschenberger's avatar
Rauschenberger committed
248
249
250
#' @param nfolds.int
#' internal folds
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
251
252
253
#' @param nfolds.ext
#' external folds
#' 
Armin Rauschenberger's avatar
Armin Rauschenberger committed
254
255
#' @inheritParams colasso
#' 
Rauschenberger's avatar
Rauschenberger committed
256
#' @examples
Rauschenberger's avatar
Rauschenberger committed
257
258
#' NA
#'
Armin Rauschenberger's avatar
Armin Rauschenberger committed
259
colasso_compare <- function(y,Y,X,plot=TRUE,nfolds.ext=5,nfolds.int=10,family="gaussian",type.measure="deviance"){
Rauschenberger's avatar
Rauschenberger committed
260
    
Armin Rauschenberger's avatar
Armin Rauschenberger committed
261
    fold <- sample(x=rep(x=seq_len(nfolds.ext),length.out=length(y)))
Rauschenberger's avatar
Rauschenberger committed
262
    pred <- matrix(data=NA,nrow=length(y),ncol=8)
Rauschenberger's avatar
Rauschenberger committed
263
    select <- list()
Armin Rauschenberger's avatar
Armin Rauschenberger committed
264
    for(i in unique(fold)){
Rauschenberger's avatar
Rauschenberger committed
265
        cat("i =",i,"\n")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
266
        fit <- colasso(y=y[fold!=i],Y=Y[fold!=i,],X=X[fold!=i,],alpha=1,nfolds=nfolds.int,type.measure=type.measure)
Rauschenberger's avatar
Rauschenberger committed
267
        for(j in seq_along(fit)){
Armin Rauschenberger's avatar
Armin Rauschenberger committed
268
          # REPLACE glmnet::predict.glmnet by stats::predict !!!
Rauschenberger's avatar
Rauschenberger committed
269
270
271
272
            pred[fold==i,j] <- glmnet::predict.glmnet(object=fit[[j]],
                                                  newx=X[fold==i,],
                                                  s=fit[[j]]$lambda.min,
                                                  type="response")
Rauschenberger's avatar
Rauschenberger committed
273
        }
Rauschenberger's avatar
Rauschenberger committed
274
        select[[i]] <- lapply(fit,function(x) which(x$beta[,x$lambda==x$lambda.min]!=0))
Rauschenberger's avatar
Rauschenberger committed
275
        pred[fold==i,8] <- mean(y[fold!=i]) # intercept-only model
Rauschenberger's avatar
Rauschenberger committed
276
    }
Rauschenberger's avatar
Rauschenberger committed
277
    colnames(pred) <- c(names(fit),"intercept")
Armin Rauschenberger's avatar
Armin Rauschenberger committed
278
279
280
281
282
283
284
285
    
    if(family=="gaussian" & type.measure=="deviance" | family=="binomial" & type.measure=="mse"){
      loss <- apply(X=pred,MARGIN=2,FUN=function(x) sum((y-x)^2))
    } else if(family=="binomial" & type.measure=="class"){
      loss <- apply(X=pred,MARGIN=2,FUN=function(x) mean(y!=x))
    } else {
      warning("Implement other loss functions!")
    }
Rauschenberger's avatar
Rauschenberger committed
286
    
Rauschenberger's avatar
Rauschenberger committed
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
    ### start temporary ###
    #stability <- numeric()
    #for(k in seq_along(fit)){
    #    matrix <- matrix(data=NA,nrow=5,ncol=5)
    #    for(i in seq_len(5)){
    #        for(j in seq_len(5)){
    #            a <- select[[i]][[k]]
    #            b <- select[[j]][[k]]
    #            matrix[i,j] <- length(intersect(a,b))/length(union(a,b))
    #        }
    #    }
    #    diag(matrix) <- NA
    #    stability[k] <- mean(matrix,na.rm=TRUE)
    #}
    #cat(stability)
    ### end temporary ###
    
Rauschenberger's avatar
Rauschenberger committed
304
305
306
307
    if(plot){
        graphics::par(mar=c(3,3,1,1))
        col <- rep(x=0,times=length(loss)-1)
        col[1] <- col[length(col)] <- 1
Armin Rauschenberger's avatar
Armin Rauschenberger committed
308
309
310
311
312
313
        graphics::plot.new()
        graphics::plot.window(xlim=c(1,length(loss)-1),ylim=range(loss))
        graphics::axis(side=1,at=seq_len(length(loss)-1),labels=names(loss)[-length(loss)])
        graphics::axis(side=2)
        graphics::box()
        graphics::points(y=loss[-length(loss)],
Rauschenberger's avatar
Rauschenberger committed
314
                       x=seq_len(length(loss)-1),
Armin Rauschenberger's avatar
Armin Rauschenberger committed
315
                       col=col+1,pch=col)
Rauschenberger's avatar
Rauschenberger committed
316
317
318
        graphics::abline(v=c(1.5,length(loss)-1.5),lty=2)
        graphics::grid()
        graphics::abline(h=loss[length(loss)],lty=2,col="red")
Rauschenberger's avatar
Rauschenberger committed
319
320
    }
    
Rauschenberger's avatar
Rauschenberger committed
321
    return(loss)
Rauschenberger's avatar
Rauschenberger committed
322
}
Armin Rauschenberger's avatar
Armin Rauschenberger committed
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394



#' @export
#' @title
#' moderated response (trial)
#'
#' @description
#' This function ...
#'
#' @inheritParams colasso
#'
#' @examples
#' NA
#' 
moderate <- function(y,Y,k=2){
  p <- ncol(Y)
  id <- which(apply(X=Y,MARGIN=2,FUN=function(x) all(y==x)))
  if(length(id)!=1){stop("Invalid input.")}
  if(k==1){
    cluster <- rep(x=1,times=p)
  } else if(k==ncol(Y)){
    cluster <- seq_len(p)
  } else {
    cluster <- stats::kmeans(x=t(Y),centers=k)$cluster
  }
  cond <- cluster == cluster[id]
  median <- apply(Y[,cond,drop=FALSE],1,median)
  message("cor ",round(cor(y,median),2))
  return(median)
}

# moderate <- function(y,Y){
#   id <- which(apply(X=Y,MARGIN=2,FUN=function(x) all(y==x)))
#   if(length(id)!=1){stop("Invalid input.")}
#   
#   n <- nrow(Y)
#   p <- ncol(Y)
#   
#   
#   kmeans <- stats::kmeans(x=t(Y),centers=3)
#   
#   ## hierarchical clustering
#   #d <- stats::dist(x=t(Y),method="euclidean")
#   ##d <- as.dist(1-abs(stats::cor(data,method="spearman")))
#   #tree <- stats::hclust(d=d,method="complete")
#   #stats:::plot.hclust(x=tree,labels=FALSE)
#   #cut <- stats::cutree(tree=tree,k=seq_len(p))
#   #cond <- apply(X=cut,MARGIN=2,function(x) x==x[id])
#   #median <- apply(X=cond,MARGIN=2,FUN=function(x) apply(Y[,x],1,median))
#   
#   ##old
#   # cluster <- matrix(data=NA,nrow=p,ncol=p)
#   # list <- list()
#   # for(i in seq_len(p)){
#   #   list[[i]] <- matrix(data=NA,nrow=n,ncol=i)
#   #   id <- stats::cutree(tree=tree,k=i)
#   #   for(j in seq_len(i)){
#   #     list[[i]][,j] <- apply(X=data[,id==j,drop=FALSE],
#   #                            MARGIN=1,function(x) median(x))
#   #   }
#   # }
#   
#   cor <- lapply(list,function(x) abs(stats::cor(y,x)))
#   id_cluster <- lapply(cor,function(x) which.max(x))
#   max <- sapply(cor,function(x) max(x))
#   id_depth <- min(which(max>0.5))
#   mod <- list[[id_depth]][,id_cluster[[id_depth]]]
#   plot(y,mod)
#   return(mod)
# }