Magnetic fields have been shown to have an effect on living tissue as early as the 1930's. Plants have been shown to have an improved growth rate when raised in a magnetic field (Mericle et al., 1964). More recently, doctors and physical therapists have used either static or fluctuating magnetic fields to aid in pain management, most comonly for broken bones. In the case study presented here, Carlos Vallbona and his colleagues sought to answer the question "Can the chronic pain experienced by postpolio patients be relieved by magnetic fields applied directly over an identified pain trigger point?" Mericle, R. P. et al. "Plant Growth Responses". in : Biological effects of magnetic fields. New York: Plenium Press 1964. p183-195. Vallbona, Carlos et. al., "Response of pain to static Magnetic fields in postpolio patients, a double blind Pilot study" Archives of Physical Medicine and Rehabilitation. Vol 78, American congress of rehabilitaion medicine, p 1200-1203. Pain_Relief = read.csv("C:/Users/eleni/Desktop/Pain_Relief.txt", sep="") View(Pain_Relief) attach(Pain_Relief) names(Pain_Relief) stats.d(Pain_Relief$Score_1)
########################################################
########################################################
my.packages=c("TeachingDemos","MASS","pwr","car","lmtest","TSA","UsingR","SMIR","bootstrap", "VGAM","class","e1071","lattice") install.packages(my.packages,repos='http://cloud.r-project.org')
stats.d=function(x,s.trim=0.05,my.adjust=2,my.factor=2,...) {# descriptive statistics
options(warn=-1)## no warnings
xname=paste(deparse(substitute(x),500),collapse="\n")
#check if data set name has been missed
if(missing(x)) stop("ΠΡΟΣΟΧΗ! Δεν ορίσατε διάνυσμα δεδομένων")
#check the data set to be numeric
if(!is.numeric(x)) stop("ΠΡΟΣΟΧΗ! Τα δεδομένα δεν είναι αριθμητικό διάνυσμα. Δοκιμάστε as.numeric()")
#check the sample size
n0=length(x)
#count the number of NAs
n=length(na.omit(x));NMs=n0-n
if(n < 2) stop("ΠΡΟΣΟΧΗ! Λιγότερες από δύο έγκυρες παρατηρήσεις")
aux=paste(" ΠΕΡΙΓΡΑΦΙΚΑ ΣΤΑΤΙΣΤΙΚΑ για ",xname);names(aux)=c("");print(aux,quote=F) ##υπολογίζει μέσο και διορθωμένες για μέσο παρατηρήσεις x=na.omit(x) mean.x=mean(x) c.x=x - mean.x ## υπολογίζει ροπές γύρω από μέσο x2=c.x^2; x2=sum(x2) x3=c.x^3; x3=sum(x3) x4=c.x^4; x4=sum(x4) ##υπολογίζει λειασμένο μέσο tmean.x=mean(x,trim=s.trim) ##υπολογίζει δειγματική διακύμανση και μέτρα μορφής s.var.x=var(x) sd.x=sqrt(s.var.x) se.mean.x=sd.x/sqrt(n) skewness=(x3/n)/(x2/n)^1.5; se.skewness <- sqrt(6/n) kurtosis=((x4/n)/(x2/n)^2)-3; se.kurtosis <- sqrt(24/n) se.s.var.x=sqrt(((kurtosis+2)s.var.x^2/n)) if(abs(mean.x/sd.x)==0){coef.var=NaN} else {coef.var=(1+1/(4n))*sd.x/mean.x} min.x=min(x); max.x=max(x) median.x=median(x) a=quantile(x,probs=c(0.25,0.75),na.rm=TRUE,names=FALSE) Lqua.x=a[1]; Uqua.x=a[2] range.x=max.x-min.x SW.prob=shapiro.test(x)[[2]]## έλεγχος Shapiro-Wilk
par(oma=c(0,0,2,0))
layout(matrix(c(1,2,3,4,5,5),nr=3,byrow=TRUE)) ## δημιουργεί 2X3 graphs
aa=hist(x,plot=FALSE,...)
x.axis.low=min(mean.x-3*sd.x,(min(x)-0.5*sd.x))
x.axis.upper=max(mean.x+3*sd.x,min(x)+0.5*sd.x)
my.xlim=c(x.axis.low,x.axis.upper)
my.height=aa$density
my.ylim=c(0,1.2*max(my.height,dnorm(mean.x,mean.x,sd.x)))
plot(aa,freq=FALSE,xlim=my.xlim,ylim=my.ylim,main="Ιστόγραμμα και Κανονική Καμπύλη",
xlab="Τιμές Δείγματος",ylab="Πυκνότητα",col="grey95");box(which="plot",lty ="solid")
aux.1=seq(x.axis.low,x.axis.upper,length.out=300)
lines(aux.1,dnorm(aux.1,mean.x,sd.x),lwd=1,lty=1)
my.jitter=jitter(x,my.factor); rug(my.jitter)
abline(h=0,lty=1)
plot(aa,freq=FALSE,xlim=my.xlim,ylim=my.ylim,main="Ιστόγραμμα και Εξομαλυσμένη Καμπύλη",
xlab="Τιμές Δείγματος",ylab="Πυκνότητα",col="grey95");box(which="plot",lty ="solid")
par(new=TRUE)
plot(density(x,adjust=my.adjust,...),type="l",xlim=my.xlim,ylim=my.ylim,ann=FALSE,bty="n")
rug(my.jitter)
abline(h=0,lty=1)
boxplot(x,pars=list(pch=8),main="Θηκόγραμμα",xlab="Τιμές Δείγματος",ylim=my.xlim,horizontal=TRUE,col="grey95")
s.x <- c.x/sd.x
qqnorm(s.x,main="qq-plot για Τυπική Κανονική",xlab="Θεωρητικά Ποσοστημόρια",ylab="Δειγματικά Ποσοστημόρια",pch=16,col="grey50")
qqline(s.x)
plot(x,main="Δείκτης - Τιμή Δείγματος",xlab="Δείκτης",ylab="Τιμές Δείγματος",pch=16,col="grey50")
lines(c(-1,n+2),c(mean.x,mean.x),type="l",lty=2)
title(main=paste("ΓΡΑΦΙΚΑ για",xname),line=-0.6,out=TRUE,cex.main=2.0,lwd=0.8)
par(mfrow=c(1,1),oma=c(0,0,0,0))## επαναφέρει σε μονό διάγραμμα
##αποτελέσματα aux.name=paste(s.trim," λειασμένος μέσος") Όνομα=c("αριθ. παρατηρήσεων","αριθ. τιμών που λείπουν", "μέσος",aux.name,"διακύμανση","συν.μεταβλητότητας", "ελάχιστο","μέγιστο","διάμεσος","Q1","Q3","εύρος","ασυμμετρία","κύρτωση(διορθωμένη)","ελ-χος Shapiro-Wilk(p-τιμή)") τιμή=c(n0,NMs,mean.x,tmean.x,s.var.x,coef.var,min.x,max.x,median.x,Lqua.x,Uqua.x,range.x,skewness,kurtosis,SW.prob) τιμή=round(τιμή,digits=6) τυπ.σφάλμα=c(NA,NA,se.mean.x,NA,se.s.var.x,NA,NA,NA,NA,NA,NA,NA,se.skewness,se.kurtosis,NA) τυπ.σφάλμα=round(τυπ.σφάλμα,digits=6) d.names.1=rep("",length(τιμή)) d.names.2=c(" ΟΝΟΜΑ"," ΤΙΜΗ","ΤΥΠ.ΣΦΑΛΜΑ") d.names=list(d.names.1,d.names.2) aux=cbind(Όνομα,τιμή,τυπ.σφάλμα) dimnames(aux)=d.names
if(abs(mean.x/se.mean.x)>2){print(aux,na.print="",quote=FALSE)} else {print(aux,na.print="",quote=FALSE) aux1='***ΠΡΟΣΟΧΗ: ο μέσος είναι σχεδόν μηδέν, η τιμή του συν.μεταβλητότητας δεν είναι έγκυρη' names(aux1)=c("") print(aux1,quote=FALSE)} }
bc.t=function(x,l,inv=FALSE,...) {# bc.t performs a Box-Cox trans on the vector x with parameter l
my.epsilon=.Machine$double.neg.eps1e2 ifelse((x<=0 & inv==FALSE),return('STOP, non-positive values'),{ bc=function(x,l){if(abs(l)<my.epsilon) x=log(x) else x=(x^l-1)/l} bc.i=function(x,l){if(abs(l)<my.epsilon) x=exp(x) else x=(lx+1)^(1/l)} if(inv==FALSE) bc(x,l) else bc.i(x,l)})}
my.distr.plot.2=function(x,factor=NULL,my.amount=NULL,...) {
xname=paste(deparse(substitute(x),500),collapse="\n") if(is.null(factor)){aux.0=1;aux.2=list(x=x)} else {if(is.factor(factor)) {aux.0=length(levels(factor));aux.2=split(x,factor)} else stop("factor should be either NULL or a factor")}
G.mean.x=mean(x,na.rm=TRUE)
my.fu.1=function(x) {mean.x=mean(x,na.rm=TRUE) sd.x=sd(x,na.rm=TRUE) max.h=dnorm(mean.x,mean.x,sd.x) c(mean.x,sd.x,max.h)}
my.fu.2=function(v) {mean.x=v[1];sd.x=v[2] x.axis.low=mean.x-3sd.x x.axis.upper=mean.x+3sd.x x.axis=seq(x.axis.low,x.axis.upper,length=100) y.axis=dnorm(x.axis,mean.x,sd.x) lines(x.axis,y.axis,lwd=2,lty=2,col="black") x1=c(mean.x,mean.x);y1=c(0,dnorm(mean.x,mean.x,sd.x)) points(x1,y1,type='h',lwd=3,lty=1,col="black")}
res.s=matrix(nrow=aux.0,ncol=3) for(i in 1:aux.0){y=aux.2[[i]];res.s[i,1:3]=my.fu.1(y)} aa=hist(x,plot=FALSE,...) my.height=aa$density max.y=max(c(my.height,res.s[,3])) min.x=min(res.s[,1]-3res.s[,2]) max.x=max(res.s[,1]+3res.s[,2])
my.ylim=c(0,1.2max.y) my.xlim=c(min.x-0.15min.x,max.x+0.15*max.x) plot(aa,freq=FALSE,ylim=my.ylim,xlim=my.xlim,main="", col="bisque",xlab=xname,ylab="πυκνότητα",...) rug(jitter(x,amount=my.amount)) for(i in 1:aux.0){v=res.s[i,1:3];my.fu.2(v)}
points(G.mean.x,0,lwd=2,pch=25,bg="red") text(x=G.mean.x,y=0.05*max.y,labels="Grand Mean",col="red")
if(aux.0==1) title(main=paste("Κατανομή ",xname),line=-1,out=TRUE,cex.main=1.5,lwd=1.2) else title(main=paste("Κατανομή και μέσοι για ",xname," και ομάδες"),line=-1,out=TRUE,cex.main=1.5,lwd=1.2)}
my.distr.plot.3=function(x,factor=NULL,my.amount=NULL,my.adjust=2,...) {
xname=paste(deparse(substitute(x),500),collapse="\n") if(is.null(factor)){k=1;aux.1=list(x=x)} else {if(is.factor(factor)) {k=length(levels(factor));aux.1=split(x,factor)} else stop("factor should be either NULL or a factor object")}
options(warn=-1)
x.matrix=matrix(nrow=k,ncol=512) y.matrix=matrix(nrow=k,ncol=512) min.x=numeric(k) max.x=numeric(k) max.y=numeric(k)
for(i in 1:k){aux.temp=density(na.omit(aux.1[[i]]),adjust=my.adjust,...) x.matrix[i,]=aux.temp$x y.matrix[i,]=aux.temp$y min.x[i]=min(aux.temp$x) max.x[i]=max(aux.temp$x) max.y[i]=max(aux.temp$y)}
G.mean.x=mean(x,na.rm=TRUE) mean.x=sapply(aux.1,mean,na.rm=TRUE)
hist.plot=hist(x,plot=FALSE,...) my.height=max(hist.plot$density)
my.ylim=c(0,1.1max(max.y,my.height)) min.x=min(na.omit(x),min.x) min.x=min.x-0.1abs(min.x) max.x=max(na.omit(x),max.x) max.x=max.x-0.1*abs(max.x) my.xlim=c(min.x,max.x)
plot(hist.plot,freq=FALSE,xlim=my.xlim,ylim=my.ylim,main="",xlab=xname,ylab="πυκνότητα",col="bisque") rug(jitter(x,amount=my.amount))
for (i in (1:k)){ lines(x.matrix[i,],y.matrix[i,],col=grey(i0.55/k),lwd=4,...) x1=c(mean.x[i],mean.x[i]) aux.x=min(abs(x.matrix[i,]-mean.x[i])) index=which(abs(x.matrix[i,]-mean.x[i])==aux.x) y1=c(0,y.matrix[i,index[1]]) points(x1,y1,type="h",col=grey(i0.55/k),lwd=3,lty=2)}
points(G.mean.x,0,lwd=2,pch=25,bg="red") text(x=G.mean.x,y=0.05*max(max.y,my.height),labels="Grand Mean",col="red")
if(k==1) title(main=paste("Κατανομή ",xname),line=-1,out=TRUE,cex.main=1.5,lwd=1.2) else title(main=paste("Κατανομή και μέσοι για ",xname," και ομάδες"), line=-1,out=TRUE,cex.main=1.5,lwd=1.2) }
my.plot.pre=function(x,y,res.reg,res.pred,new.data,tr.function=function(x){x},...){
xname=paste(deparse(substitute(x),500),collapse="\n") yname=paste(deparse(substitute(y),500),collapse="\n")
aux1=as.numeric(new.data[,1]) aux2=as.numeric(tr.function(res.pred,...)[,1]) aux3=as.numeric(tr.function(res.pred,...)[,2:3]) plot(aux1,aux2,xlab=xname,ylab=yname, pch=4,col="red",lwd=2,xlim=c(min(x,aux1,na.rm=TRUE),max(x,aux1,na.rm=TRUE)), ylim=c(min(y,aux2,aux3,na.rm=TRUE),max(y,aux2,aux3,na.rm=TRUE))) text(aux1,aux2,labels="P",pos=4,offset=0.6,cex =1.2,col="red") title(main=paste("Πρόβλεψη της ",yname," με την ",xname),line=-1,out=TRUE,cex.main=1.5,lwd=1.2) title(sub="Προβλέψεις και 95% Διαστήματα για μια Ειδική Τιμή") points(x,y)
b.0=as.numeric(res.reg[[1]])[1] b.1=as.numeric(res.reg[[1]])[2] data.aug=c(x,aux1) data.aug.min=min(data.aug,na.rm=TRUE);data.aug.max=max(data.aug,na.rm=TRUE) x.aug=seq(data.aug.min,data.aug.max,length.out=300) y=b.0+b.1*x.aug lines(x.aug,tr.function(y,...),lty=1)
for(i in 1:length(aux1)) {x.axis=rep(aux1[i],2);y=res.pred;y.axis=as.numeric(tr.function(y,...)[i,2:3]) lines(x.axis,y.axis,lwd=2,col="red")}}
plot.cum.1=function(x,y=NULL,ylab1="πιθανότητα",my.main=xname.1,d=0.1,...) {# plot.cum.1 plots a cumulative density plot
xname=paste(deparse(substitute(x),500),collapse="\n") xname.1=paste("Αθροιστική Συνάρτηση Κατανομής της ",xname)
x=na.omit(x)
if(is.null(y)){x=as.numeric(names(table(x)));y=as.numeric(table(x))} else {O.I=order(x) x=x[O.I];y=y[O.I]}
if(is.null(y)==FALSE & min(y)<=0) stop("Αρνητικές ή Μηδενικές Συχνότητες στο x")
freq.x=cumsum(y)/sum(y) max.x=max(x);min.x=min(x) range.x=max.x-min.x aux.1=d*range.x low.x=min.x-aux.1;high.x=max.x+aux.1
x.axis=c(low.x,x,high.x) y.axis=c(0,freq.x,1) plot(x.axis,y.axis,type="s",xlab=xname,ylab=ylab1,ylim=c(0,1),...)
points(x,freq.x,pch=19,...) title(main=my.main,line=-2,out=TRUE,...) }
plot.cum.2=function(x,y,delta=700,ylab1="c.d.f.",...) {#x,y are numerical vectors
max.x=max(x) min.x=min(x) n=length(x) d=numeric(n-1) points=numeric(n-1)
step=(max.x-min.x)/delta for(i in 2:n){d[i-1]=x[i]-x[i-1] points[i-1]=floor(d[i-1]/step)}
xx=matrix(ncol=max(points),nrow=n-1) yy=matrix(ncol=max(points),nrow=n-1) i=1 for(i in 1:(n-1)){j=1 while(j<=points[i]) {xx[i,j]=x[i]+(j-1)*step;yy[i,j]=y[i];j=j+1}}
k.max=sum(points) x.axis=numeric();y.axis=numeric()
back=50 for (k in 1:back){x.axis[k]=x[1]-step*(back-k);y.axis[k]=0}
for(i in 1:(n-1)){for (j in 1:points[i]){k=k+1;x.axis[k]=xx[i,j];y.axis[k]=yy[i,j]}}
forth=50 for (kk in 1:forth){k=k+1;x.axis[k]=x[n]+kk*step;y.axis[k]=y[n]}
#plot(x.axis,y.axis,type="l",xlab="x-values",ylab="cum probs") plot(x.axis,y.axis,pch=".",xlab="x-values",ylab=ylab1,...) points(x,y,pch=19)
lines(c(x[1],x[1]),c(0,y[1]),...) for (i in 2:(n)){lines(c(x[i],x[i]),c(y[i-1],y[i]),...)} }
colorie=function (x, y1, y2, ...) {###EXAMPLE of use
polygon( c(x, x[length(x):1]), c(y1, y2[length(y2):1]), ... )}
ts.dec.plot=function(res.dec,my.main=NULL,...) {## plots the results of decompose()
s.detd=res.dec$x-res.dec$trend my.frequency=frequency(res.dec$x) my.length=length(res.dec$x) season.names=season(res.dec$x)[1:my.frequency]
par(mfrow=c(2,2)) #Ραβδόγραμμα για εποχικούς δείκτες my.width=c(rep(1,my.frequency)) Heights=res.dec$figure barplot(Heights,my.width,space=0.2,col=grey(0.9),xlab="season", ylab="seasonal index", main="seasonal index", axis.lty=1,names.arg=season.names,cex.main=1,ylim=c(1.1min(Heights,0),1.1max(Heights,0))) abline(h=c(0,0),lty=1,col="black",lwd=1)
#boxplots for seasonal factor season.x=factor(rep(1:my.frequency,len=my.length),labels=season.names) boxplot(s.detd~season.x,col=grey(0.9),xlab="season",ylab="seasonal component", main="seasonal component variation",cex.main=1)
#boxplots for random boxplot(res.dec$random~season.x,xlab="season",ylab="random component", main="random component variation",col=grey(0.9),cex.main=1)
ts.plot(res.dec$x,res.dec$trend,lty=c(1,2),main="series and trend",xlab="time", ylab="series",gpars=list(cex.main=1))
par(mfrow=c(1,1))
title(my.main,outer=TRUE,line=-1,cex.main=1.2)}
show.bootstrap.results=function (obs.x, res.b, my.statistic, ...) { theta.o = my.statistic(obs.x) theta.b = mean(res.b$thetastar) bias.b = theta.b - theta.o se.b = sd(res.b$thetastar) min.b = min(res.b$thetastar) max.b = max(res.b$thetastar) names = c("obs.statistic", "boot.statistic", "bias.b", "s.error.b", "min.b", "max.b") values = c(theta.o, theta.b, bias.b, se.b, min.b, max.b) show.b = data.frame(matrix(values, ncol = 6)) colnames(show.b) = names rownames(show.b) = c("") aa = hist(res.b$thetastar, plot = FALSE) i = 2 while (theta.o > aa$breaks[i]) i = i + 1 h1 = aa$counts[i - 1] while (theta.b > aa$breaks[i]) i = i + 1 h2 = aa$counts[i - 1] plot(aa, main = "bootstrap replications", xlab = "bootstrap values, blue:observed stat, red:boot stat", ylab = "freq") lines(c(theta.o, theta.o), c(0, h1), col = "blue", lwd = 2, lty = 3) lines(c(theta.b, theta.b), c(0, h2), col = "red", lwd = 2, lty = 6) show.b }
psRsqr=function(y,my.fit,k=2,...) {### ψευδο-R^2, για μη-γραμμικά μοντέλα
e=y-my.fit s2.a=sum(e^2)/(length(e)-k) SSRes=sum(e^2) SST=sum((y-mean(my.y))^2) Q_R.2=1-(SSRes/SST) res=c(sqrt(s2.a),Q_R.2) names(res)=c('Residual standard error:','pseudo R-squared:') res}
my.panel.cor=function(x,y,digits=3,prefix="",method="spearman",cex.cor,...){
##pairs HELP, MODIFIED usr=par("usr");on.exit(par(usr)) par(usr=c(0,1,0,1)) r=cor(x,y,method=method,use ="pairwise.complete.obs") txt=format(c(r,0.123456789),digits=digits)[1] txt=paste(prefix,"ρ=",txt,sep="") if(missing(cex.cor)) {if(abs(r)<=0.5)cex.cor=3.0 else cex.cor=(2.5*(abs(r)-0.5)+3)} if(abs(r)<=0.5) aux.1=2.0 else aux.1=10*(abs(r)-0.5)+2.0 text(0.5,0.5,txt,cex=cex.cor,col=grey(sqrt(1/aux.1)))}
my.panel.line=function (x, y, col = par("col"), bg = NA, pch = par("pch"), cex=1,col.smooth ="black", span = 2/3, iter = 3,lty.smooth=2,col.line="black",...) { points(x, y, pch = pch, col = col, bg = bg, cex = cex) my.lm=lm(y~x,...) abline(my.lm,col=col.line,...) ok <- is.finite(x) & is.finite(y) lines(lowess(x[ok],y[ok],f=span,iter=iter),col=col.smooth,lty=lty.smooth,...) }
##pairs HELP
panel.hist <- function(x,...) { usr <- par("usr"); on.exit(par(usr)) par(usr = c(usr[1:2], 0, 1.5) ) h <- hist(x, plot = FALSE) breaks <- h$breaks; nB <- length(breaks) y <- h$counts; y <- y/max(y) rect(breaks[-nB], 0, breaks[-1], y, col = "cyan", ...) rug(jitter(x))}
par(mar=c(4,5,4,5)) qqnorm(Pain_Relief$Score_1) stats.d(Pain_Relief$Score_1) t.test(Pain_Relief$Score_1,alternative="two.sided",mu=9,conf.level=0.95) #το 95% διάστημα εμπιστοσύνης προκύπτει 9.356961-9.803039 t.test(Pain_Relief$Score_2,alternative="two.sided",mu=9,conf.level=0.95) #το 95% διάστημα εμπιστοσύνης για τον μέσο του δείγματος μετά τη θεραπεία 5.131603 7.028397 t.test(Pain_Relief$Change,alternative = "two.sided",mu=9,conf.level = 0.95) #το 95% διάστημαεμπιστοσύνης για την διαφορά του πριν και μετά 2.545695 4.454305
plot(Pain_Relief,Pain_Relief$Score_1) boxplot(Pain_Relief$Active,Pain_Relief$Score_1)
x=Pain_Relief$Score_1 tr.x=log(x);tr.x=bc.t(x,0)## με την bc.t() από MY.LECTURE.FUNCTIONS stats.d(tr.x)
tr.x1=bc.t(x,0.2) stats.d(tr.x1)
tr.x2=bc.t(x,0.25) stats.d(tr.x2)
library(MASS)
boxcox(x~1)## library(MASS) title('Box-Cox μετασχηματισμός για Scores1')
aux=boxcox(x~1,lambda=seq(0,0.4,1/100),plotit=T) aux$x[which(aux$y==max(aux$y))]
tr.x3=bc.t(x,0.4) stats.d(tr.x3) #απορρίπτεται πάλι η κανονικότητα με την επιλογή της τιμής λ=0.4
#δημιουργία νεας μεταβλητής για ν
Elite =rep ("No",nrow(Pain_Relief )) Elite[Pain_Relief$Score_1<8]="Yes" Elite=as.factor(Elite) Elite<-factor(Elite,labels=c(0,1))#κωδικοποιούμε την μεταβλητή ελιτ με 0 για όσους είαι >8 και 1 για όσους το score1<8 pain1=data.frame(Pain_Relief,Elite)#προσθέτουμε τη νεα μεταβλητή στο data frame
#mu.0=bc.t(25,0.4);mu.0 #t.test(tr.x,alternative="two.sided",mu=mu.0,conf.level=0.95) #αυτο δεν το καταλαβα
power.t.test(n=NULL,delta=0.15,sd=sd(tr.x3,na.rm=TRUE),sig.level=0.05, power=0.80,type=c("one.sample"),alternative=c("two.sided")) #Το μέγεθος του δείγματος που απαιτείται είναι 18 άτομα
power.t.test(n=length(na.omit(Pain_Relief$Score_1)),delta=NULL,sd=sd(tr.x,na.rm=TRUE),sig.level=0.05, power=0.80,type=c("one.sample"),alternative=c("two.sided"))
my.delta=seq(0,2.5,length.out=30) aux=power.t.test(n=length(na.omit(x)),delta=my.delta,sd=sd(na.omit(x)),sig.level=0.05, power=NULL,type=c("one.sample"),alternative=c("two.sided")) plot(my.delta,aux$power,type='l',main="Καμπύλη Ισχύος",xlab="δέλτα",ylab="ισχύς")
#_________________________________________________________________________________________________________________
no.sucs=sum(ifelse(na.omit(x)>25,1,0)) no.trials=length(na.omit(x)) binom.test(no.sucs,no.trials) #απορρίπτεται η μηδενική υπόθεση άρα η κατανονή όχι συμμετρική
#_________________________________________________________________________________________________________________
library(bootstrap) my.mean=function(x){mean(x,na.rm=TRUE)} set.seed(100) boot.res=bootstrap(Pain_Relief$Score_1,2000,my.mean)
show.bootstrap.results(Pain_Relief$Score_1,boot.res,my.mean)## από MY.LECTURE.FUNCTIONS hist(boot.res$thetastar)
stats.d(boot.res$thetastar)
set.seed(100) res.b=bcanon(Pain_Relief$Score_1,nboot=2000,my.mean,alpha=c(0.025,0.05,0.95,0.975)) res.b$confpoints #το 95% διάστημα εμπιστοσύνης boostrap διαφέρει από το προηγούμενο και είναι 9.34-9.72
my.mean=function(x){mean(x,na.rm=TRUE)} set.seed(100) boot.res=bootstrap(Pain_Relief$Score_2,2000,my.mean)
show.bootstrap.results(Pain_Relief$Score_2,boot.res,my.mean)## από MY.LECTURE.FUNCTIONS hist(boot.res$thetastar)
stats.d(boot.res$thetastar)
set.seed(100) res.b=bcanon(Pain_Relief$Score_2,nboot=2000,my.mean,alpha=c(0.025,0.05,0.95,0.975)) res.b$confpoints #_________________________________________________________________________________________________________________ library(gmodels) CrossTable(Pain_Relief$Active,Pain_Relief$Score_1)
Pain.active<-split(Pain_Relief$Score_1,Pain_Relief$Active,drop=False)#1os tropos να χωρίσουμε το Score1 σε Αctive και Inactive
Origin.F<-Pain_Relief$Active Origin.F<-factor(Origin.F,labels = c("Active","Inactive"))#2os tropos Scores1.Group<-split(x,Origin.F) Scores2.Group<-split(Pain_Relief$Score_2,Origin.F) stats.d(Scores1.Group[[1]])#από τη συνάρτηση προκύπτουν όλα τα περιγραφικά μέτρα για τους ασθενείς πρίν την ακτινοθεραπεία stats.d(Scores1.Group[[2]])#και για τους ασθενεις πριν τη θεραπεία της ομάδας ελέγχου stats.d(Scores2.Group[[1]]) stats.d(Scores2.Group[[2]])
summary(Scores1.Group) summary(Scores2.Group)
Active=rep("Active",nrow(Pain_Relief))#φτίαχνουμε μια καινουργια στήλη με την μεταβλητη active Active[Pain_Relief$Active==2]="ακτινοθεραπεία" Active
library('psych') describeBy(Pain_Relief$Score_1,Active) describeBy(Pain_Relief$Score_2,Active)
Change.group<-split(Pain_Relief$Change,Origin.F) stats.d(Change.group[[1]]) stats.d(Change.group[[2]])
library("ggplot2") par(mfrow=c(1,2)) qplot(Scores1.Group[[1]]) qplot(Scores1.Group[[2]])
stats.d(Change.group[[1]]) stats.d(Change.group[[2]])
hist(Scores2.Group[[1]]) hist(Scores2.Group[[2]]) par(mfrow=c(1,1))
#F test για την ισότητα των διακυμάνσεων στις υποομάδες του Score2 active & not var.test(Scores2.Group[[1]],Scores2.Group[[2]],alternative = "two.sided")
#####################################
wilcox.test(Change.group[[1]],Change.group[[2]],paired=FALSE)####### ελεγχος ισότητας των δυομέσων wilcox.test(Scores2.Group[[1]],Scores2.Group[[2]],paired=FALSE)
#############ζευγαρωτες παρατηρησεις (έλεγχος για την ισόττα τν δύο μέσων)###### wilcox.test(Scores1.Group[[1]] ,Scores2.Group[[1]],alternative="less",paired=TRUE) wilcox.test(Scores1.Group[[2]],Scores2.Group[[2]],alternative="less",paired=TRUE)
#__________________________________________________________________________________________________ #vizualization #__________________________________________________________________________________________________ mydatt = read.csv("C:/Users/eleni/Desktop/Pain_Relief.txt", sep="") stats.d(mydatt$Score_1) stats.d(mydatt$Score_2) stats.d(mydatt$Change) summary(mydatt) Sys.setlocale("LC_CTYPE","Greek") On_freq=table(mydatt$Active) windows(width=6,height=4,rescale="fixed") pie(On_freq,col=gray(seq(0.6,1.0,length=3)), main="Κυκλικό Διάγραμμα για Active",cex=0.9,cex.main=1.0) boxplot(mydatt$fev,col="red",main="thikoramma gia to s;ynolo tvn dedomenvn",ylab="litra ekpneomenoy aera")
qqnorm(mydatt$Score_1,main = "QQ-PLOT (Score_1)") qqline(mydatt$Score_1) qqnorm(mydatt$Score_2,main = "QQ-PLOT (Score_2)") qqline(mydatt$Score_2) qqnorm(mydatt$Change,main = "QQ-PLOT (Change)") qqline(mydatt$Change) boxplot(mydatt$Score_1,mydatt$Score_2,mydatt$Change,col="azure2", main="Θηκογράμματα",ylab="Επίπεδο Πόνου",names=c("Score_1","Score_2","Change")) hist(mydatt$Score_1,breaks=10,prob=T,col="azure2", main="Ιστόγραμμα Score_1", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες") lines(density(mydatt$Score_1)) hist(mydatt$Score_2,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Score_2", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες") lines(density(mydatt$Score_2)) hist(mydatt$Change,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Change", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες") lines(density(mydatt$Change))
###diaxvrismos
Change.Act2<-split(mydatt$Score_2,mydatt$Active)
Change.Act1<-split(mydatt$Score_1,mydatt$Active)
Change.Act<-split(mydatt$Change,mydatt$Active)
par(mfrow=c(2,2))
boxplot(Change.Act2,main="Θηκόγραμμα Score_2",names=c("Ακτινοθεραπεία", "Ομάδα Ελέγχου"),ylab="Επίπεδο Πόνου",col="seagreen4")
boxplot(Change.Act1,main="Θηκόγραμμα Score_1",names=c("Ακτινοθεραπεία", "Ομάδα Ελέγχου"),ylab="Επίπεδο Πόνου",col="seagreen4")
boxplot(Change.Act,main="Θηκόγραμμα Change",names=c("Ακτινοθεραπεία", "Ομάδα Ελέγχου"),ylab="Επίπεδο Πόνου",col="seagreen4")
par(mfrow=c(3,2),main="QQ-plots")
qqnorm(Change.Act$1
,main = "Change: Ακτινοθεραπεία")
qqline(Change.Act$1
)
qqnorm(Change.Act$2
,main = "Change: Ομάδα Ελέγχου")
qqline(Change.Act$2
)
qqnorm(Change.Act1$1
,main = "Score_1: Ακτινοθεραπεία")
qqline(Change.Act1$1
)
qqnorm(Change.Act1$2
,main = "Score_1 Ομάδα Ελέγχου")
qqline(Change.Act1$2
)
qqnorm(Change.Act2$1
,main = "Score_2: Ομάδα Ελέγχου")
qqline(Change.Act2$1
)
qqnorm(Change.Act2$2
,main = "Score_2: Ομάδα Ελέγχου")
qqline(Change.Act2$2
)
par(mfrow=c(3,2))
hist(Change.Act1$1
,breaks=5,prob=T,col="azure2", main="Ιστόγραμμα Score_1: Ακτινοθεραπεία", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες")
lines(density(Change.Act1$1
))
hist(Change.Act1$2
,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Score_1 : Ομάδα Ελέγχου", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες")
lines(density(Change.Act1$2
))
hist(Change.Act2$1
,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Score_2: Ακτινοθεραπεία", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες")
lines(density(Change.Act2$1
))
hist(Change.Act2$2
,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Score_2: Ομάδα Ελέγχου", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες")
lines(density(Change.Act2$2
))
hist(Change.Act$1
,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Change: Ακτινοθεραπεία", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες")
lines(density(Change.Act2$1
))
hist(Change.Act$2
,breaks=10,col="azure2",prob=T, main="Ιστόγραμμα Change: Ομάδα Ελέγχου", xlab="Επίπεδο Πόνου",ylab = "Συχνότητες")
lines(density(Change.Act$2
))
stats.d(Change.Act1$1
)
stats.d(Change.Act1$2
)
stats.d(Change.Act2$1
)
stats.d(Change.Act2$2
)
stats.d(Change.Act$1
)
stats.d(Change.Act$2
)