用量反応曲線から決定係数を出す

こちらの続き。
決定係数を出したかったのだが、drm関数中にはなさそうなので自作する。
 
n個のデータがあって、i=1,2,\dots,nとすると
R^2:決定係数
y_i:実測値
f_i:回帰による推定値
\widebar{y}:観測値の平均
より
R^2=1-\frac{\sum_i(y_i-f_i)^2}{\sum_i(y_i-\widebar{y})^2}
と定義されるらしい。
パラメータの数により、上の決定係数の定義は良くなる傾向があるらしく、
N:標本数
p:パラメータの数
により、以下の補正をすることもあるらしい。
R^2=1-\frac{\frac{\sum_i(y_i-f_i)^2}{(N-p-1)}}{\frac{\sum_i(y_i-\widebar{y})^2}{(N-1)}}

# データは上のリンクから

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)