###TEST NON PARAMETRICI  --  DISTRIBUTION FREE
###popolazione non normale: potenza del t.test

popolazione1=rexp(1000,1)
popolazione2=rexp(1000,2)
par(mfrow=c(2,1))
hist(popolazione1,xlim=c(0,6),freq=F, ylim=c(0,1))
xfit=seq(0,6,length=100)
yfit=dexp(xfit,1)
lines(xfit,yfit, col="red")

hist(popolazione2,xlim=c(0,6),freq=F,ylim=c(0,1.5))
xfit=seq(0,6,length=100)
yfit=dexp(xfit,2)
lines(xfit,yfit, col="red")


mean(popolazione1) ##mu1=1
mean(popolazione2)  #mu2=0.5

###confronto di potenza dei due test: t  e wilcoxon

n=6
pvalue.t=c()
pvalue.w=c()
for(i in 1:100){
camp1=sample(popolazione1,n)
camp2=sample(popolazione2,n)
pvalue.t[i]=t.test(camp1,camp2)$p.value
pvalue.w[i]=wilcox.test(camp1,camp2)$p.value
}

which(pvalue.t < 0.05)
which(pvalue.w < 0.05)


## Hollander & Wolfe (1973), 29f.
## Hamilton depression scale factor measurements in 9 patients with
##  mixed anxiety and depression, taken at the first (x) and second
##  (y) visit after initiation of a therapy (administration of a
##  tranquilizer).
x <- c(1.83,  0.50,  1.62,  2.48, 1.68, 1.88, 1.55, 3.06, 1.30)
y <- c(0.878, 0.647, 0.598, 2.05, 1.06, 1.29, 1.06, 3.14, 1.29)
wilcox.test(x, y, paired = TRUE, alternative = "greater")
                                            # approximation

## Hollander & Wolfe (1973), 69f.
## Permeability constants of the human chorioamnion (a placental
##  membrane) at term (x) and between 12 to 26 weeks gestational
##  age (y).  The alternative of interest is greater permeability
##  of the human chorioamnion for the term pregnancy.
x <- c(0.80, 0.83, 1.89, 1.04, 1.45, 1.38, 1.91, 1.64, 0.73, 1.46)
y <- c(1.15, 0.88, 0.90, 0.74, 1.21)
wilcox.test(x, y, alternative = "g")        # greater


####test di Wilcoxon (dei ranghi segnati) (senza pareggi)
###simula due campioni normali accoppiati
n=35
x=rexp(n,1)
y=rexp(n,2)
##calcola la differenza 
x-y
## ordina le differenze
w=sort(x-y)
w
sign(w)
abs(w)
##ranghi del val. ass. della differenza
rank(abs(w))
#ranghi associati alle differenze positive e a quelle negative
rpos=rank(abs(w))[w>0]
rneg=rank(abs(w))[w<0]
##somme dei ranghi
sum(rpos)
sum(rneg)
sum(rpos)+sum(rneg)
###somma dei ranghi di un campione di n elementi
sum(1:n)
n*(n+1)/2
sum(rank(w))
##statistica nel caso di n>30:  in R t e' la somma dei ranghi positivi
t=sum(rpos)
z=(t-n*(n+1)/4)/sqrt(n*(n+1)*(2*n+1)/24)
z
##Pr(Zz)
2*(1-pnorm(z))
##uso della funzione
wilcox.test(x,y, paired=TRUE, conf.int=TRUE)
t.test(x,y,paired=TRUE)

###caso di un campione: si testa un valore per la mediana
median(x)
wilcox.test(x, mu=0.62)

### comparazione di due campioni indipendenti

n1=25; n2=30
x1=rexp(n1,1)
x2=rexp(n2,2)
###unisce i dati
x=c(x1,x2)
###crea etichette
lab=c(rep(1,length(x1)), rep(2,length(x2)))
lab

##matrice dati-etichette su 2 colonne
df=cbind(x,lab)
df

###ordina matrice rispetto ai dati (colonna 1)
dford=df[order(x),]
dford
ranghi1=which(dford[,2]==1)
ranghi1
ranghi2=which(dford[,2]==2)
###somme dei ranghi
s1=sum(ranghi1)
s1
s2=sum(ranghi2)
s2
s1+s2
##controllo con valore teorico
n=n1+n2
n*(n+1)/2
sum(ranghi1)-n1*(n1+1)/2
###la statistica e' approssimativamente  Z = N(0,1)
z=(s1-n1*(n1+n2+1)/2)/sqrt(n1*n2*(n1+n2+1)/12)
z
###  Pr(Zz); calcolo del p-value
2*(1-pnorm(z))

