こちらの続き。
決定係数を出したかったのだが、drm関数中にはなさそうなので自作する。
n個のデータがあって、とすると
:決定係数
:実測値
:回帰による推定値
:観測値の平均
より
と定義されるらしい。
パラメータの数により、上の決定係数の定義は良くなる傾向があるらしく、
:標本数
:パラメータの数
により、以下の補正をすることもあるらしい。
# データは上のリンクから a<- 80/22 conc <- c(0.01,0.03,0.1, 0.3,1, 3,10, 30,100,300,1000) # Dose of isoproterenol pro0 <- c(1,1,1,7,18,22,31,33,33,34,33)*a pro10 <- c(1,1,1,6,17,26,31,33,33,33,33)*a pro30 <- c(0,0,0,4,15,24,30,32,32,32,32)*a pro100 <- c(0,1,1,1,11,21,29,32,33,33,33)*a pro300 <- c(1,1,1,1,4,15,25,30,32,33,33)*a pro1000<- c(1,1,1,1,1,5,16,26,31,32,33)*a pro3000<- c(1,1,1,1,1,1,7,18,27,31,32)*a pro.cont<- c(0,10,30,100,300,1000,3000) EC50<- c(0.8442,0.9240,1.116,1.968,4.064,11.79,30.89) library(drc) data<- rbind(pro0,pro10,pro30,pro100,pro300,pro1000,pro3000) ldata<- matrix(0,nrow(data),4) res_list<- lapply(mapply(rep,0,nrow(data)),unique) for(p in 1:nrow(data)){ res_list[[p]]<- drm(data[p,]~conc,fct = LL.4()) ldata[p,]<- res$start plot(res_list[[p]],col=p,ylim=c(0,130),type="none") # typeを指定しなければ観測データもプロットする。 if(p!=nrow(data)){ par(new=TRUE) } }
r2_coef<- function(drm_result){ if(is(drm_result) != "drc"){ return("Not drm") }else{ No_parm<- length(drm_result$fit$par) N<- nrow(drm_result$data) resp_pred<- drm_result$predres[,"Predicted values"] resp_data<- drm_result$data[,2] coef<- 1-sum((resp_data-resp_pred)^2)/sum((resp_data-mean(resp_data))^2) adj_coef<- 1-(sum(((resp_data-resp_pred)^2))/(N-No_parm-1))/(sum(((resp_data-mean(resp_data))^2))/(N-1)) } res<- list(coef,adj_coef) names(res)<- list("coef","adjusted") return(res) }
lapply(res_list,r2_coef)