Some useful R code

b<-rep(0,12)

xstar<-matrix(b,ncol=6)

library(MASS)

a<-c(1.1,1.2,1.3,1.4,1.5,1.6,2.1,2.2,2.3,2.4,2.5,2.6)

x<-matrix(a,ncol=6)

xm<-data.frame(x)

object<-lm.ridge(X1~.,xm,lambda=0.003)

coef4.1<-coef(object)

acoef<-abs(coef4.1)

for (i in 2:6)

{xstar[,i]<-x[,i]/acoef[i]

}

xstar

x[,1]

acoef

Iteration:

g1=function(x)

{a=x/(x+1)

return(a)

}

f1=function(g1,y,nlim)

{interno=0

repeat{

       interno=interno+1

         if (interno<nlim)

         {y=g1(y)

         }

         else

         {

          cat("interation=",interno,"\n")

          break

         }

       }

   return(y)

}

f1(g1,1,3)

a<-c(1,1,1,2,2,2,3,3,3)

a[-(4:6)]

1 1 1 3 3 3

g1=function(x,y)

{a=y/(x+y)

 return(a)

}

g2=function(x,y)

{b=y*(x+y)

return(b)

}

g3=function(x,y,nlim)

{iterno=0

repeat{

iterno=iterno+1

if (iterno>nlim)

{cat("iterno=",iterno,"\n")

break}

else

{x=g1(x,y)

y=g2(x,y)

}

}

return(list(x=x,y=y))

}

g3(1,1,2)

iterno= 3

$x

[1] 0.75

$y

[1] 3.375

x=1/2

y=1/2+1

x=(3/2)/(3/2+1/2)

y=(3/2)*(0.75+1.5)

y

3.375

a=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1)

length(a)

m=matrix(a,ncol=7)

write(t(m),file="C:/Users/Soemone/Desktop/08 Fall-/Research/S Xu/Simulate Large Data/a.txt",ncol=7)

An algorithm that can normalize the columns and rows of a matrix:

m<-matrix(c(1,3,5,4,7,4,6,1,5),3,3)

rowsum<-c(1,2,3)

colsum<-c(1,2,3)

while((max(rowsum)-min(rowsum))>0.0001 && (max(colsum)-min(colsum))>0.0001)

{ for (i in 1:3)

{

rowsum<-apply(m,1,sum)

m[i,]<-m[i,]*12/rowsum[i]

}

for (i in 1:3)

{

colsum<-apply(m,2,sum)

m[,i]<-m[,i]*12/colsum[i]

}

}

To get the source codes from R

In a package: use methods(package)

In a function: input the name of the function, then return, if the function is in some where, then getAnywhere(function)

fun1<-function(n)

{ if(n==1)

  list(success="The computation is successful!")

else

if (n<=0)

list(fail="We need a positive number!")

else while(n!=1)

{

if (n%%2==0) n<-(n/2)

   else n<-(n*3+1)

print(n)

}

}

fun1(3)

[1] 10

[1] 5

[1] 16

[1] 8

[1] 4

[1] 2

[1] 1

**

fun1<-function(n)

{ if(n==1)

  list(success="The computation is successful!")

else

if (n<=0)

list(fail="We need a positive number!")

else while(n!=1)

{

if (n%%2==0) n<-(n/2)

   else n<-(n*3+1)

}

print(n)

}

fun1(3)

> fun1(3)

[1] 1

How to export a data frame from R into SAS

Using the package “foreign”, and the function “write.foreign”

Suppose we have a data frame named “b” in R. If the data is not data frame, we can use “as.data.frame” to convert it into a data frame. We set up an empty file name a.txt, then we write a sas code to infile the data from a.txt. At last, we use the following codes to finish the export.

library(foreign)

a<-matrix(c(1,2,3,4,5,6,7,8,9),nrow=3)

b<-as.data.frame(a)

write.foreign(b,"C:\\Users\\Soemone\\Desktop\\a.txt","C:\\Users\\Soemone\\Desktop\\a.sas",package="SAS",dataname="one")

dataname is an optional statement, using that we can convert the name of the data in SAS into the one that we want such as “one” here. May try different kind of data file such as Excel file?

In fact, we do not need to set up the data file a.txt and the sas file a.sas. Just write a path for the data file and the sas codes file, but do not need to put anything there. Then run the codes above, it works.

R can not read Excel file xls directly unless the excel file is converted into the csv file.

Draw multiple graphs in graph.

R codes:

par(mfrow=c(2,2))

a=matrix(rep(0,16),nrow=4)

for (i in 1:4)

{plot(rnorm(a[i,]))

}

Results:

Another codes:

par(mfrow=c(2,3))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

par(mfrow=c(3,2))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

par(mfcol=c(2,3))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

par(mfcol=c(3,2))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

Outer:

x<-c(1,1,1,1,1)

y<-c(1,1,1,1,1)

pf<-function(x,y){x+y}

z<-outer(x,y,FUN="pf")

results:

 2    2    2    2    2

 2    2    2    2    2

 2    2    2    2    2

 2    2    2    2    2

 2    2    2    2    2

x<-c(1,1,1,1,1)

> y<-c(1,1,1,1,1)

> pf<-function(x,y){x+y}

> z<-outer(x,y,fun="pf")

> z

     [,1] [,2] [,3] [,4] [,5]

[1,]    1    1    1    1    1

[2,]    1    1    1    1    1

[3,]    1    1    1    1    1

[4,]    1    1    1    1    1

[5,]    1    1    1    1    1

> 

Transpose a data frame

Original data frame:

id year_mnth revenue

101 2009-01  1100

101 2009-02  1200

102 2009-02  1900

103 2009-03  2000

R codes:

dat<-read.table(“data”,header=T)

tapply(dat$revenue, list(dat$id, dat$year_mnth), sum)