###confronto con la funzione di R
wilcox.test(x1,x2)
t.test(x1,x2)


#####nota sull'uso di order()
order(x)
sort(x)
x[order(x)]
##################


###test non parametrico di comparazione tra piu' di due gruppi

###peso in Kg dei pioppi dopo il primo anno con due suoli (1 umido, 2 secco) e trattamento (controllo, fertilizzante, irrigazione fert+irrig)

ct.1=c(0.15,0.02,0.16,0.37,0.22)
fe.1=c(1.34,0.14,0.02,0.08,0.08)
ir.1=c(0.23,0.04,0.34,0.16,0.05)
fi.1=c(2.03,0.27,0.92,1.07,2.38)
cbind(ct.1,fe.1,ir.1,fi.1)


pesi.1=c(ct.1,fe.1,ir.1,fi.1)
trattamento.1=gl(4,5,20, labels=c("ctl", "fer", "irr", "fi"))
trattamento.1

summary(aov(pesi.1~trattamento.1))

kruskal.test(pesi.1~trattamento.1)


ct.2=c(0.60,1.11,0.07,0.07,0.44)
fe.2=c(1.16,0.93,0.30,0.59,0.17)
ir.2=c(0.65,0.08,0.62,0.01,0.03)
fi.2=c(0.22,2.13,2.33,1.74,0.12)
cbind(ct.2,fe.2,ir.2,fi.2)


pesi=c(ct.1,fe.1,ir.1,fi.1,ct.2,fe.2,ir.2,fi.2)
pesi
suolo=gl(2,20,40)
suolo
trattamento.2=gl(4,5,20, labels=c("ctl", "fer", "irr", "fi"))
trattamento.2
trattamento=rep(trattamento.1,2)
trattamento
data.frame(suolo,trattamento,pesi)
par(mfrow=c(1,1))
interaction.plot(suolo,trattamento,pesi)
interaction.plot(trattamento,suolo,pesi)

anova(lm(pesi~suolo))
anova(lm(pesi~trattamento))
anova(lm(pesi~suolo*trattamento))

###ricostruzione dei dati heart.rate dalla libreria ISwR
##### 9 patients with congestive heart failure (insufficienza cardiaca) before and shortly after administration of enalaprilate (0 before and 30,60,120 minutes after)
hr=c(96, 110,  89,  95, 128, 100,  72,  79, 100,  92, 106,  86,  78, 124,  98,  68,  75, 106,  86, 108,  85, 78, 118, 100,  67,  74, 104,  92, 114,  83,  83, 118,  94,  71,  74, 102)

subj=gl(9,1,36)
subj
time=gl(4,9,36,labels=c(0,30,60,120))
time
heart.rate=data.frame(subj,time,hr)
heart.rate
interaction.plot(time,subj,hr)
interaction.plot(ordered(time),subj,hr)

anova(lm(hr~subj+time))

anova(lm(hr~subj*time))###non eseguible perche' non vi sono osservazioni ripetute

## test non parametrico per analisi della varianza a due vie
### per una sola osservazione per cella
###compara rispetto a time
friedman.test(hr~time|subj,data=heart.rate)
###compara rispetto a subj
friedman.test(hr~subj|time,data=heart.rate)

########## POTENZA DEL TEST basato su ipotesi di normalita'
###livello di significativita' = probablita' di rigettare l'ipotesi quando e' vera
###alpha = prob(errore di 1a specie)
### potenza = 1- probabilita' di accettare l'ipotesi quando e' falsa
###beta = 1-prob(errore di 2a specie)
### delta = true diference in means; sd = standard deviation; sig.level=significance level; power
power.t.test(delta=0.5,sd=2,sig.level=0.05, power=0.9)
power.t.test(n=300, delta=0.5,sd=2,sig.level=0.01)
power.t.test(delta=0.5,sd=0.5,sig.level=0.01, power=0.9)

###dati accoppiati
power.t.test(delta=3,sd=2,sig.level=0.01, power=0.9,type="two.sample")
power.t.test(delta=3,sd=2,sig.level=0.01, power=0.9, type="paired")
power.t.test(n=20,delta=1)
power.t.test(power=0.90, delta=1,alternative="one.sided")
power.t.test(power=0.90, delta=1,alternative="two.sided")