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

scala = ordered(scala, levels = c("uno", "due", "tre", "quattro"))
scala

boxplot(scala ~ caratteristica,  varwidth = TRUE)
table(scala, caratteristica)
library(xtable)
xtable(table(scala, caratteristica))
xtable(table(caratteristica))

wilcox.test(scala ~ caratteristica) # errore di sintassi
wilcox.test(as.numeric(scala) ~ caratteristica)


modello = glm(caratteristica ~ as.numeric(scala), family = binomial)
summary(modello)
xtable(summary(modello))



xtable(table(caratteristica, scala))

chisq.test(table(scala, caratteristica))


library(coin)
independence_test(scala ~ caratteristica, 
					ytrafo = rank_trafo, 
					distribution = exact())

chisq_test(scala ~ caratteristica)



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

xtable(istologiaP)


# riferimento
# http://www.wiley.com/legacy/wileychi/pesarin/material.html

n = table(Y)
n
C = length(n)
C
contr = rep(1/n, n)
contr

# scambiamo (intelligentemente) i segni del secondo gruppo, 'presente'
contr[-c(1:n[1])] = - contr[-c(1:n[1])]
contr
round(contr, digits = 3)


B = 10000
T = array(0, dim = c((B+1), 1))

# matrice verticale di una colonna sola con B + 1 zeri



T[1] = X %*% contr
T[1]
mean(X[1:27]) - mean(X[28:46])
# T[1] adesso vale -1.43, ossia la diff delle medie assente - presente



for (bb in 2 : (B+1) ) {
# permutiamo casualmente il campione X
X.star = sample(X)
# calcoliamo la differenza delle medie
T[bb] = X.star %*% contr
}

T

hist(T)




t2p = function(T){
if(is.null(dim(T))){T = array(T, dim = c(length(T), 1) )}

oth = seq(1:length(dim(T)))[-1]

B = dim(T)[1]-1
p = dim(T)[2]
if(length(dim(T)) == 3){C = dim(T)[3]}

rango = function(x){
r = 1 - rank(x[-1], ties.method = "min")/B+1/B
return(c(mean(x[-1]>=x[1]),r))
}

P = apply(T, oth, rango)
return(P)
}


P = t2p(T)
P[1]