result:

2009-01 2009-02 2009-03

101    1100    1200      NA

102      NA    1900      NA

103      NA      NA    2000

Graphical images for multi-variables

1. Outline:

outline<-function(x,txt=TRUE)

{if (is.data.frame(x)==TRUE)

  x<-ax.matrix(x)

  m=nrow(x);n=ncol(x)

  plot(c(1,n),c(min(x),max(x)))

for(i in 1:m)

  {lines(x[i,],col=i)

   if (txt==TRUE)

   {k<-dimnames(x)[[1]][i]

    text(1+(i-1)%%n,x[i,1+(i-1)%%n],k)

   }

  }

}

2. Star:

      

       stars()

     

3. harmonic curve:

harmonc<-function(x)

{ if (is.data.frame(x)==TRUE)

      x<-as.matrix(x)

   t<-seq(-pi,pi,pi/30)

   m<-nrow(x);n<-ncol(x)

   f<-array(0,c(m,length(t)))

   for (i in 1:m)

   {f[i,]<-x[i,1]/sqrt(2)

    for (j in 2:n)

    { if (j%%2==0)

      f[i,]<-f[i,]+x[i,j]*sin(j/2*t)

      else

      f[i,]<-f[i,]+x[i,j]*cos(j%/%2*t)

     }

    }

 plot(c(-pi,pi),c(min(f),max(f)))

 for (i in 1:m) lines(t,f[i,],col=i)

}

data:

         x1 x2 x3 x4 x5

1 96 95 98 97 100

2 91 86 83 84 85

3 92 93 94 97 100

4 90 87 88 89 95

5 82 95 94 93 97

6 78 79 76 74 73

7 75 78 76 73 72

codes:

x<-read.table(“C:\\Users\\Soemone\\Desktop\\d.txt”)

source("C:\\Users\\Soemone\\Desktop\\a.R")

harmonc(x)

graph:

copy and paste codes from word document:

When there is a “ which is originally written using word, we need to change it to ", otherwise, there will be an error:

Error: unexpected input in "x<-read.table(“"

Can give a scalar as a value to a vector:

codes:

a<-seq(-1,1,1/10)

f<-array(0,c(4,length(a)))

print(f)

for (i in 1:4)

  {f[i,]<-1

   print(f[i,])

   }

results:

a<-seq(-1,1,1/10)

> f<-array(0,c(4,length(a)))

> print(f)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]

