如何从R语言自定义模型拟合结果中计算R-squared值

要评估fitPGH模型的拟合效果，你可以通过两种方式计算R-squared值：手动计算或修改函数自动返回结果。

方法一：拟合后手动计算 #

在现有拟合代码的基础上，添加以下步骤：

# 现有拟合代码
PAR <- kelppercent$averagePAR
Pc <- kelppercent$Kelp_Cover
myfit <- fitPGH(PAR, Pc)

# 1. 获取模型实际使用的过滤后数据（移除非有限值）
ind <- is.finite(PAR) & is.finite(Pc)
y_observed <- Pc[ind]

# 2. 计算模型对过滤后数据的预测值
y_predicted <- with(myfit, {
  ps[1] * (1 - exp(-alpha[1] * PAR[ind] / ps[1])) * exp(-beta[1] * PAR[ind] / ps[1])
})

# 3. 计算R-squared（两种等价方式）
# 方式A：利用残差平方和(SSR)与总平方和(SST)
ssr <- myfit$ssr
sst <- sum((y_observed - mean(y_observed))^2)
r_squared <- 1 - (ssr / sst)

# 方式B：观测值与预测值的相关系数平方
r_squared <- cor(y_observed, y_predicted)^2

# 输出结果
cat("R-squared值：", round(r_squared, 4), "\n")

方法二：修改fitPGH函数自动返回R-squared #

修改原函数，使其在拟合成功后自动计算并返回R-squared值，方便后续调用：

fitPGH <- function(x,                          #E 
                   y,                          #Quantum Efficiency, rETR or P
                   normalize=FALSE,            #Should curve be normalized to E (Default=TRUE for modeling Quantum Efficiency)
                   lowerlim=c(-Inf),          #Lower bounds of parameter estimates (alpha,Beta,Ps)
                   upperlim=c(Inf),  #Upper bounds of parameter estimates (alpha,Beta,Ps)
                   fitmethod=c("Nelder-Mead")) #Fitting method passed to modFit  
{
  
  #If normalize =T, assign E = 0 to very small number
  if (normalize==T)  x[x==0] <- 1e-9       
  
  #Remove NA values
  ind   <- is.finite(x) & is.finite(y)
  res   <- rep(NA,length(x))
  x     <- x[ind==T]
  y     <- y[ind==T]
  
  #Intitial Parameter Estimates
  if (normalize==T){ 
    alpha <- max(y)
    beta  <- 0
    ps    <- max(x*y)
  }
  if (normalize==F){ 
    PE    <- y/x
    alpha <-  max(PE[is.finite(PE)])
    beta  <- 0
    ps    <- max(y)
  }
  
  #Load the model
  PGH     <- function(p,x) return(data.frame(x = x, y = p[3]*(1-exp(-1*p[1]*x/p[3]))*exp(-1*p[2]*x/p[3])))
  PGH.E   <- function(p,x) return(data.frame(x = x, y = p[3]*(1-exp(-1*p[1]*x/p[3]))*exp(-1*p[2]*x/p[3])/x))
  if (normalize==F) model.1 <- function(p) (y - PGH(p, x)$y)
  if (normalize==T) model.1 <- function(p) (y - PGH.E(p, x)$y)
  
  
  #In case of non-convergence, NAs are returned
  if (class(try(modFit(f = model.1,p = c(alpha,beta,ps),method = fitmethod, 
                       lower=lowerlim,upper=upperlim, 
                       hessian = TRUE),silent=T))=="try-error"){
    fit <- list(alpha=NA,beta=NA,ps=NA,ssr=NA,residuals=rep(NA,c(length(x))),
                R_squared=NA, y_obs=NA, y_pred=NA)
  }else{
    fit <- modFit(f = model.1,p = c(alpha,beta,ps),method = fitmethod, 
                  lower=lowerlim,upper=upperlim, hessian = TRUE)
    
    # 计算预测值与R-squared
    p <- fit$par
    y_pred <- if(normalize) PGH.E(p, x)$y else PGH(p, x)$y
    sst <- sum((y - mean(y))^2)
    r_squared <- 1 - (fit$ssr / sst)
    
    fit <- list(alpha=summary(fit)$par[1,],beta=summary(fit)$par[2,],ps=summary(fit)$par[3,],
                ssr=fit$ssr,residuals=fit$residuals,model="PGH",normalize=normalize,
                R_squared=r_squared, y_obs=y, y_pred=y_pred)
  }
  
  return(fit)
  
}

修改后，拟合时直接获取R-squared：

myfit <- fitPGH(PAR, Pc)
cat("R-squared值：", round(myfit$R_squared, 4), "\n")

说明 #

R-squared的核心逻辑是1 - (残差平方和/总平方和)，反映模型解释数据变异的比例。两种计算方式等价，你可以根据需求选择。

内容的提问来源于stack exchange，提问作者SecretBeach