Gregor Gorjanc (gg)

<<>>=
1 + 1
@
\begin{Schunk}

\begin{Sinput}
> 1 + 1
\end{Sinput}

\begin{Soutput}
[1] 2
\end{Soutput}

\end{Schunk}
someFunc <- function(x, ...) ... lm(y ~ someFunc(x)) 
## Add my repository
tmp <- c("http://gregor.gorjanc.googlepages.com",           
         contrib.url(repos=getOption("repos"),
         type=getOption("pkgType")))

## Install the package
install.packages("animSci", contriburl=tmp, dep=TRUE)

## Load the package
library(package="animSci")

## Run the lactation curve examples
example(topic="lactCurves")
## Needed packages and functions
library(package="e1071")
library(package="doBy")
library(package="R2WinBUGS")

rvalue <- function(n, p, value, mean=0, sdE=0)
{
  if(!is.matrix(value)) value <- as.matrix(value)
 
  ## Frequencies
  k <- length(p)
  q <- 1 - p
  P <- p * p
  Q <- q * q
  H <- 2 * p * q

  ## Allocate matrix of arbitrary values by loci
  x <- matrix(data=1, nrow=n, ncol=k)
  y <- matrix(data=1, nrow=n, ncol=k)

  ## Simulate arbitrary values
  for(i in 1:k) {
    x[, i] <- rdiscrete(n=n, probs=c(P[i], H[i], Q[i]))
    y[, i] <- value[x[, i], i]
  }
  x <- x - 1 ## convert 1:3 to 0:2, i.e., number of A1 alleles

  ## Sum all values
  y <- rowSums(y)

  ## Add mean and normal deviate (error)
  y <- y + rnorm(n=n, mean=mean, sd=sdE)

  return(as.data.frame(cbind(y, x)))
}

prepairData <- function(x, dominance=FALSE)
{
  ret <- list()
  ret$n <- nrow(x)
  ret$y <- x$y
  if(!dominance) {
    ret$x <- x$V2 - 1 ## scale to the mean
  } else {
    ret$x <- x$V2 + 1 ## BUGS does not allow zero indexing!
  }
  ret
}

### SIMULATE DATA
###-----------------------------------------------------------------------

## --- Purely additive ---

set.seed(seed=19791123)

## Model parameters
n <- 1000
mu <- 100
sigma <- 1
value <- c(-1, 0, 1)
p <- 0.5

## Simulate
podatkiA <- rvalue(n=n, p=p, value=value, mean=mu, sdE=sigma)

par(mfrow=c(2, 1))
hist(podatkiA$y, xlab="Phenotype")

plot(podatkiA$y ~ jitter(podatkiA$V2), xlab="Genotype", ylab="Phenotype")

tmpA1 <- lm(podatkiA$y ~ podatkiA$V2)
abline(coef(tmpA1), lwd="2", col="blue")
tmpA2 <- summaryBy(y ~ V2, data=podatkiA)
points(tmpA2$y.mean ~ tmpA2$V2, pch=19, col="red")

## --- Additive + Dominance ---

set.seed(seed=19791123)

## Model parameters
n <- 1000
mu <- 100
sigma <- 1
value <- c(-1, 1, -1)
p <- 0.5

## Simulate
podatkiD <- rvalue(n=n, p=p, value=value, mean=mu, sdE=sigma)

par(mfrow=c(2, 1))
hist(podatkiD$y, xlab="Phenotype")

plot(podatkiD$y ~ jitter(podatkiD$V2), xlab="Genotype", ylab="Phenotype")

tmpD1 <- lm(podatkiD$y ~ podatkiD$V2)
abline(coef(tmpD1), lwd="2", col="blue")
tmpD2 <- summaryBy(y ~ V2, data=podatkiD)
points(tmpD2$y.mean ~ tmpD2$V2, pch=19, col="red")

### ESTIMATE MODEL PARAMETERS AND DIC
###-----------------------------------------------------------------------

modelA <- function()
{
  ## --- Additive model (regression on gene content) ---
  ## Phenotypes
  for(i in 1:n) {
    y[i] ~ dnorm(mu[i], tau2)
    mu[i] <- a + b * x[i]
  }
  ## Location prior
  a ~ dnorm(0, 1.0E-6)  
  b ~ dnorm(0, 1.0E-6)
  ## Variance prior
  lsigma ~ dunif(-10, 10) ## sigma between 4.539993e-05 and 2.202647e+04
  tau2 <- pow(sigma, -2)
  sigma <- exp(lsigma)
}

