coef(lme.strat)
fixef(lme.strat)
fake
print(lme.no.strat)
print(lme.strat)
fake.rs <- melt(fake, id=c("person1",  "treatment.ind", "time"), measure="y")
cast(fake.rs,  trreatment.ind ~ ., function(x) c(M=mean(x), SD=sd(x) , N=length(x))
str(fake.rs)
ranef(lme.strat)
var(lme.strat[[1]][1])
lme.strat[[1]][1]
lme.strat[[1]]
lme.strat[1]
ranef(lme.strat$placebo.ind)
xx <- ranef(lme.strat)
str(xx)
xx@.Data
xx@.Data[[1]][1:5]
xx@.Data$placebo.ind
xx@.Data$[[1]]
xx@.Data[[1]]
xxx<-xx@.Data[[1]]
str(xxx)
xxx[1, 1:5]
xxx[1:5,1]
xx <- ranef(lme.strat)@.Data[[1]][1:5,1]
xx
var(ranef(lme.strat)@.Data[[1]][1:5,1])
var(ranef(lme.strat)@.Data[[2]][1:5,1])   # variance of treatment.ind
var(ranef(lme.strat)@.Data[[2]][6:10, 1])   # variance of treatment.ind
cast(fake.rs,  treatment.ind ~ ., function(x) c(M=mean(x), SD=sd(x) , N=length(x)))
fake.rs <- melt(fake, id=c("person1",  "treatment.ind", "time"), measure="y")#
cast(fake.rs,  treatment.ind ~ ., function(x) c(M=mean(x), SD=sd(x) , N=length(x)), margins=TRUE)
J = 10; K = 100#
set.seed(500)#
fake = Strat.mean.var.simple (J,K)#
#
fake.rs <- melt(fake, id=c("person1",  "treatment.ind", "time"), measure="y")#
cast(fake.rs,  treatment.ind ~ ., function(x) c(M=mean(x), SD=sd(x) , N=length(x)), margins=TRUE)
cast(fake.rs,  treatment.ind ~ ., function(x) c(M=mean(x), VAR=var(x) , N=length(x)))
J = 10; K = 100
J = 10; K = 5#
set.seed(500)#
#
    time <- rep(seq(0,1, ,length=K), J) # K measurements#
    person <- rep(1:(J/2), each = K)#
    treatment <- rep(0:1, each = J/2)#
    treatment.ind <- rep(0:1, each = (J/2)*K)#
    person1 <- rep(1:J, ,each = K)#
    placebo.ind.1 <- treatment.ind < 1#
    placebo.ind = ifelse( placebo.ind.1, 1, 0)
   data.frame (y, time, person1, treatment.ind, placebo.ind)
   data.frame (time, person1, treatment.ind, placebo.ind)
    sigma.y.true = 1.2
    a.true.P = rnorm (J/2, mu.a.true.P, sigma.a.true.P)#
    a.true.T = rnorm (J/2, mu.a.true.T, sigma.a.true.T)
    mu.a.true.P = 4.8#
    sigma.a.true.P = 2.2#
#
# ... M and SD of patients in treatment group#
    mu.a.true.T = 8.8#
    sigma.a.true.T = 4.2#
 #
 # ...  sd of residual#
    sigma.y.true = 1.2
   a.true.P = rnorm (J/2, mu.a.true.P, sigma.a.true.P)#
    a.true.T = rnorm (J/2, mu.a.true.T, sigma.a.true.T)
    y.P <- rnorm( (J/2)*K,  a.true.P[person],  sigma.y.true)#
    y.T <- rnorm( (J/2)*K,  a.true.T[person],  sigma.y.true)
y.P
y.T
a.true.P[person]
person
a.true.P
    y.T <- rnorm( (J/2)*K,  a.true.T[person],  sigma.y.true)
atrue.T
a.true.T
a.true.T[person]
?glmer
str(cbpp)
?nlmer
(m1 <- lmer(cbind(incidence, size - incidence) ~ period + (1 | herd),#
            family = binomial, data = cbpp))
RSiteSearch("glmer lmer")
x=data.frame(resp=rbinom(100,1,0.8),  treat=sample(letters[7:8],100,TRUE),  pat=as.factor(sample(letters[1:4],100,TRUE))  )  
x=data.frame(resp=rbinom(100,1,0.8),  treat=sample(letters[7:8],100,TRUE),  pat=as.factor(sample(letters[1:4],100,TRUE))
x=data.frame(resp=rbinom(100,1,0.8),  treat=sample(letters[7:8],100,TRUE),  pat=as.factor(sample(letters[1:4],100,TRUE))))
x=data.frame(resp=rbinom(100,1,0.8), treat=sample(letters[7:8],100,TRUE),pat=as.factor(sample(letters[1:4],100,TRUE))))
x=data.frame( resp=rbinom(100,1,0.8), treat=sample(letters[7:8],100,TRUE), pat=as.factor(sample(letters[1:4],100,TRUE))))
x
x=data.frame( resp=rbinom(100,1,0.8), treat=sample(letters[7:8],100,TRUE), pat=factor(sample(letters[1:4],100,TRUE))))
x=data.frame( resp=rbinom(100,1,0.8), treat=sample(letters[7:8],100,TRUE), pat=as.factor(sample(letters[1:4],100))))
x=data.frame( resp=rbinom(100,1,0.8), treat=sample(letters[7:8],100,TRUE))
x$pat <- sample(letters[1:4],100,TRUE)
str(x)
lmer(resp~treat+(1|pat),family=binomial,data=x)
lmer(resp~treat+(1|pat),data=x)
hist(x$resp)
lmer(resp~treat+(1|pat),family=poisson,data=x)
glmer(resp~treat+(1|pat),family=binomial,data=x)
f <- c(rep(1:10,2))#
gain <- rnorm(20, mean = 30 )#
phase <- rnorm(20, mean = 10)#
coh <- runif(20)#
sys <- c(rep("system1",10), rep("system2",10))#
sys_df <- data.frame( Frequency = f, Gain = gain, Phase = phase,#
Coherence = coh, System = sys)
library(reshape)#
library(ggplot2)#
melted_sys_df <- melt(sys_df, id.var = c("Frequency", "System") )#
ggtheme(theme_bw)
ggplot(melted_sys_df, aes(x = Frequency, y = value, colour = System)) +#
        geom_line() + facet_grid(variable ~ .) + scale_x_log10()
