I am working on a dataset where I need to fit a non-linear model to observed data and perform a permutation test to check the significance of the fit. The dataset consists of measurements related to mean values (mean.val), variance (var.val), and a calculated variable (VR.val). The goal is to fit the model VR.val ~ (var.val / (mean.val^2)) - (1 / mean.val) and then perform a permutation test to determine if the observed relationship is significant.

I am encountering an error when trying to fit the non-linear model using the nls function in R. Specifically, the error message is:

Error in model fitting: dims 
 do not match the length of object [72]

Here is the code I am using to fit the model. You can run it online here: https://rdrr.io/snippets/

# Load necessary libraries
library(dplyr)
library(nls2)

# Create the dataset
df <- tibble(
  mean.val = c(4.0908376, 5.1765852, 88.6562224, 2.2583648, 0, 0, 1.7031457, 11.3987297, 1.8400311, 0.4573476, 6.7048651, 48.6677506, 0.3508331, 15.8838267, 5.4980611, 4.1961417, 0, 27.9447267, 161.4052832, 65.2066269, 229.9773003, 80.5368571, 0.8070257, 44.522054, 26.7971708, 18.9939169, 42.9711176, 223.5922348, 2.6659812, 75.3845185, 229.0444097, 16.3551315, 35.483228, 14.5498903, 0, 96.8871753, 8.9939657, 23.3706346, 47.894185, 44.3451732, 0, 12.0112976, 56.6164317, 24.1342521, 25.5272703, 16.7003554, 3.5707173, 11.8552761, 109.6624751, 3.1851156, 165.5078472, 25.8925659, 0, 7.0753578, 128.7662599, 6.1903777, 16.0177727, 17.1843838, 0, 4.689188, 0, 0, 5.3534193, 0, 0, 0, 4.7981103, 9.7982937, 0.5867907, 23.1950658, 0, 0),
  var.val = c(4.697147e+01, 1.071881e+02, 1.168461e+04, 1.066658e+01, 0, 0, 1.188319e+00, 4.664088e+02, 6.717956e+00, 8.366671e-01, 1.579417e+02, 2.025821e+02, 4.923355e-01, 6.317895e+02, 4.581204e+01, 7.043042e+01, 0, 1.050463e+03, 3.413311e+04, 5.354310e+03, 5.974960e+04, 1.030781e+04, 2.605162e+00, 1.250515e+03, 2.038420e+03, 8.177155e+02, 3.283674e+03, 1.939515e+05, 2.842982e+01, 1.656980e+04, 1.549785e+05, 8.067766e+02, 3.935897e+03, 1.989293e+02, 0, 2.794197e+04, 1.369383e+02, 1.387005e+03, 2.509730e+03, 2.796107e+03, 0, 2.678243e+02, 3.825565e+03, 8.742454e+02, 9.923611e+02, 4.794730e+02, 5.100009e+01, 8.238411e+01, 1.597242e+04, 4.057985e+01, 9.698418e+04, 2.681700e+03, 0, 2.630001e+01, 4.747850e+04, 1.532831e+02, 7.726944e+02, 8.294296e+02, 0, 1.720190e+01, 0, 0, 1.146364e+02, 0, 0, 0, 5.712195e+01, 6.559322e+01, 1.377293e+00, 2.031250e+03, 0, 0),
  VR.val = c(2.56233945, 3.80682246, 1.4753261, 1.6486021, 0, 0, -0.1774834, 3.50193516, 1.4407373, 1.81347905, 3.36416676, 0.06498251, 1.14964199, 2.44120316, 1.33363381, 3.76168584, 0, 1.30939697, 1.30401283, 1.24393769, 1.12535676, 1.576778, 2.76088209, 0.60840734, 2.8013587, 2.21394245, 1.75503535, 3.8750636, 3.62490359, 2.90250193, 2.94978082, 2.95495357, 3.09787915, 0.87094969, 0, 2.96630615, 1.58167963, 2.49664515, 1.0732319, 1.3993234, 0, 1.77313892, 1.17580472, 1.45951317, 1.48368978, 1.65926691, 3.71994423, 0.50181467, 1.31905401, 3.68603965, 3.53445099, 3.96137888, 0, 0.38402702, 2.855705, 3.83845897, 2.94921243, 2.75054803, 0, 0.56905734, 0, 0, 3.8132035, 0, 0, 0, 2.27278903, 0.58115744, 2.29581483, 3.73236742, 0, 0)
)

# Vectorized function to calculate VR values
calculate_VR <- function(mean_val, var_val) {
  ifelse(mean_val == 0 | var_val == 0, 0, (var_val / (mean_val^2)) - (1 / mean_val))
}

# Recalculate VR values to follow the theoretical relationship
df <- df %>%
  dplyr::mutate(VR.val = calculate_VR(mean.val, var.val))

# Fit non-linear model using nls
fit_model <- function(data) {
  tryCatch({
    model <- nls(VR.val ~ I(var.val / (mean.val^2)) - I(1 / mean.val), data = data,
                 start = list(mean.val = mean(data$mean.val), var.val = mean(data$var.val)),
                 algorithm = "port")
    model_summary <- summary(model)
    
    # Calculate RMSE
    predicted <- predict(model)
    residuals <- data$VR.val - predicted
    rmse <- sqrt(mean(residuals^2))
    
    list(
      coefficients = coef(model),
      RMSE = rmse
    )
  }, error = function(e) {
    cat("Error in model fitting:", e$message, "n")
    return(list(coefficients = NA, RMSE = NA))
  })
}

# Fit the model to observed data
observed_fit <- fit_model(df)

I am not very sure how to correct this. I tried many approaches and crosschecking the structure but nothing seems to work.