[1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

[2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

[3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

[4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

     [,15] [,16] [,17] [,18] [,19] [,20] [,21]

[1,]     0     0     0     0     0     0     0

[2,]     0     0     0     0     0     0     0

[3,]     0     0     0     0     0     0     0

[4,]     0     0     0     0     0     0     0

>

> for (i in 1:4)

+   {f[i,]<-1

+    print(f[i,])

+    }

 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

or:

a<-seq(-1,1,1/10)

f<-array(0,c(4,length(a)))

print(f)

for (i in 1:4)

  {f[i,]<-sin(a)

   print(f[i,])

   }

results:

> a<-seq(-1,1,1/10)

> f<-array(0,c(4,length(a)))

> print(f)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]

[1,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

[2,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

[3,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

[4,]    0    0    0    0    0    0    0    0    0     0     0     0     0     0

     [,15] [,16] [,17] [,18] [,19] [,20] [,21]

[1,]     0     0     0     0     0     0     0

[2,]     0     0     0     0     0     0     0

[3,]     0     0     0     0     0     0     0

[4,]     0     0     0     0     0     0     0

>

> for (i in 1:4)

+   {f[i,]<-sin(a)

+    print(f[i,])

+    }

 [1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

 [7] -0.38941834 -0.29552021 -0.19866933 -0.09983342  0.00000000  0.09983342

[13]  0.19866933  0.29552021  0.38941834  0.47942554  0.56464247  0.64421769

[19]  0.71735609  0.78332691  0.84147098

 [1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

 [7] -0.38941834 -0.29552021 -0.19866933 -0.09983342  0.00000000  0.09983342

[13]  0.19866933  0.29552021  0.38941834  0.47942554  0.56464247  0.64421769

[19]  0.71735609  0.78332691  0.84147098

 [1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

 [7] -0.38941834 -0.29552021 -0.19866933 -0.09983342  0.00000000  0.09983342

[13]  0.19866933  0.29552021  0.38941834  0.47942554  0.56464247  0.64421769

[19]  0.71735609  0.78332691  0.84147098

 [1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

 [7] -0.38941834 -0.29552021 -0.19866933 -0.09983342  0.00000000  0.09983342

[13]  0.19866933  0.29552021  0.38941834  0.47942554  0.56464247  0.64421769

[19]  0.71735609  0.78332691  0.84147098

Gibbs sampling:

R codes from Bayesian Statistics: An introduction.

nr<-50

m<-5000

k<-10

n<-16

alpha<-2

beta<-4

lambda<-16

maxn<-24

h<-rep(0,n+1)

for (i in 1:m)

{ 

  pi<-runif(1)

  for (j in 1:k)

    {

     y<-rbinom(1,n,pi)

     newalpha<-y+alpha

     newbeta<-n-y+beta

     pi<-rbeta(1,newalpha,newbeta)

    }

  for (t in 0:n)

   {

     if (t==y)

        h[t+1]<-h[t+1]+1

    }

}

print(h)

barplot(h)

results:

> print(h)

 [1] 200 372 508 573 537 551 470 428 376 300 235 184 141  67  37  16   5

> barplot(h)

R codes translated from computational statistics handbook with MATLAB:

k=1000

m=500

a=2

b=4

x=rep(0,k)

y=rep(0,k)

n=16

x[1]=rbinom(1,n,0.5)

y[1]=rbeta(1,x[1]+a,n-x[1]+b)

for (i in 2:k)

{ x[i]=rbinom(1,n,y[i-1])

  y[i]=rbeta(1,x[i]+a,n-x[i]+b)

}

fhat=rep(0,17)

for (j in 1:17)

 { fhat[j]=mean(dbinom(j-1,n,y[500:k]))}

barplot(fhat)

Array a_1,...a_n, want to get sum( a_j: a_j>K) and sum_{1<=j<=n} {max(a_j,K)}

a <- c(3,6,2,5,7,9)
K <- 4
## Question 1
sum(a[a > K])
## Question 2
sum(apply(rbind(K, a), 2, max))

Or

##Question 1:

sum(t(a)%*%(a>k))

##Question 2:

sum(t(a)%*%(a>K))+sum(t(rep(K))%*%(a<K))

Multiple lines put in one plot in R:

a<-scan()

0.6566        2.1185        1.2320       

0.647         2.0865        1.2325       

0.6532        2.1060        1.2287       

0.6487        2.1290        1.2313       

0.6594        2.1285        1.2341       

0.6577        2.1070        1.2343       

0.6579        2.1345        1.2340       

0.6734        2.1705        1.2362       

0.675         2.1845        1.2372       

0.6592        2.1550        1.2340       

0.6647        2.1710        1.2305

      

a<-as.matrix(a)

dim(a)

a<-matrix(a,nrow=11)

dim(a)

colnames(a)<-c("A","B","C")

y<-1:11*5

y

matplot(y,a,type="l")

Notes: When you read in the multiple data, you should type a<-scan() first, then paste the data in the script, then first submit “a<-scan()”, then select one data line, then submit it. Do this line by line. At last, select a blank line, then submit it. In this way we can read in all the data.

Matplot and legend

R codes:

matplot(a,type="o",pch=c(1,2,3))

legend(1, 9, c("one","two","three"),

       pch =c(1,2,3), col =1:3)

results:

Another codes:

matplot(a,type="o",pch=c(1,2,3),col=5:7)

legend(1, 9, c("one","two","three"),

       pch =c(1,2,3), col =5:7)

Results:

Another codes:

matplot(a,type="o",pch=c(1,2,3),col=1:3)

legend(1, 9, c("one","two","three"),

       pch =c(1,2,3), col =5:7)

Without dots in the plot.

data:

b

1 2 3

4 5 6

7 8 1

2 3 .

codes:

a<-read.table("C:\\Users\\Soemone\\Desktop\\b.txt")

matplot(a,col=1:3,type="l")

legend(1,7,c("one","two","three"),lty=1:3,col=1:3)

results:

Another codes:

a<-read.table("C:\\Users\\Soemone\\Desktop\\b.txt")

matplot(a,col=4:6,type="l",lty=4:6)

legend(1,7,c("one","two","three"),lty=4:6,col=4:6)

results:

When “.” is used to indicate the missing value, there will be a warning message when we draw the plot. But if “NA” is used to indicate the missing value, there will be no warning message.

Another codes:

a<-read.table("C:\\Users\\Soemone\\Desktop\\b.txt")

matplot(a,col=4:6,type="l",lty=c(1,1,1),lwd=rep(1.8,3))

legend(1,7,c("one","two","three"),lty=c(1,1,1),col=4:6,lwd=rep(1.8,3))

results:

How to calculate a factorial in R:

prod(1:n)=gamma(n+1)

An R code in the paper.

a<-scan()

643.3944  595.0646  577.5084

589.3512  574.781   566.2284

588.3362  577.3106  560.6669

588.9926  578.6509  560.6045

587.5358  578.1017  560.3376

588.4713  578.6467  560.7271

589.0297  578.5529  560.7797

588.5336  578.4009  560.7044

589.0542  579.102   560.6257

588.28    578.4213  560.3738

588.2308  579.052   560.2325

c<-matrix(a,nrow=11,byrow=T)

b<-c(0.95,0.85,0.75,0.65,0.55,0.45,0.35,0.25,0.15,0.05,0)

length(b)

matplot(b,c,type="o",pch=c(1,2,3),col=rep(1,3),xlab="Level",ylab="PE")

legend(0.1,640, c("Sample size=350","Sample size=500","Sample size=650"),

       pch =c(1,2,3), col =c(1,1,1))

Program the indicator variables.

ff <- factor(sample(letters[1:5], 25, replace=TRUE))

ff

 [1] e b d e a e e a e d a c a e c e a b c e a e e e b

Levels: a b c d e

diag(nlevels(ff))[ff,]

      [,1] [,2] [,3] [,4] [,5]

 [1,]    0    0    0    0    1

 [2,]    0    1    0    0    0

 [3,]    0    0    0    1    0

 [4,]    0    0    0    0    1

 [5,]    1    0    0    0    0

 [6,]    0    0    0    0    1

 [7,]    0    0    0    0    1

 [8,]    1    0    0    0    0

 [9,]    0    0    0    0    1

[10,]    0    0    0    1    0

[11,]    1    0    0    0    0

[12,]    0    0    1    0    0

[13,]    1    0    0    0    0

[14,]    0    0    0    0    1

[15,]    0    0    1    0    0

[16,]    0    0    0    0    1

[17,]    1    0    0    0    0

[18,]    0    1    0    0    0

[19,]    0    0    1    0    0

[20,]    0    0    0    0    1

[21,]    1    0    0    0    0

[22,]    0    0    0    0    1

[23,]    0    0    0    0    1

[24,]    0    0    0    0    1

[25,]    0    1    0    0    0

>

> model.matrix(~ ff - 1)

   ffa ffb ffc ffd ffe

1    0   0   0   0   1

2    0   1   0   0   0

3    0   0   0   1   0

4    0   0   0   0   1

5    1   0   0   0   0

6    0   0   0   0   1

7    0   0   0   0   1

8    1   0   0   0   0

9    0   0   0   0   1

10   0   0   0   1   0

11   1   0   0   0   0

12   0   0   1   0   0

13   1   0   0   0   0

14   0   0   0   0   1

15   0   0   1   0   0

16   0   0   0   0   1

17   1   0   0   0   0

18   0   1   0   0   0

19   0   0   1   0   0

20   0   0   0   0   1

21   1   0   0   0   0

22   0   0   0   0   1

23   0   0   0   0   1

24   0   0   0   0   1

25   0   1   0   0   0

attr(,"assign")

[1] 1 1 1 1 1

attr(,"contrasts")

attr(,"contrasts")$ff

[1] "contr.treatment"

How to find the rank of a matrix?

qr()$rank

example:

 (a<-diag(c(1,1,1)))

     [,1] [,2] [,3]

[1,]    1    0    0

[2,]    0    1    0

[3,]    0    0    1

qr(a)$rank

[1] 3

If you wish to know the index positions of TRUE elements of a logical vector x, then use which(x).

example:

x<-c(1,1,1,2,3,3,3,4,4,4)

which(x%%2==0)

4  8  9 10

 x<-c(1,1,1,2,3,3,3,4,4,4)

 a=which(x%%2==0)

 (x[a])

[1] 2 4 4 4

Textbook: Introduction to scientific programming and simulation with R:

Ex 3.1

x.values <- seq(-2, 2, by = 0.1)

n <- length(x.values)

y.values <- rep(0, n)

for (i in 1:n) {

x <- x.values[i]

if (x<0||x==0)

{y=-x^3}

else if ((0<x && x<1)||x==1)

   {y=x^2}

   else {y=sqrt(x)}

y.values[i] <- y

}

plot(x.values,y.values,type="l")

Ex 3.2

af<-function(x,n)

{h=0

for (i in 0:n)

{h=h+x^i}

show(h)

}

af(2,3)

15

Ex 3.4

af<-function(x,n)

{h=0

i=0

while (i<=n)

{h=h+x^i

i=i+1}

show(h)

}

af(2,3)

15

Using vector operations

x=2; n=3

(sum(rep(x,(n+1))^(0:n)))

15

Ex. 3.5

a<-c(1,1)

b<-c(cos(1.5),-sin(1.5),sin(1.5),cos(1.5))

theta<-matrix(b,ncol=2,byrow=T)

(c<-theta%*%a)

           [,1]

[1,] -0.9267578

[2,]  1.0682322

Ex 3.6

a<-c(1.1,1.2,1.3,1.4,1.5,1.6,1.7)

h=1

for (i in 1:length(a))

{h=h*a[i]}

(ave=h^(1/length(a)))

[1] 1.385527

x<-c(1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8)

h=0

for ( i in 1:length(x))

{h=h+(1/x[i])

}

a=h^(-1)

print(a)

0.1766348

(b=(sum(x^(-1)))^(-1))

0.1766348

Ex 3.7

x<-c(1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8)

sum=0

for (i in 1:length(x))

{if (i%%3==0)

{sum=sum+x[i]}

}

print(sum)

2.9

Ex 3.10

x<-c(1.1,0.7,0.4,0.3,1.9,1.8,1.7,2.1)

x.min=x[1]

for (i in 2:length(x))

{if (x[i]<x.min||x[i]==x.min)

{x.min=x[i]}

}

print(x.min)

0.3

Ex 3.13

t<-seq(0,10,0.001)

 r<-sqrt(t)

 theta<-2*pi*t

 xt=r*cos(theta)

 yt=r*sin(theta)

 plot(xt,yt,type="l")

Ex 2.1.(d)

iff<-function(a)

{

if (a>0)

{show(floor(((a-floor(a))*100)%%10))}

else if(a<0)

{show(floor(((ceiling(a)-a)*100)%%10))}

else {show(0)}

}

a<--2.358

iff(a)

5

Ex 2.2

(a)   a<-c(1:8,7:1)

(b)  rep(1:5,1:5)

(c)   matrix(rep(1,9),ncol=3)-diag(c(1,1,1))

(d)  matrix(c(0,2,3,0,5,0,7,0,0),ncol=3,byrow=T)

Ex.2.3

polar<-function(a)

{b<-rep(0,2)

r=sqrt(a[1]^2+a[2]^2)

if (a[1]==0||a[2]==0)

{theta=0}

else if (x>0||x==0)

{theta=asin(a[2]/r)}

else if (x<0)

{theta=-asin(a[2]/r)+pi}

b[1]=r;b[2]=theta;show(b)

}

a<-c(1,2)

polar(a)

2.236068 1.107149

Ex 2.4

a<-1:100

a[a%%2!=0]

 [1]  1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

[26] 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99

a<-1:100

 a[a%%2!=0|a%%3!=0]

 [1]   1   2   3   4   5   7   8   9  10  11  13  14  15  16  17  19  20  21  22

[20]  23  25  26  27  28  29  31  32  33  34  35  37  38  39  40  41  43  44  45

[39]  46  47  49  50  51  52  53  55  56  57  58  59  61  62  63  64  65  67  68

[58]  69  70  71  73  74  75  76  77  79  80  81  82  83  85  86  87  88  89  91

[77]  92  93  94  95  97  98  99 100

 a<-1:100

 a[a%%2!=0|a%%3!=0|a%%7!=0]

 [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19

[20]  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38

[39]  39  40  41  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58

[58]  59  60  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77

[77]  78  79  80  81  82  83  85  86  87  88  89  90  91  92  93  94  95  96  97

[96]  98  99 100

Ex 2.5

 queue <- c("Steve", "Russell", "Alison", "Liam")

 queue[queue!="Alison"]

[1] "Steve"   "Russell" "Liam"

 Generate a upper triangle matrix:

R codes:

GenerateW <- function(Dimension)

{

  result <- matrix(0, nrow=Dimension, ncol=Dimension)

  for( i in 1:Dimension)

  {

    result[i,i:Dimension] <- 1

  }

  return(result)

}

GenerateW(3)

     [,1] [,2] [,3]

[1,]    1    1    1

[2,]    0    1    1

[3,]    0    0    1

> 

How to use sink function in R:

Example:

sink("R001.txt")

x <- c(2,-6,-4,8,5,4,1,3,4,-9,0,1)

A <- matrix(x, ncol=3)

A

A.prima <- t(A)

A.prima

dim(A)

dim(A.prima)

sink()

Metropolis sampler:

1. Start with some initial value

2. Given this initial value, draw a candidate value  from some proposal distribution.

3. Compute the ratio  of the density at the candidate and initial points, .

4. With probability min(, accept the candidate point, else retain the current point.

5. Return to step 2.

Using R to do a  Metropolis sampler:

R codes:

sig1=0.5

sig2=0.1

sig3=10

n=500

x1=rep(0,n)

x2=x1

x3=x1

x1[1]=-10

x2[1]=0

x3[1]=0

for(i in 2:n)

{y=rnorm(1,0,1)*sig1+x1[i-1]

u=runif(1,0,1)

alpha=dnorm(y,0,1)/dnorm(x1[i-1],0,1)

if (u<=alpha)

{x1[i]=y}

else

x1[i]=x1[i-1]

}

plot(x1,type="l")

Plot:

Metropolis-Hastings sampler:

1. Start with some initial value

2. Given this initial value, draw a candidate value  from some proposal distribution.

3. Let  denote the probability of  given the current value is . Compute the ratio  of the density at the candidate and initial points, .

4. With probability min(, accept the candidate point, else retain the current point.

5. Return to step 2.

R codes which generated a Metropolis-Hastings sampler for a scaled inverse Chi- square distribution:

The MCMC codes deal with generating a Metropolis-Hastings sampler for a scaled inverse Chi-square distribution.

Recall that we can write the probability distribution for a scaled inverse Chi-square as:

Where n is the degrees of freedom, a is the scaling parameter. Here we will assume n=5,a=4. First, we need to define a function f(x,n,a)  in R:

f<-function(x,n,a) x^(-n/2)*exp(-a/(2*x))

First we need to specify a vector for the chain. Let’s call it met and give it length 5000.

met<-numeric(5000)

Lets now start the chain, say at value 1. We will call the current value in the chain last,

last<-1

We use a normal Chi-square distribution with 5 degrees of freedom. Then the codes will be:

for( i in 1:5000)

{ cand<-rchisq(1,5)

  alpha<-(f(cand,5,4)/f(last,5,4))*(dchisq(cand,5)/dchisq(last,5))

if  (runif(1)<min(alpha,1))

last<-cand

met[i]<-last}

plot(met,type=”l”)

Plot a histogram, and control the width of the columns using “lwd”.

a<-c(0,0,1,0,0,1.5,0,0,0,1.3,0,0,2.1)

plot(a,lwd=2,type="h")

a<-c(0,0,1,0,0,1.5,0,0,0,1.3,0,0,2.1)

plot(a,lwd=5,type="h")

Using R codes to write the logistic simulation:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\WinBugs\\bino\\gen.csv",sep=",",header=FALSE))

beta<-read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\WinBugs\\bino\\map1.csv",sep=",",header=TRUE)

bet<-beta$effect

n<-numeric(nrow(x))

r<-numeric(nrow(x))

p<-numeric(nrow(x))

for (i in 1:nrow(x))

{n[i]=rpois(1,32)

p[i]=x[i,]%*%bet

p[i]=exp(p[i])/(1+exp(p[i]))

r[i]=rbinom(1,n[i],p[i])

}

out<-matrix(c(n,r),ncol=2)

out

write(t(out),file="\\Users\\Soemone\\Desktop\\Guelph\\WinBugs\\bino\\simulationphe.csv",sep=",",ncolumns=2)

“write” can also output character strings.

Paste Greek

a<-c(1,2,3,2,3,4,5)

plot(a,main=expression(sigma==2))

or

a<-c(1,2,3,2,3,4,5)

plot(a,main=~sigma~"=2")

Result:

An example:

z1<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\example data sets\\highauto.csv",sep=",",header=FALSE))

z2<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\example data sets\\lowauto.csv",sep=",",header=FALSE))

par(mfrow=c(2,1))

plot(z2,type="h",lwd=3.7,xlab="Lag",ylab="Autocorrelation",ylim=c(0,1))

abline(h=0)

text(25,0.9,labels="(a)")

plot(z1,type="h",lwd=3.7,xlab="Lag",ylab="Autocorrelation",ylim=c(0,1))

text(25,0.9,labels="(b)")

abline(h=0)

Layout

nf<-layout(matrix(c(1,1,2,3), 2,2, byrow = TRUE))

layout.show(nf)

results:

Curve and histogram

x <- rnorm(1000);

hist(x, freq = FALSE, col = "white");

curve(dnorm, col = 2, add = TRUE)

Density estimate:

plot(density(x.norm),main="Density estimate of data")

Layout:

layout(matrix(c(1,1,2,3), 2,2, byrow = TRUE))

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\BayesianReview\\damage\\post-sample.csv",sep=",",header=TRUE))

y=x[,2]

par(omi=c(0.02,0.02,0.02,0.02))

plot(y,type="l", main="Diagnostics for beta",ylab="beta",xlab="iteration")

par(omi=c(0.02,0.02,2,2))

acf(y,ylab="Autocorrelation",main="",lag.max=50)

par(omi=c(0.1,1.3,1.3,0.1))

plot(density(y),type="l",ylab="Posterior Density",xlab="beta",main="")

layout(matrix(c(1,1,2,3), 2,2, byrow = TRUE))

x<as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\gamma.csv",sep=",",header=TRUE))

print(x)

y=x[,6]

z=x[,1]

plot(z,y,type="l", main="Diagnostics for beta",ylab="beta",xlab="iteration",axes = FALSE)

axis(1,at=c(10000,20000,30000,40000,50000),lab=c("10000","20000","30000","40000","50000"))

axis(2,at=c(1.0,1.5,2.0),lab=c("1.0","1.5","2.0"))

box()

# par(omi=c(0.02,0.02,2,2))

acf(y,ylab="Autocorrelation",main="",lag.max=50,ylim=c(-1,1))

# par(omi=c(0.1,1.3,1.3,0.1))

plot(density(y),type="l",ylab="Posterior Density",xlab="beta",main="")

hist(y,freq=FALSE,main="")

lines(density(y))

Histogram and density curve:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\BayesianReview\\wheat\\gamma22.csv",sep=",",header=TRUE))

y=x[,2]

z=x[,1]

hist(y,prob=TRUE,xlab="QTL effect",ylab="prob",main="",breaks=14)

lines(density(y,adjust=2),lty=1)

box()

text(-0.4,0.63,labels="(b)")

Column:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\interval.csv",sep=",",header=TRUE))

y=x[,6]

z=x[,1]

plot(z,y,type="h",lwd=3.7,xlab="",ylab="QTL effect",ylim=c(-1.0,1.5))

text(0,1.5,labels="(a)")

abline(h=0)

abline(v=c(73.51310938,109.4142975,174.3300143,192.46167),lty=3)

require(stats)  # normally loaded

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure3-2.csv",sep=",",header=FALSE))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome  Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(b)")

abline(h=0)

lines(x6,x1,type = "l",lty=1)

lines(x6,x2,type="l",lty=1)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=2)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(50,4.66,labels="100%")

text(145,3.36,labels="94%")

text(205,-2.44,labels="73%")

text(290,-1.78,labels="57%")

text(250,2.44,labels="100%")

text(360,3.36,labels="100%")

text(610,1.3,labels="12%")

text(630,-1.3,labels="15%")

text(780,0.97,labels="12%")

text(820,1.93,labels="99%")

text(905,4.01,labels="100%")

text(990,2.44,labels="83%")

text(1100,-1.5,labels="40%")

text(1220,-1.2,labels="34%")

text(1305,-2.44,labels="55%")

text(1335,1.78,labels="40%")

text(1425,1.2,labels="36%")

text(1405,-1.93,labels="35%")

text(1800,0.91,labels="12%")

text(2300,1.09,labels="26%")

abline(h=0)

text(0,6.7,labels="(b)")

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure7.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE")

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure10-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-5,5),xlim=c(0,350),axes = FALSE)

axis(1,at=c(0,50,100,150,200,250,300,350),lab=c("0","50","100","150","200","250","300","350"))

axis(2,at=c(-5,-4,-3,-2,-1,0,1,2,3,4,5),lab=c("-5","-4","-3","-2","-1","0","1","2","3","4","5"))

box()

text(0,5,labels="(a)")

abline(h=0)

abline(v=c(73.6,109.5,174.4,192.5),lty=2)

lines(x6,x2,type = "l",lty=2)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=1)

lines(x6,x5,type="l",lty=1)

Figure 1

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(a)")

abline(h=0)

Figure 2

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure2.csv",sep=",",header=T))

hist(x,prob=TRUE,xlab="QTL effect",ylab="Density",main="",breaks=37)

lines(density(x),lty=1)

box()

Figure 3

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure3-2.csv",sep=",",header=F))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome  Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(b)")

abline(h=0)

lines(x6,x1,type = "l",lty=1)

lines(x6,x2,type="l",lty=1)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=2)

Figure 4

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure4-3.csv",sep=",",header=F))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome  Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(c)")

abline(h=0)

lines(x6,x1,type = "l",lty=1)

lines(x6,x2,type="l",lty=1)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=2)

Figure 5 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(50,4.66,labels="100%")

text(145,3.36,labels="84%")

text(205,-2.44,labels="48%")

text(290,-1.78,labels="28%")

text(250,2.44,labels="100%")

text(360,3.36,labels="99%")

text(610,1.3,labels="1%")

text(630,-1.3,labels="5%")

text(780,0.97,labels="3%")

text(820,1.93,labels="98%")

text(905,4.01,labels="99%")

text(990,2.44,labels="63%")

text(1100,-1.5,labels="17%")

text(1220,-1.2,labels="12%")

text(1305,-2.44,labels="28%")

text(1335,1.78,labels="18%")

text(1425,1.2,labels="15%")

text(1405,-1.93,labels="17%")

text(1800,0.91,labels="8%")

text(2300,1.09,labels="12%")

abline(h=0)

text(0,6.7,labels="(a)")

Figure 5 (b)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(50,4.66,labels="100%")

text(145,3.36,labels="94%")

text(205,-2.44,labels="73%")

text(290,-1.78,labels="57%")

text(250,2.44,labels="100%")

text(360,3.36,labels="100%")

text(610,1.3,labels="12%")

text(630,-1.3,labels="15%")

text(780,0.97,labels="12%")

text(820,1.93,labels="99%")

text(905,4.01,labels="100%")

text(990,2.44,labels="83%")

text(1100,-1.5,labels="40%")

text(1220,-1.2,labels="34%")

text(1305,-2.44,labels="55%")

text(1335,1.78,labels="40%")

text(1425,1.2,labels="36%")

text(1405,-1.93,labels="35%")

text(1800,0.91,labels="12%")

text(2300,1.09,labels="26%")

abline(h=0)

text(0,6.7,labels="(b)")

Figure 6

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure6-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="False Positive Rate",ylim=c(0,1),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(0,0.2,0.4,0.6,0.8,1.0),lab=c("0.0","0.2","0.4","0.6","0.8","1.0"))

box()

abline(h=0)

abline(h=0.1,lty=2)

text(0,1,labels="(b)")

Figure 7

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure7.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1))

Figure 8 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure8-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-6,4),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-6,-4,-2,0,2,4),lab=c("-6","-4","-2","0","2","4"))