modelD <- function()
{
  ## --- Additive + Dominance model (genotype as a factor) ---
  ## Phenotypes
  for(i in 1:n) {
    y[i] ~ dnorm(mu[i], tau2)
    mu[i] <- g[x[i]]
  }
  ## Location prior
  g[1] ~ dnorm(0, 1.0E-6)
  g[2] ~ dnorm(0, 1.0E-6)
  g[3] ~ dnorm(0, 1.0E-6)  
  ## Variance prior
  lsigma ~ dunif(-10, 10) ## sigma between 4.539993e-05 and 2.202647e+04
  tau2 <- pow(sigma, -2)
  sigma <- exp(lsigma)
}

tmpAA <- prepairData(x=podatkiA)
tmpAD <- prepairData(x=podatkiA, dominance=TRUE)
tmpDA <- prepairData(x=podatkiD)
tmpDD <- prepairData(x=podatkiD, dominance=TRUE)

initsA <- function() list(a=rnorm(n=1, mean=mu), b=rnorm(n=1), lsigma=0)
initsD <- function() list(g=rnorm(n=3, mean=mu),               lsigma=0)

(fitAA <- bugs(model=modelA, data=tmpAA, inits=initsA,
               n.chains=3, n.burnin=5000, n.iter=10000, n.thin=10,
               bugs.seed=19791123, parameters=c("a", "b", "sigma"),
               bugs.directory="C:/Programs/BUGS/WinBUGS14_orig"))
##            mean  sd   2.5%    25%  50%  75%  97.5% Rhat n.eff
## a          99.9 0.0   99.9   99.9  100  100  100.0    1  1500
## b           1.0 0.0    0.9    1.0    1    1    1.1    1  1500
## sigma       1.0 0.0    1.0    1.0    1    1    1.0    1   400
## deviance 2828.3 2.4 2826.0 2826.0 2828 2829 2834.0    1  1500
## pD = 2.9 and DIC = 2831.2

(fitAD <- bugs(model=modelD, data=tmpAD, inits=initsD,
               n.chains=3, n.burnin=5000, n.iter=10000, n.thin=10,
               bugs.seed=19791123, parameters=c("g", "sigma"),
               bugs.directory="C:/Programs/BUGS/WinBUGS14_orig"))
##            mean  sd   2.5%    25%    50%  75%  97.5% Rhat n.eff
## g[1]       98.9 0.1   98.8   98.9   98.9   99   99.1    1   660
## g[2]       99.9 0.0   99.8   99.9   99.9  100  100.0    1  1500
## g[3]      101.0 0.1  100.8  100.9  101.0  101  101.1    1  1500
## sigma       1.0 0.0    1.0    1.0    1.0    1    1.0    1  1300
## deviance 2829.2 2.9 2826.0 2827.0 2829.0 2831 2837.0    1  1500
## pD = 4.0 and DIC = 2833.2

## pD seem OK and DIC agrees with simulated data - it is worse
## for AD model for data following additive model

(fitDA <- bugs(model=modelA, data=tmpDA, inits=initsA,
               n.chains=3, n.burnin=5000, n.iter=10000, n.thin=10,
               bugs.seed=19791123, parameters=c("a", "b", "sigma"),
               bugs.directory="C:/Programs/BUGS/WinBUGS14_orig"))
##            mean  sd   2.5%    25%    50%    75%  97.5% Rhat n.eff
## a         100.0 0.0   99.9  100.0  100.0  100.0  100.1    1   570
## b           0.0 0.1   -0.1    0.0    0.0    0.1    0.2    1  1500
## sigma       1.4 0.0    1.3    1.4    1.4    1.4    1.5    1  1300
## deviance 3512.9 2.6 3510.0 3511.0 3512.0 3514.0 3520.0    1  1500
## pD = 3.0 and DIC = 3515.9

(fitDD <- bugs(model=modelD, data=tmpDD, inits=initsD,
               n.chains=3, n.burnin=5000, n.iter=10000, n.thin=10,
               bugs.seed=19791123, parameters=c("g", "sigma"),
               bugs.directory="C:/Programs/BUGS/WinBUGS14_orig"))
