線形モデルの予測値を導くためには、それなりに長いスクリプトを書く必要があり面倒です（参考）。また、lmerのようにランダム効果を入ってくると、信頼区間を推定するためのスクリプトが一気に煩雑になります。そこで、予測の平均値およびその信頼区間を計算する関数predictRを作りました。第一引数にモデルオブジェクトをいれ、第二引数に予測に用いたい説明変数Xの名前を指定するだけで、Xの最小値から最大値まで動かしたときの予測値と信頼区間を返します。なお、重回帰モデルの場合、第二引数以外の変数は平均（連続値）もしくは一番名前の若いグループ名（因子型）に固定されます。オフセット項についてもoffset=Tとすることでユニット当たりの予測値を返します。
関数の実装には、以下のスクリプトをRコンソールにコピペします。
predictR <- function(m, f_var, linkinv = T, ci.level=0.95, offset=F,
Interaction = F, int_var=NULL, pct = NULL,
n=100){
## Preparation for a new data set---
# model class
CL <- class(m)[1]
# data frame
dat <- model.frame(m)
# Random effect
if(CL=="lmerMod"|CL=="glmerMod"){
REID <- NULL
for(i in 1:length(names(ranef(m)))){
tmp <- which(colnames(dat)==names(ranef(m))[i])
REID <- c(REID, tmp)
}
dat <- dat[,-REID]
}
# offset term
fixID <- which(colnames(dat)!="(offset)")
ofsID <- which(colnames(dat)=="(offset)")
dat <- dat[,fixID]
if(CL=="lmerMod"|CL=="glmerMod"){
dat.check <- length(names(dat)[-1]) - length(names(fixef(m))[-1])
}else{
dat.check <- length(names(dat)[-1]) - length(names(coefficients(m))[-1])
}
if(offset==T){
if(length(ofsID)==0|dat.check>0){ stop("Specify offset term as e.g., 'offset = log(x)' in the model object") }
}else{
if(length(ofsID)!=0|dat.check>0){ stop("Should be 'offset=T' and specify as e.g., 'offset = log(x)' in the model object") }
}
# turn "character" or "logical" into "factor"
dat <- as.data.frame(lapply(dat,function(x){
if(is.character(x)|is.logical(x)){
as.factor(x)
}
else{x}
}))
# choose focal variable
fID <- which(names(dat)==f_var)
if(is.factor(dat[,fID])){
stop("focal variable must be numeric")
}else{
X1 <- list(seq(min(dat[,fID]),max(dat[,fID]), length = n))
names(X1) <- f_var
}
# Interction terms
if(Interaction == T){
iID <- which(names(dat)==int_var) # column ID of interction terms
if(length(iID)>1){ stop("Only one interaction term is allowed") }
if(is.null(int_var)==T){ stop("Specify interaction terms") }
if(is.factor(dat[,iID])==1){#if factor
if(length(iID)==1){ X2 <- list(unique(dat[,iID])); names(X2) <- int_var }
}else{#if numeric
if(length(iID)==1){ X2 <- list(quantile(dat[,iID],pct)); names(X2) <- int_var }
if(is.null(pct)==T){ stop("Specify percentiles of interaction term (ex. pct=c(0.2,0.8))") }
}
X <- c(X1,X2)
X <- expand.grid(X)
}else{
if(length(int_var)>0){ stop("Should be 'Interaction = T'") }
X <- X1
}
# make new data set
if(Interaction==T){
dat_r <- as.data.frame(dat[,-c(1,fID,iID),drop=F])
if(dim(dat_r)[2]==0){
new_dat <- X
}else{
fixed <- rep(lapply(dat_r,function(x) if(is.numeric(x)) mean(x) else factor(levels(x)[1],levels = levels(x))))
new_dat <- cbind(X,as.data.frame(fixed))
}
}else{
dat_r <- as.data.frame(dat[,-c(1,fID),drop=F])
if(dim(dat_r)[2]==0){
new_dat <- X
}else{
fixed <- rep(lapply(dat_r,function(x) if(is.numeric(x)) mean(x) else factor(levels(x)[1],levels = levels(x))))
new_dat <- cbind(X,as.data.frame(fixed))
}
}
## Get predicted values---
# interval
cr <- ci.level + (1-ci.level)/2
CR <- qnorm(cr,0,1)
# inverse link function
ilink <- ifelse(linkinv==T, family(m)$linkinv, identity)
if(Interaction==T){
output <- list(NULL)
Index <- ifelse(is.factor(new_dat[,int_var]), length(unique(new_dat[,int_var])), length(pct))
for(i in 1:Index){
st <- I(1+(i-1)*n); en <- i*n
if(CL == "lm"|CL == "glm"|CL == "negbin"){
fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
mm <- model.matrix(fo, new_dat[st:en,])
us.fit <- mm%*%coefficients(m)
se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
fit <- ilink(us.fit)
uci <- ilink(us.fit+CR*se)
lci <- ilink(us.fit-CR*se)
y.tmp <- data.frame(fit=fit,uci=uci,lci=lci)
}
if(CL == "lmerMod"| CL == "glmerMod"){
fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
mm <- model.matrix(fo, new_dat[st:en,])
us.fit <- mm%*%fixef(m) # unscaled fitted value
se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
fit <- ilink(us.fit)
uci <- ilink(us.fit + CR*se)
lci <- ilink(us.fit - CR*se)
y.tmp <- data.frame(fit=fit,uci=uci,lci=lci)
}
output[[i]] <- data.frame(y.tmp, new_dat[st:en,])
}#i
u_int_var <- unique(new_dat[,int_var])# unique values or name for interaction terms
if(is.factor(new_dat[,int_var])==T){
names(output) <- as.character(u_int_var)
}else{
names(output) <- as.factor(paste(pct*100,"%","_",int_var,sep=""))
}
return(output)
}else{
if(CL == "lm"|CL == "glm"|CL == "negbin"){
fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
mm <- model.matrix(fo, new_dat)
us.fit <- mm%*%coefficients(m) # unscaled fitted value
se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
fit <- ilink(us.fit)
uci <- ilink(us.fit + CR*se)
lci <- ilink(us.fit - CR*se)
Y <- data.frame(fit=fit,uci=uci,lci=lci)
}
if(CL == "lmerMod"| CL == "glmerMod"){
fo <- as.formula(paste("~",as.character(formula(m,fixed.only=T)[3]))) # formula
mm <- model.matrix(fo, new_dat)
us.fit <- mm%*%fixef(m) # unscaled fitted value
se <- sqrt(diag(mm%*%tcrossprod(vcov(m),mm)))
fit <- ilink(us.fit)
uci <- ilink(us.fit + CR*se)
lci <- ilink(us.fit - CR*se)
Y <- data.frame(fit=fit,uci=uci,lci=lci)
}
output <- data.frame(Y,new_dat)
return(output)
}
}
引数は以下の通りです。
m: モデルオブジェクト(lm, glm, glm.nb, lmer)
f_var: 予測を行いたい説明変数の名前
linkinv: 逆リンク関数による応答変数の変換（デフォルトT）
ci.level: 信頼区間のレベル（デフォルト0.95）
offset: offset項の有無（デフォルトF）
Interaction: f_varとの交互作用項の有無（デフォルトF）
int_var: 交互作用項の名前（ひとつだけ指定可能）
pct: 交互作用項が連続値の場合、予測に用いる分位点（デフォルトNULL）
n: 予測に用いる説明変数の長さ（デフォルト100）
#lm, glm, glm.nb, lmer, glmer
#データを用意する
> D
y1 x1 x2 x3
1 1 -0.2469068 19 C
2 3 0.7600840 18 T
3 3 0.6211508 16 C
4 3 -0.4965345 18 T
5 4 0.5435650 19 C
6 5 0.5004225 23 T
m <- glm(y1~x1+x2+x3, D, family=poisson)
fit <- predictR(m, f_var="x1")
> head(fit)
fit uci lci x1 x2 x3
1 2.351697 3.517449 1.185944 -1.610692 19.92 C
2 2.352722 3.501309 1.204135 -1.568597 19.92 C
3 2.353747 3.485274 1.222221 -1.526502 19.92 C
4 2.354772 3.469348 1.240196 -1.484408 19.92 C
5 2.355798 3.453539 1.258057 -1.442313 19.92 C
6 2.356823 3.437849 1.275796 -1.400218 19.92 C
#fit: 予測値の平均
#uci,lci: ci.levelで設定した信頼区間の上限・下限
#x1,x2,x3: 予測値を導く際の説明変数の代入値
#matplotで図示
matplot(fit[,"x1"], fit[,1:3], type="l", col="black")
#交互作用項（int_var）が因子の場合
m <- glm(y1~x1*x3, D, family=poisson)
fit <- predictR(m, f_var="x1", Interaction = T, int_var = "x3")
> head(fit$C,5)
fit uci lci x1 x3
1 1.295155 2.358462 0.7112373 -1.721629 C
2 1.318219 2.377811 0.7307986 -1.676699 C
3 1.341694 2.397378 0.7508794 -1.631769 C
4 1.365586 2.417169 0.7714919 -1.586839 C
5 1.389905 2.437190 0.7926487 -1.541909 C
> head(fit$T,5)
fit uci lci x1 x3
101 1.087035 1.653718 0.7145378 -1.721629 T
102 1.119484 1.692855 0.7403141 -1.676699 T
103 1.152902 1.732966 0.7669992 -1.631769 T
104 1.187318 1.774078 0.7946234 -1.586839 T
105 1.222761 1.816220 0.8232178 -1.541909 T
#交互作用項（int_var）が連続値の場合：int_varの分位点を指定
m <- glm(y1~x1*x2, D, family=poisson)
fit <- predictR(m, f_var="x1", Interaction = T, int_var = "x2", pct=c(0.2, 0.8)
> head(fit$`20%_x2`,5)
fit uci lci x1 x2
1 1.000534 1.678020 0.5965770 -1.721629 16
2 1.026491 1.707934 0.6169346 -1.676699 16
3 1.053121 1.738428 0.6379697 -1.631769 16
4 1.080443 1.769517 0.6597036 -1.586839 16
5 1.108473 1.801214 0.6821577 -1.541909 16
> head(fit$`80%_x2`,5)
fit uci lci x1 x2
101 1.248705 2.032481 0.7671726 -1.721629 23
102 1.277472 2.063563 0.7908334 -1.676699 23
103 1.306902 2.095166 0.8152064 -1.631769 23
104 1.337010 2.127301 0.8403113 -1.586839 23
105 1.367811 2.159981 0.8661687 -1.541909 23
※交互作用のないモデルでInteraction = Tとすると、int_varで指定した変数の相加効果のみが反映されます。
※式中にlog(x)のような関数が入るとエラーがでます。一度別名変数として格納することで実行可能です。
自作関数になりますので、誤り等がある可能性があります。各自の責任のもとでご利用ください。エラーなどありましたらコメントいただけると幸いです。
SPONSOR LINK