www = "http://www.biostatisticaumg.it/dataset/doublekinetics.csv"
doublekinetics = read.csv(www, header = TRUE)
# doublekinetics = read.csv(file.choose(), header = TRUE)
attach(doublekinetics)
head(doublekinetics)
tail(doublekinetics)
str(doublekinetics)


library(xtable)
xtable(head(doublekinetics))
xtable(tail(doublekinetics))


# Figura 1
library(lattice)
xyplot(cumulative ~ time , type = "b", groups = sample , xlab = "time [minute]", ylab = "cumulative %")


# Figura 3
par(mfrow = c(1,2))
V = 100
K = 60
x = 0:1000
y = V * x / (K + x)
plot(x, y, "l", ylim = c(0, 100), main = "cinetica di Michaelis e Menten")
abline(h = V, lty = 2)
points(K, V/2, pch = 19)
text(K, V/2, "(K, V/2)", pos = 4)
Y = x/y
a = K/V
b = 1/V
plot(x, Y, "l",  main = "linearizzazione", ylab = "y/x")
points(K, 2*K/V, pch = 19)
text(K, 2*K/V, "(K, 2K/V)", pos = 4)


linearized = time/cumulative


par(mfrow = c(1,2))
plot(time, cumulative, pch = 21, col = "black", bg = "green", main = "cumulative")
points(time[time < 121], cumulative[time < 121], pch = 21, col = "black", bg = "yellow")
plot(time, linearized, pch = 21, col = "black", bg = "green", main = "linearized")
points(time[time < 121], linearized[time < 121], pch = 21, col = "black", bg = "yellow")


linearizedS = scale(linearized)
timeS = scale(time)


sbagliato = lm(linearizedS ~ timeS * pH)
summary(sbagliato)
xtable(summary(sbagliato))


library(lme4)


mixed0 = lmer(linearizedS ~ timeS * pH + (timeS | sample))
summary(mixed0)


xtable(coef(summary(mixed0)))
xtable(as.data.frame(VarCorr(mixed0)))


mixed1 = lmer(linearizedS ~ timeS * pH + (1 | sample))
summary(mixed1)

xtable(coef(summary(mixed1)))
xtable(as.data.frame(VarCorr(mixed1)))


plot(timeS, linearizedS, pch = 21, col = "black", bg = "green", main = "linearized")
points(timeS[timeS < -0.12], linearizedS[timeS < -0.12], pch = 21, col = "black", bg = "yellow")


plot(timeS, linearizedS, pch = 21, col = "black", bg = "green", main = "linearized and rescaled")
abline(3.3667, 4.2311, col = "yellow", lwd = 4)
abline(3.3667-3.8533, 4.2311-3.5858, col = "green", lwd = 4)
points(timeS[timeS < -0.12], linearizedS[timeS < -0.12], pch = 21, col = "black", bg = "yellow")
points(timeS[timeS >= -0.12], linearizedS[timeS >= -0.12], pch = 21, col = "black", bg = "green")