bode <-  ggplot(sys_df, aes(x = Frequency, colour = System)) +#
                  geom_line() + scale_x_log10()#
bode_gain <- bode + aes(y = Gain)#
bode_phase <- bode + aes(y = Phase)
bode_phase
bode_gain
n = data.frame(id = c(1,2,3,4), parent = c(0,1,2,2), value =#
c(5,5.5,7,3), date = c(1,2,3,3.5))#
#
plot(n$date, n$value)#
#
do.call(#
  segments,#
  with(#
    merge(n,n,by.x="parent", by.y="id"),#
    data.frame(x0=date.x, y0=value.x, x1=date.y, y1=value.y)#
  )#
)
df <- data.frame(id = c(1,2,3,4), parent = c(0,1,2,2), value =#
c(5,5.5,7,3), date = c(1,2,3,3.5))#
#
selfmerge <- merge(df, df, by.x="parent", by.y="id")#
#
qplot(date, value, data=df) +#
geom_segment(data=selfmerge, aes(x=date.x, xend=date.y, y=value.x,#
yend=value.y))
selfmerge
plot(n$date, n$value)#
#
do.call(#
  segments,#
  with(#
    merge(n,n,by.x="parent", by.y="id"),#
    data.frame(x0=date.x, y0=value.x, x1=date.y, y1=value.y)#
  )#
)
library(sem)
citations(sem)
citation(sem)
?sem
rm(list=ls())#
#
# #
library(lme4)#
set.seed(1234321)#
#
load("taf.X.rda")
dir()
library(TeachingDemos)#
op <- par(mfrow = c(3,3), ## split region#
          oma = c(5,0,4,0) + 0.1, ## create outer margin#
          mar = c(5,4,2,2) + 0.1) ## shrink some margins #
