Analisi di regressione non lineare (Venables, Ripley, 2002)

# se il modello non e' linearizzabile si deve ricorrere a tecniche iterative
# Per stimare Y ~ b0 + b1 * 2^(-X/theta) si usa nls
library(MASS)   # ci serve perche' contiene i dati wtloss
data(wtloss)

plot(wtloss)    # perdita di peso (Y) in funzione del tempo (X)
library(nls)    # questa e' la libreria che ci serve


model <- nls(Weight ~ b0 + b1 * 2^(-Days/th), data=wtloss, start=list(b0=90,b1=95,th=120), trace =TRUE)

# i parametri e l'espressione della funzione possono essere qualsiasi
# si devono specificare i valori iniziali
# si puo' anche specificare il gradiente della funzione (si veda VR, 2002)

model
summary(model)

Analisi della sopravvivenza (Dalgaard, 2002)

library(survival)

read.table("melanom.txt") -> melanom
names(melanom)
attach(melanom)

Surv(days,status==1)    # trsforma i dati in oggetti utili all'analisi di sopravvivenza
?Surv                   # sempre meglio dare un'occhiata all'help di Surv

survfit(Surv(days,status==1))

surv.all <- survfit(Surv(days,status==1)) # oggetti di tipo survfit
summary(surv.all)

plot(surv.all)      # disegno della curva di sopravvivenza

surv.bysex <- survfit(Surv(days,status==1) ~ sex)              # analisi per sesso

plot(surv.bysex)
plot(surv.bysex,conf.int=TRUE,col=c("blue","red"))

survdiff( Surv(days,status==1) ~ sex)   # e' un test del Chi-Quadrato

# stratificando anche per lo stato di avanzamento della malattia
# l'effetto della variabile sesso diminuisce

survdiff( Surv(days,status==1) ~ sex + strata(ulc))

coxph(Surv(days,status==1) ~ sex)

summary(coxph(Surv(days,status==1) ~ sex))

Analisi in componenti principali

read.table("nci.txt") -> nci
nci

cor(nci)
library(mva)

princomp(nci,cor=TRUE) -> pc

summary(pc)
plot(pc)      # Comp1 <-> nocivita', Comp2 <-> filtro/senza filtro

loadings(pc) # autovettori delle componenti principali
pc$scores    # coordinate delle osservazioni sulle componenti principali
str(pc)

Analisi fattoriale

# sopra abbiamo sospettato di due tipi di gruppi di variabili, chiediamo quindi 2 fattori
factanal(nci,factors=2) 

Analisi dei gruppi

hc <- hclust(dist(nci))      # costruisce i gruppi a partire da una matrice di distanze
plot(hc)                     # disegna il dendrogramma
cutree(hc,3)       # taglia il dendrogramma in modo da ottenere 3 gruppi
cutree(hc,4)

rect.hclust(hc,k=3,border="blue")
rect.hclust(hc,k=4)

hclust(dist(nci),method="centroid") -> hc1	      
hclust(dist(nci),method="single") -> hc2	      
par(mfrow=c(1,3))
plot(hc)
plot(hc1)
plot(hc2)
par(mfrow=c(1,1))

Analisi discriminante

data(iris3)
Iris <- data.frame(rbind(iris3[,,1], iris3[,,2], iris3[,,3]),
                        Sp = rep(c("s","c","v"), rep(50,3)))
table(Iris$Sp)
train <- sample(1:150, 75)
table(Iris$Sp[train])
z <- lda(Sp ~ ., Iris, prior = c(1,1,1)/3, subset = train)
plot(z)
predict(z, Iris[-train, ])$class

R, TeX, HTML, PDF e PostScript

Migliorare i grafici (Iacus-Masarotto, 2003)

load("dati1.rda")
plot(table(dati$Z), main = "lwd=1")
plot(table(dati$Z),lwd=10, main = "lwd=10")

str(colors())
palette()

plot(1,type="n")
abline(h=0.6, lty=2)
abline(h=0.8, lty=3)
abline(h=0.9, lty=4)
abline(h=1.2, lty=5)
abline(h=1.4, lty=6)

par(mfrow=c(4,2))
plot(dati$Z)
plot(dati$Z,type="l",main="type=\"l\"")
plot(dati$Z,type="b",main="type=\"b\"")
plot(dati$Z,type="c",main="type=\"c\"")
plot(dati$Z,type="o",main="type=\"o\"")
plot(dati$Z,type="h",main="type=\"h\"")
plot(dati$Z,type="s",main="type=\"s\"")
plot(dati$Z,type="S",main="type=\"S\"")
par(mfrow=c(1,1))

x <- c(1,2,5,9,10,11)
y <- c(1,7,5,4,3,1)
plot(x,y, main="Titolo", sub="sottotitolo", ylab="valori di y")

plot(x,y, main="Titolo", sub="sottotitolo", axes=FALSE,  xlab="numeri ics")
axis(2)
axis(1,x,c("uno","due","cinque","nove", "dieci","undici"))

x <- c(1,2,5,9,10,11)
y <- c(1,7,5,4,3,1)
plot(x,y)
abline(lm(y~x),lty=3)
text(4,3,"La retta di regressione")
text(6,4,expression(y[i]==hat(beta)[0]+hat(beta)[1]*x))

curve(dnorm(x),-5,5,lty=3)
curve(dt(x,df=1),-5,5,add=TRUE)
legend(-4.5,0.3,legend=c("normale","t Student"),lty=c(1,3))

x <- c(1,2,5,9,10,11)
y <- c(1,7,5,4,3,1)
z <- c(2,4,4,3,2,3)
plot(x,y,type="n",ylab="y,z")
points(x,y,pch=5,cex=2)
points(x,z,pch=8,cex=2)
lines(x,y)
lines(x,z,lty=3)
legend(6,6.5,legend=c("uomini","donne"),pch=c(5,8), lty=c(1,3))

curve(dnorm(x),-3,3,axes=FALSE,ylab="", xlab="",ylim=c(0,.5), main="I retini")
axis(1,c(-3,-1,0,1,3),c("","-z",0,"z",""))
vals <- seq(-3,-1,length=100)
x <- c(-3, vals, -1, -3)
y <- c(0, dnorm(vals),0, 0)
polygon(x,y,density=20,angle=45)
vals <- seq(1,3,length=100)
x <- c(1, vals, 3, 1)
y <- c(miny, dnorm(vals),miny, miny)
polygon(x,y,density=20,angle=45)
abline(h=miny)
text(0,0.45,expression(Phi(-z)== 1-Phi(z)))
lines(c(0,0),c(0,dnorm(0)))

vals <- seq(-3,-1,length=100)
x <- c(-3, vals, -1, -3)
y <- c(0, dnorm(vals),0, 0)
polygon(x,y,density=20,angle=45)