#####################


www = "http://www.biostatisticaumg.it/dataset/percussive.txt"
percussive = read.table(www, header = TRUE)
attach(percussive)
time = time - 12
treatment = factor(treatment)
treatment1 = treatment[1:70]
levels(treatment)[1] = "hfpv"
levels(treatment)[2] = "control"
str(percussive)
head(percussive)
tail(percussive)


#########  la bella figura "Trellis"

library(lattice)
xyplot(pafi ~ time | treatment, type = "b", groups = subject,
                xlab = "time (hour)", ylab = expression(PaO[2] / FiO[2]))





#### non c'e' dubbio sulle Estimates, i dubbi sorgono sugli Std. Errors:

ancovasbagliatocross = lm(pafi ~ time * treatment)
summary(ancovasbagliatocross)
ancovasbagliatoplus = lm(pafi ~ time + treatment)
summary(ancovasbagliatoplus)
AIC(ancovasbagliatocross, ancovasbagliatoplus)

###  scegliamo ancovasbagliatocross?

ancovasbagliatocross



#### importiamo il pacchetto "lme4"


library(lme4)


effettifissi = pafi ~ time * treatment                    ### sbagliato

effettimisti1 = pafi ~ time * treatment + (time | subject)   ### novita' della notazione

modellomisto1 = lmer(effettimisti1)
summary(modellomisto1)


### Gli Effetti Fissi

####   confrontiamo gli standard error 
####   e confrontiamo i consuntivi t
####   e confrontiamo i p-value



summary(ancovasbagliatocross)$coef
summary(modellomisto1)$coef


### Gli Effetti Casuali

summary(ancovasbagliatocross)
summary(modellomisto1)






##### intermezzo di Massimino  sulla normale bivariata,
###  ovvero la questione del cammello e del dromedario

https://it.padlet.com/massimo_borelli/myqvutrl9ooe




www = "http://www.biostatisticaumg.it/dataset/raul.csv"
raul = read.csv(www, header = TRUE)
attach(raul)
head(raul)

Ldhuno Ldhtre   Esito

xsane = Ldhuno[Esito == "benigno"]
ysane= Ldhtre[Esito == "benigno"]
length(xsane)
length(ysane)


mx = mean(xsane)		# 28.88712
sx = sd(xsane)			# 2.99404
my = mean(ysane)		# 22.23852
sy = sd(ysane)			# 2.344704
r = cor(xsane, ysane)	#-0.215333

n=50000
# genero due N(0,1) indipendenti
n1=rnorm(n)
n2=rnorm(n)

# combino
n3=r*n1+sqrt(1-r^2)*n2

# allora
xx = mx+sx*n1
yy = my+sy*n3
# sono le simulazioni cercate 
par(mfrow=c(1,2))
plot(xsane, ysane, col = "chartreuse4")
plot(xx,yy, col = "chartreuse4")


xx = trunc(mx+sx*n1)
yy = trunc(my+sy*n3)
# sono le simulazioni cercate 
par(mfrow=c(1,2))
plot(xsane, ysane)
abline(lm(ysane ~ xsane))
plot(xx,yy)
abline(lm(yy ~ xx))


##### fine dell'intermezzo: la normale bivariata non e' due normali appiccicate assieme



#### Attività  6.1 : tocca a voi, vi interrogo



library(lme4)

www = "http://www.biostatisticaumg.it/dataset/frattale.csv"
frattale = read.csv(www, header = TRUE)
attach(frattale)
str(frattale)
frattale



ipotesi1 = BC ~ DBC - 1
ipotesi2 = BC ~ DBC - 1 + (DBC - 1 | StudyID)


modello1 = lm(ipotesi1)
summary(modello1)

modello2 = lmer(ipotesi2)
summary(modello2)


# Per favore spiegate:

# perché gli errori standard dei due modelli sono diversi, 
#  cosa rappresentano le deviazioni standard degli effetti casuali, 
# perché non appare la parola Corr nel summary
# con rnorm fate qualche simulazione di dati ulteriore





### 6.4.1    Criterio di Informazione di Akaike


effettimisti1 = pafi ~ time * treatment + (time | subject)
effettimisti2 = pafi ~ time + treatment + (time | subject)

modellomisto2 = lmer(effettimisti2)
summary(modellomisto2)

###  !!! sbagliato
AIC(modellomisto1)
AIC(modellomisto2)




### giusto!!!


modellomisto1ml = lmer(effettimisti1, REML = FALSE)
modellomisto2ml = lmer(effettimisti2, REML = FALSE)
AIC(modellomisto1ml, modellomisto2ml)


#### raffinatezza:



effettimisti3 = pafi ~ time * treatment + (1 | subject) + (time - 1 | subject)

modellomisto3 = lmer(effettimisti3)
summary(modellomisto3)


###############################################
ANOVA  "solo effetti fissi"
###############################################




www = "http://www.biostatisticaumg.it/dataset/fresher.csv"
fresher = read.csv( www, header = TRUE )
attach(fresher)
str(fresher)
levels(gym)


relation3 = weight ~ gym
anovamodel = aov(relation3)
summary(anovamodel)


model3 = lm(relation3)
AIC(linearmodel)


boxplot(relation3)


gymOS = gym
levels(gymOS)
levels(gymOS)[2] = "occasionalsporty"
levels(gymOS)[3] = "occasionalsporty"
levels(gymOS)

relation3bis = weight ~ gymOS
model3bis = lm(relation3bis)
AIC(model3bis, model3)






###############################################
ANOVA  mixed model
###############################################




www = "http://www.biostatisticaumg.it/dataset/densitometry.txt"
densitometry = read.table(www, header = TRUE)
attach(densitometry)
names(densitometry)
str(densitometry)
tail(densitometry)

# subject measure     region


table(subject)

intensity = sqrt(measure)



situazione  = intensity ~ region
boxplot(situazione)
tapply(intensity, region, mean)

boxplot(situazione, horizontal = TRUE)



sbagliato = lm(intensity ~ region)
summary(sbagliato)


rel1 = intensity ~ region + (region | subject)
mixed1 = lmer(rel1)
summary(mixed1)


regionStemCo = region
levels(regionStemCo)
levels(regionStemCo)[1]
levels(regionStemCo)[1] = "brainstemcortex"
levels(regionStemCo)
levels(regionStemCo)[3] = "brainstemcortex"
levels(regionStemCo)




rel2 = intensity ~ regionStemCo + (regionStemCo | subject )
mixed2 = lmer(rel2)
summary(mixed2)



rel0 = intensity ~ 1 + (1 | subject )
mixed0ML = lmer(rel0, REML = FALSE)
mixed2ML = lmer(rel2, REML = FALSE)
AIC(mixed0ML, mixed2ML)


MLmixed2 = lmer ( intensity ~ regionStemCo + (regionStemCo | subject ) , REML = FALSE)