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


library(xtable)
xtable(head(mieloma))
xtable(tail(mieloma))


table(patient)
xtable(t(as.matrix(table(patient))))

table(time)
xtable(t(as.matrix(table(time))))


table(patient)[2]
xtable(t(as.matrix(table(patient))))


# Figura 1
library(lattice)
xyplot(paraprotein ~ time | chain, type = "p",
   groups = patient, xlab = "time [minute]",
   ylab = "paraproteins [ng/ml]")


max(time)
max(paraprotein)



# per ragioni di convergenza, riscaliamo le variabili


fives = time/5  # cinquine di minuti
response = paraprotein/150 # median paraproteins
table(fives)
summary(response)


table(subject)  #  da 1 a 9

fives[1:5]
response[1:5]



110/5    # 22


stimaA = numeric(22)
stimaB = numeric(22)
stimaC = numeric(22)


logistica = function(x, a, b, c)
{0 + (a-0)/(1+exp(-b*(x-c)))}


for(i in 1:22)
{
	y = response[c(5*i - 4, 5*i - 3, 5*i - 2, 5*i - 1, 5*i )]
	x = fives[c(5*i - 4, 5*i - 3, 5*i - 2, 5*i - 1, 5*i )]
	regressione = nls( y ~ logistica(x, a, b, c), start =  list(a = 8.5, b = 0.5, c = 11))
	stimaA[i] = coef(regressione)[[1]]
	stimaB[i] = coef(regressione)[[2]]
	stimaC[i] = coef(regressione)[[3]]
}


stimaA
stimaB
stimaC


y = response[1:5]
x = fives[1:5]
regressione = nls( y ~ logistica(x, a, b, c), start =  list(a = 8.5, b = 0.5, c = 11))
regressione


patient22 = patient[time == 0]
chain22 = chain[time == 0]


data.frame(patient22, chain22, stimaA, stimaB, stimaC )
xtable(data.frame(patient22, chain22, stimaA, stimaB, stimaC ))

par(mfrow = c(1,2))
boxplot(stimaA ~ chain22, main = "plateau A")
boxplot(stimaB ~ chain22, main = "speed B")


moltomale = lm(stimaB ~ chain22 )
summary(moltomale)


library(lme4)
modelloB1 = lmer(stimaB ~ chain22 + (1| patient22))
summary(modelloB1)

modelloB0 = lmer(stimaB ~ 1 + (1| patient22))
ml.modelloB1 = lmer(stimaB ~ chain22 + (1| patient22), REML = FALSE)
ml.modelloB0 = lmer(stimaB ~   1     + (1| patient22), REML = FALSE)
anova(ml.modelloB1, ml.modelloB0)
# 0.004773 **


library(INLA)
formula = stimaB ~ chain22 + f(patient22, model = "iid")
output = inla( formula, family = "gaussian", data = data.frame(patient22, chain22, stimaA))
summary(output)