box()

text(0,4,labels="(a)")

abline(h=0)

abline(v=c(85.6,146.6,212.3,282.2),lty=2)

lines(x6,x1,type = "l",lty=2)

lines(x6,x2,type="l",lty=2)

lines(x6,x3,type="l",lty=1)

lines(x6,x4,type="l",lty=1)

Figure 8 (b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure8-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1),ylim=c(30,55))

text(0,55,labels="(b)")

Figure 9 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure9-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-1,2),xlim=c(0,1100),axes = FALSE)

axis(1,at=c(0,100,200,300,400,500,600,700,800,900,1000,1100),lab=c("0","100","200","300","400","500","600","700","800","900","1000","1100"))

axis(2,at=c(-1,0,1,2),lab=c("-1","0","1","2"))

box()

text(0,2,labels="(a)")

abline(h=0)

abline(v=c(140,320.8,486.1,632.4,836.2,946.4),lty=2)

lines(x6,x2,type = "l",lty=2)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=1)

lines(x6,x5,type="l",lty=1)

Figure 9 (b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure9-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1),ylim=c(18,21),axes = FALSE)

axis(1,at=c(0,0.2,0.4,0.6,0.8,1),lab=c("0.0","0.2","0.4","0.6","0.8","1.0"))

axis(2,at=c(18,19,20,21),lab=c("18","19","20","21"))