plot(1:10, main = "a", pch = 1:2, col= 1:2) #
plot(1:10, main = "b", pch = 1:2, col= 1:2) #
tmp1 <- cnvrt.coords( 0.5, 0, input='plt' )$tdev # save location for#
mtext#
plot(1:10, main = "c", pch = 1:2, col= 1:2) #
plot(1:10, main = "d", pch = 1:2, col= 1:2) #
plot(1:10, main = "e", pch = 1:2, col= 1:2) #
plot(1:10, main = "f", pch = 1:2, col= 1:2) #
plot(1:10, main = "g", pch = 1:2, col= 1:2) #
plot(1:10, main = "h", pch = 1:2, col= 1:2) #
plot(1:10, main = "i", pch = 1:2, col= 1:2) #
## title #
mtext("My Plots", side = 3, outer = TRUE, font = 2, line = 1, cex = 1.2,#
#
	at=tmp1$x) #
## draw legend #
par(xpd=NA)#
tmp2 <- cnvrt.coords( tmp1$x, 0.05, input='tdev' )$usr # get location#
for legend#
legend(tmp2$x, tmp2$y, legend = c("Type 1", "Type 2"), #
 	pch = 1:2, col = 1:2, ncol = 2, xjust=0.5, yjust=0.5)#
par(op)
for legend {legend(tmp2$x, tmp2$y, legend = c("Type 1", "Type 2"), #
 	pch = 1:2, col = 1:2, ncol = 2, xjust=0.5, yjust=0.5)}#
par(op)
for legend {legend(tmp2$x, tmp2$y, legend = c("Type 1", "Type 2"),  pch = 1:2, col = 1:2, ncol = 2, xjust=0.5, yjust=0.5)}
install.packages("Matrix", repos = "http://r-forge.r-project.org", type="source")
install.packages("Matrix")
install.packages("lme4", repos = "http://r-forge.r-project.org", type="source")
library(lme4, lib.loc=.libUsr)
library(gridBase)
?grid
med = median(sleepstudy$Reaction)#
sleepstudy$bin = (sleepstudy$Reaction > med) + 0#
mod = lmer(bin ~ Days + (1|Subject) + (0+Days|Subject),#
            data = sleepstudy, family = binomial)#
#
ilog = function(x) { 1/(1 + exp(-x)) }#
boxplot(ilog(fitted(mod)) ~ bin, data = sleepstudy)
fitted(mod)
str(sleepstudy)
(mod <- lme(measure ~ nitro*year,random= ~ 1|bloc/nitro))#
anova(mod)#
summary(aov(measure ~ nitro*year + Error(bloc/nitro)))
(mod <- lmer(measure ~ nitro*year,random= ~ 1|bloc/nitro))#
anova(mod)#
summary(aov(measure ~ nitro*year + Error(bloc/nitro)))
RSiteSearch("npk")
aov1 = aov(yield ~  N+P+K + Error(block), npk)
summary(aov1)
(mod <- lmer(measure ~ nitro*year,random= ~ 1|bloc/nitro), npk)#
anova(mod)#
summary(aov(measure ~ nitro*year + Error(bloc/nitro), npk))
(mod <- lmer(measure ~ nitro*year,random= ~ 1|bloc/nitro), npk)#
anova(mod)#
summary(aov(yield ~ nitro*year + Error(bloc/nitro), npk))
str(npk)
(mod <- lmer(measure ~ nitro*year,random= ~ 1|bloc/nitro), npk)#
anova(mod)#
summary(aov(yield ~ N*P*K + Error(1/nitro), npk))
example(aov)
as.matrix(summary(npk.aov)[[1]])exampl
str(npk.aov)
as.matrix(summary(npk.aov)[[2]])
summary(npk.aov)
summary(npk.aov)[[1]]
as.matrix(summary(npk.aov)[[1]])
effects(npk.aov[[3]])
effects(npk.aov[[2]])
effects(fits[[2]])
effects(fitted(npk.aov)[[2]])
effects(fits(npk.aov)[[2]])
example(aov) #
as.data.frame.summary.aovlist <- function(x) { #
   if(length(x) == 1) { #
     as.data.frame(x[[1]]) #
   } else { #
     lapply(unlist(x, FALSE), as.data.frame) #
   } #
} #
x1 <- summary(npk.aov) #
x2 <- summary(npk.aovE) #
as.data.frame(x1) #
as.data.frame(x2)
library("glmmADMB", lib.loc=.libUsr)
?glmm.admb
data(epil2)#
  glmm.admb(y~Base*trt+Age+Visit,random=~Visit,group="subject",data=epil2,family="nbinom")
