Методы регрессии для данных с интервалом для х и у - PullRequest
Новая
       66

Методы регрессии для данных с интервалом для х и у

0 голосов
/ 07 февраля 2020

У меня есть набор данных, в котором как зависимые (y), так и независимые переменные (x) представляют собой интервалы, а не точечные данные. Интервалы проистекают из некоторой неопределенности, то есть истинное значение находится где-то в пределах интервала (для простоты, предполагая равномерную вероятность в пределах интервала).

Теперь я хотел бы провести некоторый регрессионный анализ (на самом деле смешанный модель, которая также включает в себя случайные величины). Я хотел бы рассчитать средние коэффициенты регрессии и некоторые оценки достоверности. Кроме того, я хотел бы получить интервалы прогнозирования для нового набора данных для построения отношения y~x.

Поэтому я подумал, что симуляция Монте-Карло, то есть повторная выборка точек из их интервалов и последующие (повторяющиеся) модели, может быть способом решения этой задачи; Точнее говоря, я думал об использовании lmer() из R-пакета lme4 для определения симуляции как дополнительных случайных факторов ... Однако я не уверен, действительно ли это правильный способ сделать это.

Поскольку мои реальные данные более сложны и могут включать в себя различные типы моделей (GLMM), я ищу довольно общий способ (такой как моделирование Монте-Карло и / или метод начальной загрузки и т. Д. c), который может применяется также в других типах моделей.

Вот некоторый код для создания примера набора данных: 

library(MASS)
library(ggplot2)
library(randomcoloR)

# Create dataframe with x and y
set.seed(999)
Sigma <- matrix(c(8,5,5,5),2,2)
df <- data.frame(mvrnorm(n = 30, mu=c(10, 10), Sigma))
colnames(df) <- c("y_mean","x_mean")

# Create to groups and increase y-value for one group
df$group <- sample(c("A","B"),nrow(df),replace=T)
df$y_mean[df$group=="B"] <- df$y_mean[df$group=="B"]+4
# Id/Subject
df$id=seq(1,nrow(df))

#Color for each id within each group
df$group_id <- paste(df$group,df$id,sep="_")
df$col <- "black"
df$col[grepl(".*A",df$group_id)] <- randomColor(count = sum(grepl(".*A",df$group_id)), hue = c("red"))
df$col[grepl(".*B",df$group_id)] <- randomColor(count = sum(grepl(".*B",df$group_id)), hue = c("blue"))


# Expand from x-y-points to ranges/intervals
set.seed(999)
delta_x <- runif(nrow(df),1,5)
delta_y <- runif(nrow(df),1,5)
df$x_min <- df$x_mean - delta_x/2
df$x_max <- df$x_mean + delta_x/2
df$y_min <- df$y_mean - delta_y/2
df$y_max <- df$y_mean + delta_y/2

Теперь для запуска моделирования для выборки из данных интервала и построения отдельных линий регрессии (см. plot of intervals): 

# Multiple simulations / Monte Carlo
R <- 20 #number of simulations
sim_ls <- list()
for(i in 1:R){
  x_sim_i <- apply(df, 1, function(x) runif(1,min=as.numeric(x["x_min"]),max=as.numeric(x["x_max"])))
  y_sim_i <- apply(df, 1, function(x) runif(1,min=as.numeric(x["y_min"]),max=as.numeric(x["y_max"])))
  sim_ls[[i]] <- data.frame(x=x_sim_i,
                            y=y_sim_i,
                            rep=paste("rep",i,sep=""),
                            id=df$id,
                            group=df$group,
                            col=df$col)
}
sim_df <- data.frame(do.call(rbind,sim_ls))
sim_df$rep <- as.factor(sim_df$rep)
sim_df$id <- as.factor(sim_df$id)
sim_df$group <- as.factor(sim_df$group)

# Plot as rectangles of intervals
# and regression lines for single simulations
ggplot() +
  geom_rect(data=df, 
            aes(xmin = x_min, xmax =x_max, ymin = y_min, ymax = y_max), 
            colour = NA, fill = df$col, alpha=0.3)+
  geom_point(data=df, aes(x= x_mean,y=y_mean),shape=1,size=3)+
  #scale_fill_viridis(discrete=T) +
  geom_smooth(data = sim_df, aes(x=x,y=y,group=rep), method='lm', formula= y~x,
              fill="grey",alpha=0.05,colour="grey20")

А теперь появилась идея использовать lmer с имитацией в качестве повторяющихся случайных факторов (вдохновлено https://github.com/lme4/lme4/issues/388#issuecomment -231398937 ) и создать plots with prediction/confidence intervals: 

library(lme4)
lmer_sim <- lmer(y~x+group+(1|rep),data=sim_df)
newDat_sim <- data.frame(x=seq(5,16,0.1),rep=999999,id=999999,group="B")
summary(newDat_sim)

## https://github.com/lme4/lme4/issues/388#issuecomment-231398937
## param, RE, and conditional
b1_sim <- bootMer(lmer_sim,FUN=function(x) {simulate(x,newdata=newDat_sim,re.form=~0,
                                                     allow.new.levels=TRUE,cond.sim=TRUE)[[1]]},
                  nsim=100,seed=101)
## param and RE (no conditional)
b2_sim <-  bootMer(lmer_sim,FUN=function(x) {simulate(x,newdata=newDat_sim,re.form=~0,
                                                      allow.new.levels=TRUE,cond.sim=FALSE)[[1]]},
                   nsim=100,seed=101)
## param only
b3_sim <- bootMer(lmer_sim,FUN=function(x) {predict(x,newdata=newDat_sim,re.form=~0)},
                  nsim=100,seed=101)

#### Confidence and prediction intervals for *unobserved* levels
bootsum <- function(x,ext="_1") {
  d <- t(data.frame(apply(x$t,2,
                        function(x) c(mean(x),quantile(x,c(0.025,0.975))))))
  colnames(d) <- paste0(c("bpred","lwr","upr"),ext)
  rownames(d) <- NULL
  return(d)
}
dd_sim <- cbind(newDat_sim,
                 bootsum(b1_sim),
                 bootsum(b2_sim,"_2"),
                 bootsum(b3_sim,"_3"))

# Plot data and prediction/confidence intervals
ggplot()+
  geom_rect(data=df, 
            aes(xmin = x_min, xmax =x_max, ymin = y_min, ymax = y_max), 
            colour = NA, fill=df$col, alpha=0.3)+
  geom_point(data=sim_df,aes(x=x,y=y),colour=sim_df$col,alpha=0.5)+
  ## param, RE, and conditional
  geom_ribbon(data=dd_sim,colour=NA, 
              aes(x=x,ymin=lwr_1,ymax=upr_1),alpha=0.2)+
  ## param and RE (no conditional)
  geom_ribbon(data=dd_sim,colour=NA,
              fill="blue",
              aes(x=x,ymin=lwr_2,ymax=upr_2),alpha=0.2)+
  ## param only
  geom_ribbon(data=dd_sim,colour=NA,
              fill="red",
              aes(x=x,ymin=lwr_3,ymax=upr_3),alpha=0.3)
...