Modeling a Physiological Process

[Brandon Hurr]

Dear D-RUGgers, 

I'll preface with the fact that I'm into R for the pretty pictures and data munging and am not a statistician. My eyes typically roll into my skull when I see math equations. But, I see this as a learning opportunity.

I'm reviewing a manuscript that models weight loss of a commodity postharvest under a treated condition. They used excel (😕) to do a simple linear regression of weight loss over time. I'm fine with that and they attempt to show goodness of fit with r^2 and something called the mean relative percentage deviation modulus. 

They have multiple temperatures at which they are measuring the effect of the treatment. 

Something like this (note: example is only a single "treatment"):

require(dplyr)

require(broom)

require(ggplot2)

df1 <- data.frame(time = seq(0, 60, 4), temp = rep(c(15,20,25,30), each=16), pctwgtloss = c(seq(0, 15, length.out= 16), seq(0, 16, length.out= 16), seq(0, 18, length.out= 16), seq(0, 20, length.out= 16)), error = runif(16*4, min = -2, max = 2))

arrhenius <- 

df1 %>%

mutate(pctloss = pctwgtloss + error) %>% #add in some error to the data

select(-pctwgtloss, -error) %>% #get rid of error columns

mutate(pctloss = ifelse(time 0= 0, 0, pctloss)) %>% #can't have anything but 0 at time 0

group_by(temp) %>% 

do(fit = lm(pctloss ~ 0 + time, data=.)) %>% #linear model forced through origin

tidy(fit) %>% #clean up model output with broom

select(temp, estimate) %>% #select model parameters for next level of model

ungroup() %>% #get rid of groupings

mutate(temp = 1/(temp + 273.15)) %>% #convert T into kelvin

mutate(estimate = -log(estimate)) %>% #convert k in to -ln(k)

do(fit2 = lm(estimate ~ temp, data=.)) %>% #natural log of k estimate ~ 1/T

tidy(fit2) %>% #clean up model output with broom

select(estimate) 

#save away second model parameter estimates

K = exp(-arrhenius$estimate[1])

E = arrhenius$estimate[2] * 8.31 #unsure why, but Ea contains R

finalplot <- 

data.frame(temp = seq(15, 30, 1)) %>% #all possible temps between 15 and 30 in 1 C increments

mutate(temp = temp + 273.15)  #convert T into kelvin

finalplot %>%

mutate(t = 9/K * exp(E/8.31/temp)) %>% #calculate t for each temp

mutate(temp = temp - 273.15) %>% #mutate temp back to C from K

ggplot(., aes(x=temp, y=t)) + #plot it!

geom_point() + #points

stat_smooth() + #smooth line so you can see the slight curve

labs(x="Temp (C)", y="Time to Unacceptable (h)") #pretty labels for reading

What they do is they fit the data (per temperature * treatment) to a linear model, take the slope coefficient, plot those [-log(slope)] by 1/temperature use those to determine the activation energy Ea (slope) and the K value (e^-intercept) for the two treatments according to the Arrhenius Equation. 

This was put into a model (equation) that predicts the storability of the commodity. 

M (weightloss cutoff) = K * e^(-Ea/RT)*t 

where T is temperature in Kelvin, R is a constant (8.31 kJ/mol(K), t is time, and K and Ea are derived from the slope and intercept above. 

In this example with a 9% weight loss cutoff, 

t = 9/55.85 * e^(1560*8.31/(8.31*T))  #I really don't understand how Ea needs R, but you don't get back the right numbers without it. 

This allows you to predict under a given treatment when you would expect that commodity to fail (not pass inspection) based upon the weight loss parameter alone. 

I don't see anything inherently wrong with what they've done, but wanted some feedback from people who might know better. By the time you get to the second model there is no "error", and the measure of goodness of fit is the r^2 and mean relative percentage deviation modulus, which looks an awful lot like an SSE. 

Any ideas/suggestions on how to improve the modeling/stats? 

Thanks

Brandon