data(epil2)#
  glmm.admb(y~Base*trt+Age+Visit,random=~1,group="subject",data=epil2,family="nbinom")
install.packages("ggplot2")
point
library(ggplot2)
?geom_smooth
?scale_fill_identity
?plotmatrix
plotmatrix(mtcars[, 1:3])
plotmatrix(mtcars[, 1:3]) + geom_smooth(method="lm")
plotmatrix(mtcars[, 1:3]) + geom_smooth(method="smooth")
plotmatrix(mtcars[, 1:3]) + geom_smooth()
p <- ggplot(diamonds, aes(x=carat, y=..density..)) + geom_histogram(binwidth=0.2)
p <- ggplot(diamonds, aes(x=carat, y=..mean..))
p <- ggplot(diamonds, aes(x=carat, y=..mean..)) + geom_line()
p
str(p)
?geom_line
fph <- 0.4 #
Sigh <- sqrt(0.0002) #
Sigi <- sqrt(0.04) #
#
reH <- rnorm(90, fph, Sigh)  ## hospid effects
reH <- rnorm(90, fph, Sigh)  ## hospid effects
set.seed(7658943)
fph <- 0.4 #
Sigh <- sqrt(0.0002) #
Sigi <- sqrt(0.04)
ls
ls()
?rnorm
reh <- rnorm(90, fph, Sigh)  ## hospid effects
rnorm(90, fph, Sigh)  ## hospid effects
rnorm(90, fph, Sigh)
reh <- rnorm(90, fph, Sigh)
dta <- within(expand.grid(hospid = 1:90, empid = 1:80), #
            fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi))
reH <- rnorm(90, fph, Sigh)
dta <- within(expand.grid(hospid = 1:90, empid = 1:80), fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi))
dta
(fm1 <- lmer(fpi1 ~ (1|hospid), dta))
?within
library(lme4)#
set.seed(7658943) #
#
fph <- 0.4 #
Sigh <- sqrt(0.0002) # Sigma between hospitals#
Sigi  <- sqrt(0.04)     # Sigma between individuals (residual)#
#
reH <- rnorm(90, fph, Sigh)  # random effect for Hospital#
data <- within(expand.grid(hospid = 1:90, empid = 1:80), #
                     fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi)) #
#
(fm1 <- lmer(fpi1 ~ (1|hospid), data))
library(lme4,lib.loc=.libUsr)
str(Dyestuff)
VarCorr(fm1)
(fm1 <- lmer(Reaction ~ Days + Subject + (Days|Subject), sleepstudy))
(fm2 <- lmer(Reaction ~ Days + Subject + (Days|Subject), sleepstudy))
print(fm2 <- lmer(Reaction ~ Days + Subject + (Days|Subject), sleepstudy),cor=FALSE)
ranef(fm1)
ranef(fm2)
lmer(score ~ Machine + (1|Worker/Machine), Machine)
library(MEMSS)
lmer(score ~ Machine + (1|Worker/Machine), Machines)
str(Machines)
Machines
269*64*8*3
1242/269
library(lme4)#
set.seed(7658943)#
reH <- rnorm(90, sd = sqrt(0.0002))#
dta <- within(data.frame(hospid = gl(90,80)),#
              fpi1 <- reH[hospid] + rnorm(length(hospid), 0.4, sqrt(0.04)))#
