r:notesmichel
Notes de Michel
Code
#coursR
#lancer Tinn-R
#sélectonner les commandes et effectuer par R avec Ctrl+D
#ou menu/R/send to R/selection
#programmé avec menu/R/hotkeys
#il faut parfois convertir une variable alpha en num et l'inverse
#a<-scan (taper enter pour chaque entree) introduit des num depuis l'écran
#b<-scan(what="") introduit des alpha
#b<-as.numeric(b) convertit en chiffres des chiffres introduits en scan(what="")
#b<-as.character(b) convertit en facteurs des chiffres introduits en num
##1 INTRODUCTION
getwd()
#donne le répertoire par défaut
setwd("g:/r")
#assigne un répertoire
ls()
#affiche les tables de .RData
load(".RData")
#Charge le fichier RData
str(lbw)
#affiche les 19 1eres lignes de lbw
rm(ndf)
#supprime ndf de RData. La suppression définitive se fait en sauvant l'image .RData
save.image()
loadhistory(file = ".Rhistory")
savehistory(file = ".Rhistory")
#ctrl L nettoie l'écran
## PACKAGES
library()
install.packages("foreign")
library(foreign)
remove.packages("foreign", lib=.Library)
library(help=foreign)
## TRIER
lbw <- lbw[order(lbw$PDMR) , ]
## EXECUTER UNE SYNTAXE
source("ndf.R")
sink(ndf)
#envoie la sortie vers le fichier ndf<<<;
sink()
#arrette la sortie
#on peut sauver l'historique des syntaxes avec menu/fichier/sauver l'historique
# fichier .rhistory
##2 DESCRIPTIVES
##2.1 STATISTIQUES
table(lbw$FUM)
#affiche les fréquences absolues
#On ne doit pas spécifier lbw si on attache le fichier
attach(lbw)
detach(lbw)
#detache le fichier lbw
table(lbw$FUM)/length(lbw$FUM)
#affiche les fréquences relatives
shapiro.test(lbw$PDBB)
#teste la normalité
summary(lbw$PDBB)
#decrit les variables (min p25 p50 p75 max moy)
quantile(lbw$PDBB,c(.33,.66))
#quantile
var(lbw$PDBB)
#variance
sd(lbw$PDBB)
#déviation standard
##2.2 GRAPHIQUES
barplot(table(lbw$FUM)/length(lbw$FUM) )
#ne pas oublier de fermer le graphe)
hist(lbw$PDBB,breaks=10)
boxplot(lbw]PDBB)
qqnorm(lbw$PDBB)
plot(density(lbw$PDBB))
##3 TRANSFORMATIONS
#Transformer avec value labels. variable label ne fonctionne pas même avec Hmisc.
#operateurs + - * / ^
#log exp sqrt trunc
#< > >= <= != (different de ) ==
# & (and) |(or)
lbw$lbw1<-9
lbw$lbw1[lbw$PDBB<2500] <-1
lbw$lbw1[lbw$PDBB>=2500] <-2
table(lbw$lbw1)
#fréquences absolues
table(lbw$lbw1)/length(lbw$lbw1)
#fréquences relatives
lbw$lbw1 <- NULL
#supprime la variable
#recode
library(car)
lbw$pdmrr <- recode(lbw$PDMR, " lo:110=1; 111:120=2; 121:130=3; 131:hi=4 ")
#value labels
lbw$pdmrr <- factor(lbw$pdmrr,levels = c(1,2,3,4),labels = c("<=110", "111-120" ,"121-130",">130"))
lbw$NSE <- factor(lbw$NSE,levels = c(1,2,3),labels = c("high", "mid", "low"))
lbw$lbw1 <- factor(lbw$lbw1,levels = c(1,2),labels = c("<2500", ">=2500"))
#transforme pdmrr et NSE en factor
table(lbw$pdmrr)/length(lbw$pdmrr)
#valeurs manquantes
x=read.fwf(file="h:/r/test.txt",widths=c(4,4),col.names=c("id","v1"))
#affiche NA pour var num ou <NA> pour var alpha
on doit introduire des blancs
avec read delim, le blanc alpha n'est pas affiché <NA>
x<-read.delim(file="e:/r/mis.raw",header=F)
les missing avec x <-read.table(ndf.txt) peuvent etre NA(num) ou blancs(alpha " ")
x <-read.table("e:/r/x.txt")
#on transforme un 9 en NA avec x$id[x$id==9] <-NA ou x$v1[x$v1=="b"] <-NA
mean(x$id,na.rm=TRUE) ,na.rm=TRUE permet d'afficher la moyenne meme si valeurs manquantes
#x$V2[x$V2==""] <-NA remplace un blanc par un <NA>
#dates
## 4 COMPARER DES PROPORTIONS
attach(lbw)
table(FUM,lbw1)
#chiffres absolus
prop.table(table(FUM,lbw1),1)
#row
prop.table(table(lbw1,FUM),2)
#col
summary(table(lbw1,FUM))
#chi2
library(epitools)
round(epitab(FUM, lbw1, method = "oddsratio")$tab, 2)
#or et rr
#le 2 signifie 2 digits > la virgule
round(epitab(FUM, lbw1, method = "riskratio")$tab, 2)
#analyse stratifiée
table(FUM,lbw1,NSE)
#autre façon +complexe lapply(split(lbw,f=lbw[,"NSE"]), function(x) table(x[,"FUM"],x[,"lbw1"],))
#%colonnes
prop.table(table(FUM,lbw1,NSE),c(2,3))
#autre façon +complexe
lapply(split(lbw,f=lbw[,"NSE"]), function(x) prop.table(table(x[,"FUM"],x[,"lbw1"]),2))
#%rangees
prop.table(table(FUM,lbw1,NSE),c(1,3))
mantelhaen.test(FUM,lbw1,NSE)
#calculator
hep=rbind(c(17,35),c(8,53),c(22,491))
hep<-matrix(c(17,8,22,35,53,491) ,nrow=3)
chisq.test(hep)
## 5 comparaison de moyennes
tapply(PDBB, FUM, mean, na.rm=T)
var.test(ndf$PDBB~ndf$FUM)
bartlett.test(ndf$PDBB~ndf$FUM)
levene.test(ndf$PDBB~ndf$FUM) ne marche qu'avec library car et FUM factor
si variances inegales
>t.test(lbw$PDBB~lbw$FUM)
si variances egales
t.test(lbw$PDBB~lbw$FUM,var.equal=TRUE)
wilcox.test(PDBB,FUM)
si deux groupes apparies
t.test(PDBB,PDMR,paired=T) exemple non valide
si deux groupes apparies non parametriques
pairwise.wilcox.test(PDBB,FUM,paired=T)
#stratifie
avec aggregate
aggregate(PDBB, by=list(FUM,NSE),FUN=mean, na.rm=TRUE)
##anova
tapply(PDBB, NSE, mean, na.rm=T)
result=aov(lbw$PDBB~lbw$NSE)
result=aov(PDBB~NSE,data=lbw)
summary(result)
coef(result) donne les contrastes
mean(ndv[ndv1==1]) donne la moyenne de ndv pour ndv=1, permet de calculer les moyennes
TukeyHSD.aov(result)
#donne contrastes
kruskal.test(PDBB,NSE)
## REGRESSION CORRELATION
cor(PDMR,PDBB)
plot(PDMR, PDBB)
result=lm(PDBB~PDMR,data=lbw)
summary(result)
plot(result)
abline(result)
## REGRESSION MULTIPLE
result=lm(WGT~AGE+HGT,data=regmul)
result=lm(WGT~AGE+HGT+(AGE*HGT),data=regmul)
#? collinearité?
## LOGISTIQUE
log <- glm(LBW ~ factor(FUM), family=binomial,data=ndf)
summary(log)
#donne OR
exp(coef(log))
#-2LL
logLik(log)
#tester la difference entre modeles?
1-pchisq(diff,df)
log <- glm(LBW ~ factor(FUM)+factor(NSE)+factor(FUM)*factor(NSE)
, family=binomial)
#l'interaction se définit par factor(FUM)*factor(NSE)
#result#OR=0.49 (comme epi) -2ll=229,8 p=0,0276 (comme epi)
#logistique lbw1~NSE donne 2,33 (mid) et 1,89 (low) -2ll=229,7
#?goodness of fit?
## SURVIE
library(survival)
attach(leukemia)
res <- survfit(Surv(TIME, STATUS)~ TREAT, conf.type="none",data=leukemia)
summary(res)
plot(res, main="Kapklan Meier",xlab="Time", ylab="Survival Probability")
#test du logrank
survdiff(Surv(TIME, STATUS)~ TREAT,data=leukemia)
## COX
#une interaction se definit par +TREAT:GENDER et se teste par la différence de LR test entre le modèle avec et le modèle sans interaction.
#Une analyse par strate s'effectue avec +strata(var)
res <- coxph(Surv(TIME, STATUS)~TREAT, data=leukemia,method=c("breslow"))
summary(res)
#HR=4,52
res <- coxph(Surv(TIME, STATUS)~TREAT+LOGWBC, data=leukemia,method=c("breslow"))
summary(res)
#HR=3,65 et LR=43,4
res <- coxph(Surv(TIME, STATUS)~TREAT+LOGWBC+LOGWBC:TREAT, data=leukemia,method=c("breslow"))
summary(res)
#LR=43,8 La diff de 0,8 n'est pas significative p=0,37 et l'interaction non plus p=0,51
1-pchisq(0.8,1)
#les 3 l suivantes testent l'interaction de TREAT et de LOGWBC avec le temps ns p=0,98
res <- coxph(Surv(TIME, STATUS)~TREAT, data=leukemia,method=c("breslow"))
time.res<-cox.zph(res)
time.res
#test graphique de l'interaction avec le temps résidus de schoenfeld: la droite doit rester lineaire
plot(time.res[1])
abline(h=0, lty=1)
#1 donne le plot pour TREAT et 2 pour LOGWBC
#Quand pha n'est pas vérifié pour GENDER par ex, il faut ajuster pour le temps
#en utilisant strata
#Mais il faut avoir verifié que les HR pour TREAT ne different pas dans chaque strate de GENDER
#sur base des IC
## IMPORTER
# Lire un fichier dbf, stata ou spss
lbw<-read.dbf("f:/r/lbw.dbf")
ttp_w<-read.spss("e:/r/cours/ttp_w.sav")
#Attention, certaines variables de fichier spss sont en format ordinal ou nominal et ne seront pas lues par R en survie notamment.
# Lire un fichier ASCII
#ascii fixed !il faut un carriage return a la derniere ligne
lbw=read.fwf(file="h:/r/lbw.dat",widths=c(5,3,4,3,3,3,3 ,3,3,4),col.names=c("ID","AG","PDMR","NSE","FUM","TRAV","HTA","IU","CONS","PDBB"))
x=read.fwf(file="g:/r/cours/lbw.dat",widths=c(5,3,4,3,3,3,3,3,3,4),col.names=c("id","agm","pm","nse","fum","trav","hta","iu","cons","pbb"))
# Lire sur le WEB
hmohiv<-read.table("http://www.ats.ucla.edu/stat/R/examples/asa/hmohiv.csv", sep=",", header = TRUE)
#lire fichier delim par tab
x<-read.delim(file="g:/r/mis.raw",header=F)
#introduire comme excel data.entry(x=c(NA))
a<-scan()
a<-scan(what="")
## EXPORTER
write.dta(lbw, file="e:/r/cours/toto.dta")
write.table(lbw,file="e:/r/cours/toto.txt") il y a un header avec le nom des variables
#on le relit avec
x=read.table("h:/r/cours/toto.txt")
##add
x=read.fwf(file="f:/r/file1.txt",widths=c(3,2),col.names=c("id","v1"))
y=read.fwf(file="f:/r/file2.txt",widths=c(3,2),col.names=c("id","v1"))
res<-rbind(x,r)
res
##merge
x=read.fwf(file="f:/r/file1.txt",widths=c(3,2),col.names=c("id","v1"))
y=read.fwf(file="f:/r/file3.txt",widths=c(3,2),col.names=c("id","v2"))
res<-merge(x,y,by="id",all=T)
res<-merge(x,y,by.x="id",by.y="id",all=T)
res
all=T:affiche tout
all=F:affiche commun
all.x=T:affiche les records de x
all.Y=T:affiche les records de Y
r/notesmichel.txt · Last modified: 2017/02/01 09:34 by carl