box()

text(0,21,labels="(b)")

Figure 10 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure10-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-5,5),xlim=c(0,350),axes = FALSE)

axis(1,at=c(0,50,100,150,200,250,300,350),lab=c("0","50","100","150","200","250","300","350"))

axis(2,at=c(-5,-4,-3,-2,-1,0,1,2,3,4,5),lab=c("-5","-4","-3","-2","-1","0","1","2","3","4","5"))

box()

text(0,5,labels="(a)")

abline(h=0)

abline(v=c(73.6,109.5,174.4,192.5),lty=2)

lines(x6,x2,type = "l",lty=2)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=1)

lines(x6,x5,type="l",lty=1)

Figure 10 (b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure10-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1),ylim=c(415,455),axes = FALSE)

axis(1,at=c(0,0.2,0.4,0.6,0.8,1),lab=c("0.0","0.2","0.4","0.6","0.8","1.0"))

axis(2,at=c(415,425,435,445,455),lab=c("415","425","435","445","455"))

box()

text(0,455,labels="(b)")

library(scatterplot3d)

xx<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\Compare QTL Methods\\analysis-inter\\sim-poi-blup-03.csv",sep=",",header=T))

x<-xx[,1]

y<-xx[,2]

z<-xx[,3]

scatterplot3d(y,x,z,color="blue",type="h",lwd=3,xlim=c(1,161),ylim=c(1,161))