dotplot(reorder(hospid, fpi1) ~ fpi1, dta)#
(fm1 <- lmer(fpi1 ~ 1|hospid, dta, verb = TRUE))
install.packages("Matrix",repos="http://r-forge.r-project.org")
install.packages("Matrix",repos="http://r-forge.r-project.org", type="source")
library(Matrix)
install.packages("lme4",repos="http://r-forge.r-project.org", type="source")
sessionInfo()
install.packages("Matrix",repos="http://r-forge.r-project.org", type="source", lib.loc=.libUsr)
?install.packages
install.packages("Matrix",repos="http://r-forge.r-project.org", type="source", lib=.libUsr)
library(NRAIA)
?NRAIA
library(NRAIA, helo)
library(MASS)
round(cor(Animals),2)
round(cor(log(Animals)),2)
plot(Animals)
plot(log(Animals))
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))#
(fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy))
vc <- VarCorr(fm1)
vc
varcomps<-c(unlist( lapply(vc, diag) ), # random intercept variances#
              attr(vc,"sc")^2)            # residual variance
str(varcomps)
fm1
varcomps
(m1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),#
            family = binomial, data = cbpp))
library(lattice)
xyplot(val ~ loc | mouse, data = df,#
       groups=valtype,#
       type=c("p","l"),#
       distribute.type = TRUE,#
#
       col=c("black", "blue"))
df
24921/88579
22372/80631
22192/80138
88579-24921
(88579-24921)/12
library(nlme)
?predict.lme
fm3 <- lme(distance ~ age*Sex, data = Orthodont, random = ~ 1)#
df3.1 <- with(Orthodont, data.frame(age=seq(5, 20, 5),#
                    Subject=rep(Subject[1], 4),#
                    Sex=rep(Sex[1], 4)))
summary(fm3)
library(gmodels)
ci(fm3)
detach(package:nlme)
str(Orthodont)
xx <- data(nlme::Orthodont)
xx <- data(nlme:Orthodont)
xx <- data(Orthodont)
ci(fm1)
install.packages("multcomp")
library(multcomp)
library(lme4)
?lmer
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
glht(fm1)
?glht
vignette("generalsiminf", package = "multcomp"
)
library(lattice)#
library(lme4)#
#
# ######### Chunk 1 #######################
# Generate a set of data for a 2 x 5 design, only subjects as random factor#
set.seed(1000)#
rFreq <- rnorm(2, mean=5, sd=0.5)#
rCond <- rep(rFreq / 1.1, 2)#
nbase <- rnorm(240, 10, 0.1)#
data <- within(expand.grid(freq=LETTERS[1:2], cond=letters[1:4],#
                          id=factor(1:30)), {#
                             n <- rFreq[as.numeric(freq)] + rCond[as.numeric(cond)] + nbase#
                          })#
data <- data[order(data$freq, data$cond), ]#
## Simulate an interaction#
data$n[data$freq == "A"] <- rev(data$n[data$freq == "A"])#
## Have a look#
xyplot(n ~ cond | freq, data=data, groups=id,#
       type="b", pch=19, cex=0.3)#
