Extracting coefficients of non-linear equation across categorical variable

Question

This question appears to be off-topic because it focuses on programming, debugging, or performing routine operations, or it asks about obtaining datasets. You could try the support links we maintain or the Open Data site instead.

If you are asking about how to debug some code or carry out a task in a programming language, this is not an on-topic question. If you're asking about how to obtain data, that is not an on-topic question. If you feel your question is truly about statistics as described in the help center, please [edit] to clarify.

Closed 25 mins ago.

(Private feedback for you)

I am trying to extract the coefficients k1,2 and 3 from the equation:

Mt=M1exp(−k1⋅CDI⋅t)+M2exp(−k2⋅CDI⋅t)+M3exp(−k3⋅CDI⋅t)

where Mt is mass at time t, m1 is the initial labile carbon content (in %), m2 is holocellulose content, m3 is lignin content, and finally, CDI

is climate decomposition index (see below image).

I have managed to extract the coefficients from a single sight using the following code:

eqtn <- function(m1, k1, cdi, t, m2, k2, m3, k3){(m1 * exp(-k1 * cdi * t)+
                                             m2 * exp(-k2 * cdi * t)+
                                             m3 * exp(-k3 * cdi * t))}

nls(mass_remaining_percent ~ eqtn(scf_mean_initial, k1, cdi_mean, days_between, 
                                  holocellulose_mean_initial, k2, lignin_mean_initial, k3),
    start = list(k1 = 0.0007, k2 = 0.0005, k3 = 0.0001), data = a.3_pooled_data)

Does anyone know how I can apply this over categories of sight (field code) and extract the coefficients/ R2? I know that it needs splitting into many small datasets, then applying the model, extracting data and recombining but I can't figure out how to do it

Thanks in advance

I have just managed to write the above code, fitting the model to a singular line of code, I know that I can do this manually, fitting to many smaller datasets but it will take an age and I would like to learn how to do it using the apply family or the broom package.

Answer 1

For coefficients, use either coef to return a vector or summary.nls(...)$coefficients to return a matrix which also includes t-stats and p-values on results of nls(). To run model by groupings, consider by (the object-oriented wrapper of tapply):

eqtn <- function(m1, k1, cdi, t, m2, k2, m3, k3) {
    (m1 * exp(-k1 * cdi * t) +
     m2 * exp(-k2 * cdi * t) +
     m3 * exp(-k3 * cdi * t))
}

model <- mass_remaining_percent ~ eqtn(
    scf_mean_initial, k1, cdi_mean, days_between, 
    holocellulose_mean_initial, k2, lignin_mean_initial, k3
)

run_model_by_group <- function(sub) {
    results <- nls(
        model,
        start = list(k1 = 0.0007, k2 = 0.0005, k3 = 0.0001), 
        data = sub
    )

    return(summary(results)$coefficients)
}

# LIST OF COEFFICIENT MATRICES BY FIELD CODE
field_code_coefficient_matrices <- by(
    a.3_pooled_data,
    a.3_pooled_data$Field.Code,
    FUN = run_model_by_group
)