library(lattice)

xx<-read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\Compare QTL Methods\\analysis-inter\\sim-poi-blup-03.csv",sep=",",header=T)

wireframe(effect~marker1*marker2,data=xx, color="blue",type="h",lwd=3)

Compare two methods:

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\Compare QTL Methods\\analysis-nointer\\simul2\\pseudo\\em-jeffreys-blup-poi-01.csv",sep=",",header=F))

x1=x[,7]

x2=x[,8]

plot(x2,x1,type="h",lwd=3,col=1,xlab="Genome Position(centiMorgan)",ylab="QTL effect",xlim=c(0,2400),ylim=c(-4,6),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-4,-2,0,2,4,6),lab=c("-4","-2","0","2","4","6"))

box()                                                                       

abline(h=0)

text(0,6,labels="(d)")

Another graph codes:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\PRESSBINO.csv",sep=",",header=TRUE))

x1=x[,1]

x2=x[,2]

x3=x[,3]

plot(x1,x2,type="o",pch=1,col=1,xlab="T Values",ylab="PRESS",ylim=c(14000,16500),xlim=c(-2,7),axes = FALSE)

axis(1,at=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7),lab=c("-2","-1.5","-1","-0.5","0","0.5","1","1.5","2","2.5","3","3.5","4","4.5","5","5.5","6","6.5","7"))