#
tapply(data$n, list(data$cond, data$freq), mean)
print(n.lmer2 <- lmer(n ~ freq*cond + (1 | id), data), cor=FALSE)
library("lme4");#
#
#
bdf_lmer1<-lmer(bdi~bdi.pre+time+treatment+drug+length+(1|subject),method="ML",na.action=na.omit,data=bdf_long)#
#
#
#
summary(bdf_lmer1)#
#
#
#
bdf_lmer2<-lmer(bdi~bdi.pre+time+treatment+drug+length+(time|subject),method="ML",na.action=na.omit,data=bdf_long)#
#
#
#
summary(bdf_lmer2)#
#
#
#
anova(bdf_lmer1, bdf_lmer2)
?bdf_long
library(reshape)#
library(ggplot2)#
#
Year <- 2007:1996
all <-      c(329, 421, 393, 372, 276, 254, 262, 219, 227, 204, 208, 165)#
D.abs <-    c(25, 30, 25, 22, 25, 13, 16, 10, 7, 5, 5, 11)#
UK.abs <-   c(39, 53, 48, 48, 41, 27, 21, 26, 22, 16, 16, 12)#
NL.abs <-   c(16, 20, 17, 25, 14, 11, 7, 4, 2, 6, 3, 3)#
CA.abs <-   c(31, 42, 46, 40, 27, 27, 33, 31, 23, 19, 18, 18)#
ANZ.abs <-  c(21, 18, 26, 19, 13, 15, 16, 14, 9, 10, 6, 4)  # Australia, New Zealand#
F.abs <-    c(8, 14, 14, 14, 7, 6, 11, 3, 5, 7, 4, 2) #
AS.abs <-   c(5, 5, 2, 4, 1, 0, 1, 0, 0, 1, 0, 0)           # Austria, Switzerland#
SC.abs <-   c(6, 4, 6, 3, 7, 2, 6, 4, 1, 3, 3, 1)           # SE, FI, NO, DK, IL#
SE.abs <-   c(12, 20, 11, 5, 10, 9, 4, 6, 4, 4, 0, 4)       # Italy, Spain, P, Greece#
ASIA.abs <- c(7, 15, 7, 4, 5, 5, 6, 3, 1, 4, 2, 0)          # China, Taiwan, Japan, India#
IS.abs <-   c(5, 13, 20, 9, 14, 11, 4, 5, 5, 4, 7, 2)       # Israel#
#
d <- data.frame(Year, all, D.abs, UK.abs, NL.abs, CA.abs, ANZ.abs, F.abs, AS.abs,#
	    SC.abs, SE.abs, ASIA.abs,IS.abs,#
	    D=D.abs/all, UK=UK.abs/all, NL=NL.abs/all, CA=CA.abs/all, ANZ=ANZ.abs/all,#
	    F=F.abs/all, AS=AS.abs/all, SC=SC.abs/all, SE=SE.abs/all, ASIA=ASIA.abs/all, IS=IS.abs/all)#
#
d$US <- 1-rowSums(d[,14:24])	#
ix <- which(d$Year >= 1997 & d$Year <=2006)  # 10-year slot#
d <- d[ix, ]#
d.rs <- melt(d, id="Year", measure=c("D", "UK", "CA", "ANZ", "IS"), var="Country")#
qplot(Year, value, data=d.rs, colour=Country, geom=c("point", "smooth"), method=lm, formula=y~poly(x,2),#
       ylab="Proportion of Publications")
d.rs <- melt(d, id="Year", measure=c("D.abs", "UK.abs",  "CA.abs", "ANZ.abs", "F.abs", "NL.abs"), var="Country")#
qplot(Year, value, data=d.rs, colour=Country, geom=c("point", "smooth"), method=lm, formula=y~poly(x,2))
setwd('/Users/kliegl/documents/manuscripts/APS')
d.rs <- melt(d, id="Year", measure=c("D", "UK", "CA", "ANZ", "US", "IS"), var="Country")#
qplot(Year, value, data=d.rs, colour=Country, geom=c("point", "smooth"), method=lm, formula=y~poly(x,2),#
       ylab="Proportion of Publications")
?read
?read.table
a <- read.delim2("PsyCountryTime.txt")
nrow(a)
str(a)
hist(a$Citable.Documents)
a <- read.delim2("PsyCountryTime.txt", dec=",")
a <- read.delim1("PsyCountryTime.txt", dec=",")
setwd('/Users/kliegl/documents/manuscripts/APS/SCIImago')
a <- read.delim("PsyCountryTime.txt", dec=",")
a
