0

I have following dataframe:

df <- structure(list(y = c(0.82, 0.77, 0.46, 0.7, 0.82, 0.92, 0.84, 0.88, 0.86, 0.92, 0.91, 0.96, 0.91, 0.92, 0.89, 0.95, 0.95, 0.88, 0.92, 0.88, 0.94, 0.72, 0.9, 0.95, 0.96, 0.92, 0.94, 0.93, 0.93, 0.94, 0.93, 0.89, 0.94, 0.94, 0.91, 0.88, 0.96, 0.91, 0.9, 0.95, 0.83, 0.95, 0.92, 0.91, 0.86, 0.94, 0.93, 0.83, 0.87, 0.76), x = c(0, 0.03, 0.07, 0.1, 2.2, 2.18, 2.33, 2.48, 2.63, 2.77, 2.92, 3.07, 3.22, 3.37, 3.52, 3.66, 3.81, 3.96, 4.11, 4.16, 4.21, 4.26, 4.31, 4.36, 4.41, 4.46, 4.51, 4.55, 4.6, 4.65, 4.7, 4.75, 4.8, 4.85, 4.9, 4.96, 5.01, 5.07, 5.12, 5.18, 5.24, 5.29, 5.35, 5.4, 5.46, 5.51, 5.57, 5.27, 4.98, 4.68), z = c(1.54, 1.48, 1.51, 1.05, 1.29, 0.6, 1.03, 0.95, 0.98, 0.89, 0.81, 0.91, 0.31, 0.69, 0.17, 0.48, 0.51, 0.74, 0.79, 0.77, 0.69, 0.5, 0.75, 0.85, 0.77, 0.7, 0.66, 1.02, 0.69, 0.51, 0.63, 0.45, 0.46, 0.7, 0.74, 0.68, 0.72, 0.84, 0.5, 0.62, 0.32, 0.74, 0.52, 0.65, 1.07, 0.96, 1.03, 1.41, 1.88, 0.83)), row.names = c(1L, 2L, 3L, 4L, 43L, 44L, 45L, 46L, 47L, 48L, 49L, 50L, 51L, 52L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L, 69L, 70L, 71L, 72L, 73L, 74L, 75L, 76L, 77L, 78L, 79L, 80L, 81L, 82L, 83L, 84L, 85L, 86L, 87L, 88L), class = "data.frame")

I need to fit equation y ~ ax^bz^c for the data given above. I tried the code below:

Log transforms the above equation to log(y) ~ log(a) + blog(x) + clog(z), so, one can fit a linear model as below:

x <- log(df$x)
y = log(df$y)
z = log(df$z)
x[!(is.finite(x))] <- NA
y[!(is.finite(y))] = NA
z[!(is.finite(z))] = NA

# Model fitting
m <- lm(y~x+z, na.action = na.exclude)

coeff <- list(a = coef(m)[1], b = coef(m)[2], c = coef(m)[3])
# Prediction 
y_pred <- coeff[[1]] + coeff[[2]]*x + coeff[[3]]*z   # or predict(m)
CORR_1 <- cor(y,y_pred, use = "pairwise.complete.obs")

Converting back to original scale

x <- df$x
y = df$y
z = df$z
y_pred <- exp(coeff[[1]])*x^coeff[[2]]*z^coeff[[3]]
CORR_2 <- cor(y,y_pred, use = "pairwise.complete.obs")

I was expecting CORR_1 and CORR_2 to be same but their values are different. Why is that so? What is the best way to fit y ~ ax^bz^c?

Artem
  • 3,304
  • 3
  • 18
  • 41
raghav
  • 533
  • 2
  • 11
  • Maybe start with `m2 <- nls(y ~ a * x^b * z^c, df, start = list(a=1, b=1, c=1))`. – Rui Barradas Nov 19 '21 at 21:00
  • Why do results differ for linear model fitting for log(y) ~ log(a) + blog(x) + clog(z) and fitting using your model? Fundamentally both are same. – raghav Nov 19 '21 at 21:16
  • Only the intercept is really different. The differences in the other parameters estimates have to do with `nls` convergence. – Rui Barradas Nov 19 '21 at 22:57
  • There's an algebraic error. log(a^b^c) = b ^ c * log(a) -> loglog(a^b^c) = c*log(b) + log(log(a)). – Artem Nov 21 '21 at 20:20

0 Answers0