##            mean  sd   2.5%    25%    50%  75%  97.5% Rhat n.eff
## g[1]       98.9 0.1   98.8   98.9   99.0   99   99.1    1  1500
## g[2]      100.9 0.1  100.8  100.9  100.9  101  101.0    1  1500
## g[3]       99.0 0.1   98.8   98.9   99.0   99   99.1    1  1500
## sigma       1.0 0.0    1.0    1.0    1.0    1    1.0    1   530
## deviance 2829.2 2.9 2826.0 2827.0 2829.0 2831 2837.0    1  1500
## pD = 4.1 and DIC = 2833.3

## pD seems OK and DIC also agrees with the model --> dominance
## model is better for this dataset
library(package="MASS")
Sigma <- matrix(c(10, 10, 10, 20), nrow=2)
mu <- c(100, 100)
tmp <- mvrnorm(n=200, mu=mu, Sigma=Sigma, empirical=TRUE)
plot(tmp, xlab="X variable (unit)", ylab="Y variable (unit)") 
par(bty="l", pty="m", mar=c(3, 3, 1, 1), mgp=c(1.75, 0.75, 0))
plot(tmp, xlab="X variable (unit)", ylab="Y variable (unit)")
podatki <- read.table(file="corn.txt", na.strings=".", header=TRUE)
cols <- c(rgb(red=204, blue=0,   green=0,   max=255),
          rgb(red=0,   blue=153, green=0,   max=255),
          rgb(red=0,   blue=0,   green=204, max=255),
          rgb(red=204, blue=0,   green=153, max=255))
par(bty="l", pty="m", mar=c(5, 4, 1, 1))
matplot(x=podatki$YR, y=podatki[, c("IHP", "ILP", "IHO", "ILO")], type="l", lty=1,
        xlab="Year", ylab="Concentration (%)", col=cols[c(1, 1, 2, 2)], lwd=2)
legend("topleft", c("Protein", "Fat"), lty=1,
       lwd=2, col=cols[c(1, 2)], bty="n")

[[1]]
0.7071068

[[2]]
1.224745*x

[[3]]
-0.7905694 + 2.371708*x^2

[[4]]
-2.806243*x + 4.677072*x^3

[[5]]
0.7954951 - 7.954951*x^2 + 9.280777*x^4
x <- 1:100
leg4 <- as.matrix(as.data.frame(polynomial.values(polynomials=leg4coef,
                                                  x=scaleX(x, u=-1, v=1))))
colnames(leg4) <- c("leg0", "leg1", "leg2", "leg3", "leg4")
leg4 <- leg4[, 2:ncol(leg4)]
Legendre <- function(x, n, normalized=TRUE, intercept=FALSE, rescale=TRUE)
{
  ## Create a design matrix for Legendre polynomials
  ## x - numeric
  ## n - see orthopolynom
  ## normalized - logical, see orthopolynom
  ## intercept - logical, add intercept
  tmp <- legendre.polynomials(n=n, normalized=normalized)
  if(!intercept) tmp <- tmp[2:length(tmp)]
  polynomial.values(polynomials=tmp, x=x, matrix=TRUE)
}

polynomial.values <- function(polynomials, x, matrix=FALSE)
{
  ## Changed copy of polynomial.vales from orthopolynom in order
  ## to add matrix argument
  require(polynom)
  n <- length(polynomials)
  if(!matrix) {
    values <- vector(mode="list", length=n)
  } else {
    values <- matrix(ncol=n, nrow=length(x))
  }
  j <- 1
  while(j <= n) {
    if(!matrix) {
      values[[j]] <- predict(polynomials[[j]], x)
    } else {
      values[, j] <- predict(polynomials[[j]], x)
    }
    j <- j + 1
  }
  values
}

2009-12-19

2009-12-13

2009-11-20

2009-11-16

2009-11-15

2009-11-12

2009-10-13

2009-10-07

2009-10-02

2009-09-22

2009-09-21

2009-09-20

2009-09-10

2009-09-04

2009-09-03

2009-08-26

2009-08-08

2009-08-07

2009-07-30

2009-07-16

2009-07-14

2009-07-06

2009-07-03

2009-06-27

2009-06-19

2009-06-16

2009-06-12

2009-06-11

2009-06-09

2009-06-04

2009-05-27

2009-05-25

2009-05-02

2009-04-30

2009-04-23

2009-04-21

2009-04-19

2009-04-08

2009-04-03

2009-04-01

2009-03-31

2009-03-24

2009-03-21

2009-03-20

2009-03-17

2009-03-14

2009-03-13

2009-03-12

2009-03-11

2009-03-10

2009-03-06

2009-02-22

2009-02-20

2009-02-16

2009-02-15

2009-02-11

2009-02-10

2009-02-06