axis(2,at=c(14000,14500,15000,15500,16000,16500),lab=c("14000","14500","15000","15500","16000","16500"))

lines(x1,x3,type = "o",pch=2,col=2)

box()

legend(5,16250, c("Gneneral","Diag"),

       pch =c(1,2), col =1:2)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\sim-map-1.csv",sep=",",header=T))

x1=as.matrix(x[,3])

x2=as.matrix(x[,5])

plot(x1,x2,type="h",lwd=3.7,xlab="",ylab="QTL effect",ylim=c(-2,2),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2),lab=c("-2","-1.5","-1","-0.5","0","0.5","1.0","1.5","2.0"))

box()

text(0,2,labels="(b)")

abline(h=0)

Time and Date:

x <- as.Date("31DEC2000", format = "%d%b%Y")

y <- as.Date("31DEC1999", format = "%d%b%Y")

diff <- x - y

new.dt <- y + as.difftime(diff, unit="days")

Clean R space

rm(list = ls(all = TRUE))

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\bina.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

plot(x1,x2,type="h",lwd=8.0,lty=1,col="grey",xlab="Genome  Position(centiMorgan)",ylab="QTL effect",ylim=c(-1.5,3),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-1.5,0,1.5,3),lab=c("-1.5","0","1.5","3.0"))

box()

lines(x1,x3,type = "h",lwd=3.7,lty=1,col="black")

Figure 1

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

plot(x4,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-2.0,2.0),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400),lab=c("0","200","400","600","800","1000","1200","1400","1600","1800","2000","2200","2400"))

axis(2,at=c(-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0),lab=c("-2.0","-1.5","-1.0","-0.5","0.0","0.5","1.0","1.5","2.0"))

box()

abline(h=0)

abline(h=-2.0)

x10<-(x5!=0)

x11<-rep(-1.8,length(x4))

x12<-x10*x11

x13<-rep(-2.0,length(x4))

x14<-x10*x13

segments(x4,x14,x4,x12)

legend(1000,2,c("True"),lty=1,col=1,lwd=3.7)

Figure 2

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

plot(x4,x8,type="h",lwd=3.7,col=6,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-2.0,2.0),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400),lab=c("0","200","400","600","800","1000","1200","1400","1600","1800","2000","2200","2400"))

axis(2,at=c(-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0),lab=c("-2.0","-1.5","-1.0","-0.5","0.0","0.5","1.0","1.5","2.0"))

box()

abline(h=0)

abline(h=-2.0)

x11<-(x5!=0)

x12<-rep(-1.8,length(x4))

x13<-x11*x12

x14<-rep(-2.0,length(x4))

x15<-x11*x14

segments(x4,x15,x4,x13)

legend(1000,2,c("Overdispersion"),lty=1,col=6,lwd=3.7)

Figure 3

a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure3.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

plot(x5,x6,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-2.0,2.0),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0),lab=c("-2.0","-1.5","-1.0","-0.5","0.0","0.5","1.0","1.5","2.0"))

box()

abline(h=-2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x11<-(x3!=0)

x12<-rep(-1.8,length(x3))

x13<-x11*x12

x14<-rep(-2.0,length(x3))

x15<-x11*x14

lines(x5,x8,type="h",lwd=3.7,col=2)

segments(x5,x15,x5,x13)

legend(270,2,c("Expectation","Overdispersion"),lty=1,col=c(4,2),lwd=3.7)

b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure3.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

plot(x5,x7,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="LOD score",ylim=c(-20,600),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(0,100,200,300,400,500,600),lab=c("0","100","200","300","400","500","600"))

box()

abline(h=0)

abline(h=-20)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x11<-(x3!=0)

x12<-rep(-20,length(x3))

x13<-x11*x12

x14<-rep(0,length(x3))

x15<-x11*x14

lines(x5,x9,type="h",lwd=3.7,col=2)

segments(x5,x15,x5,x13)

legend(270,600,c("Expectation","Overdispersion"),lty=1,col=c(4,2),lwd=3.7)

Figure 4

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure4.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

bins=seq(-1.5,31.5,by=3)

hist(x3,breaks=bins,col="gray",main="",xlab="Number of seeded spikelets",axes = FALSE)

axis(1,at=seq(0,30,by=3),lab=c("0","3","6","9","12","15","18","21","24","27","30"))

axis(2,at=c(0,10,20,30,40,50,60,70,80,90),lab=c("0","10","20","30","40","50","60","70","80","90"))

box()

Figure 5

a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure5.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1

b=x5-2.9

plot(b,x6,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-0.2,1.4),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-0.2,0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4),lab=c("-0.2","0.0","0.2","0.4","0.6","0.8","1.0","1.2","1.4"))

box()

abline(h=-0.2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(-0.1,length(x3))

x14<-x12*x13

x15<-rep(-0.2,length(x3))

x16<-x12*x15

 lines(a,x8,type="h",lwd=3.7,col=2)

 lines(x5,x10,type="h",lwd=3.7,col=1)

segments(x5,x16,x5,x14)

legend(270,1.4,c("Expectation","Overdispersion","Mixture"),lty=1,col=c(4,2,1),lwd=3.7)

b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure5.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1

b=x5-2.9

plot(b,x7,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="LOD score",ylim=c(-2,14),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(0,2,4,6,8,10,12,14),lab=c("0","2","4","6","8","10","12","14"))

box()

abline(h=-2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(-1,length(x3))

x14<-x12*x13

x15<-rep(-2,length(x3))

x16<-x12*x15

 lines(a,x9,type="h",lwd=3.7,col=2)

 lines(x5,x11,type="h",lwd=3.7,col=1)

segments(x5,x16,x5,x14)

legend(270,14,c("Expectation","Overdispersion","Mixture"),lty=1,col=c(4,2,1),lwd=3.7)

Figure 6

a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure6.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1.2

plot(a,x6,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-0.5,0.8),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-0.4,-0.2,0.0,0.2,0.4,0.6,0.8),lab=c("-0.4","-0.2","0.0","0.2","0.4","0.6","0.8"))

box()

abline(h=-0.5)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(-0.45,length(x3))

x14<-x12*x13

x15<-rep(-0.5,length(x3))

x16<-x12*x15

 lines(x5,x8,type="h",lwd=3.7,col=2)

segments(x5,x16,x5,x14)

legend(270,0.8,c("Expectation","Overdispersion"),lty=1,col=c(4,2,1),lwd=3.7)

 

b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure6.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1.3

plot(a,x7,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-2,70),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(0,10,20,30,40,50,60,70),lab=c("0","10","20","30","40","50","60","70"))

box()

abline(h=-2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(0,length(x3))

x14<-x12*x13

x15<-rep(-2,length(x3))

x16<-x12*x15

 lines(x5,x9,type="h",lwd=3.7,col=2)

segments(x5,x16,x5,x14)

legend(270,70,c("Expectation","Overdispersion"),lty=1,col=c(4,2),lwd=3.7)