############################################################################################################
# Loading the first data set for testing.
############################################################################################################
myd<-read.csv("Simulated Choice data revised simplified.dat")
# Wileyto analysis to derive k for second stage
myd$ip1<-1-myd$LLmagn/myd$SSmagn
myrWileyto<-by(myd, myd$Subject, function(x) coef(glm(Waited~ip1+LLdelay-1, data=x, family=binomial)))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(myrWileyto)
## Length Class Mode
## 1 2 -none- numeric
## 2 2 -none- numeric
## 3 2 -none- numeric
## 4 2 -none- numeric
## 5 2 -none- numeric
## 6 2 -none- numeric
## 7 2 -none- numeric
## 8 2 -none- numeric
## 9 2 -none- numeric
## 10 2 -none- numeric
## 11 2 -none- numeric
## 12 2 -none- numeric
## 13 2 -none- numeric
## 14 2 -none- numeric
## 15 2 -none- numeric
## 16 2 -none- numeric
## 17 2 -none- numeric
## 18 2 -none- numeric
## 19 2 -none- numeric
## 20 2 -none- numeric
## 21 2 -none- numeric
## 22 2 -none- numeric
## 23 2 -none- numeric
## 24 2 -none- numeric
## 25 2 -none- numeric
## 26 2 -none- numeric
## 27 2 -none- numeric
## 28 2 -none- numeric
## 29 2 -none- numeric
## 30 2 -none- numeric
## 31 2 -none- numeric
## 32 2 -none- numeric
## 33 2 -none- numeric
## 34 2 -none- numeric
## 35 2 -none- numeric
## 36 2 -none- numeric
## 37 2 -none- numeric
## 38 2 -none- numeric
## 39 2 -none- numeric
## 40 2 -none- numeric
## 41 2 -none- numeric
## 42 2 -none- numeric
## 43 2 -none- numeric
## 44 2 -none- numeric
## 45 2 -none- numeric
## 46 2 -none- numeric
## 47 2 -none- numeric
## 48 2 -none- numeric
## 49 2 -none- numeric
## 50 2 -none- numeric
# The following commands perform a logistic regression on each subject and extracts the parameters.
# This looks complicated because I'm extracting information for a bunch of individual regressions.
# The multilevel approach is much more straightforward.
library(plyr)
fitList <- dlply(myd, .(Subject), glm, formula=Waited~ip1+LLdelay-1, family=binomial)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
jack<-as.data.frame(t(sapply(fitList,coef)))
jack$k<-jack$LLdelay/jack$ip1
jack$Subject<-as.numeric(row.names(jack))
jack$logk<-log(jack$k+1) # problem with some negative k values; this captures all but one of them
## Warning in log(jack$k + 1): NaNs produced
# Setting it up and testing for a group comparison
jack$Group<-ifelse(jack$Subject < 26, 1, 0)
jackr<-lm(logk~Group, data=jack)
summary(jackr)
##
## Call:
## lm(formula = logk ~ Group, data = jack)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6740 -0.3327 -0.1004 0.3904 2.2022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.3345 0.1568 2.133 0.038384 *
## Group 0.8611 0.2291 3.758 0.000491 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7839 on 45 degrees of freedom
## (3 observations deleted due to missingness)
## Multiple R-squared: 0.2389, Adjusted R-squared: 0.2219
## F-statistic: 14.12 on 1 and 45 DF, p-value: 0.0004909
library(lattice)
densityplot(~logk|Group, data=jack)
wilcox.test(k~Group, data=jack) # nonparametric test of full set of 50 k values
##
## Wilcoxon rank sum exact test
##
## data: k by Group
## W = 220, p-value = 0.07408
## alternative hypothesis: true location shift is not equal to 0
###### Multilevel regression of the same data - the recommended analysis
library(lme4)
## Loading required package: Matrix
myr<-glmer(Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1+(log(LLmagn/SSmagn)+log(LLdelay+1)-1|Subject), data=myd, family=binomial)
summary(myr)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## Waited ~ log(LLmagn/SSmagn) + log(LLdelay + 1) - 1 + (log(LLmagn/SSmagn) +
## log(LLdelay + 1) - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 3371.8 3404.4 -1680.9 3361.8 4995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.3671 -0.3322 -0.0404 0.2722 8.3049
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject log(LLmagn/SSmagn) 1.819 1.349
## log(LLdelay + 1) 1.672 1.293 0.17
## Number of obs: 5000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## log(LLmagn/SSmagn) 2.3247 0.2103 11.057 <2e-16 ***
## log(LLdelay + 1) -1.9261 0.1981 -9.723 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## l(LL/S
## lg(LLdly+1) 0.022
# Generating plots like those that appear in the paper.
plot(betaM~I(-betaD), data=myd, xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]), xlim=c(-4, .4), ylim=c(0,5))
plot(coef(myr)$Subject[[1]]~coef(myr)$Subject[[2]], cex=jack$logk+.5, xlab=expression(beta["Delay"]), ylab=expression(beta["Magnitude"]), xlim=c(-4, .4), ylim=c(0,5))
# four points aren't being plotted due to NaN or strong negative k
coefSub <- coef(myr)$Subject[c(1,4,9),]
points(coefSub[[1]]~coefSub[[2]], pch = 3)
jill<- data.frame(aggregate(myd$betaM, by=list(myd$Subject), FUN=mean),aggregate(myd$betaD, by=list(myd$Subject), FUN=mean))
colnames(jill)<-c("Subject", "betaM", "S", "betaD")
jill<-jill[c(1,2,4)]
jill$fit.betaM <- coef(myr)$Subject[[1]]
jill$fit.betaD <- coef(myr)$Subject[[2]]
# Plot estimated by actual regression weights
par(mfrow=c(1,2))
plot(jill$fit.betaM~jill$betaM, xlim=c(0,5), ylim=c(0,5), xlab=expression(beta["Magnitude"]), ylab=expression(paste("Fitted ", beta["Magnitude"])))
plot(function(x) x, 0,5, add=T, lty=2)
plot(jill$fit.betaD~ I(-jill$betaD), xlim=c(-4, .25), ylim=c(-4,.25), xlab=expression(beta["Delay"]), ylab=expression(paste("Fitted ", beta["Delay"])))
plot(function(x) x, -4,.25, add=T, lty=2)
# Same plots but for fits of each individual rather than a multilevel model.
myrindivid<-by(myd, myd$Subject, function(x) coef(glm(Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1, data=x, family=binomial)))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
summary(myrindivid)
## Length Class Mode
## 1 2 -none- numeric
## 2 2 -none- numeric
## 3 2 -none- numeric
## 4 2 -none- numeric
## 5 2 -none- numeric
## 6 2 -none- numeric
## 7 2 -none- numeric
## 8 2 -none- numeric
## 9 2 -none- numeric
## 10 2 -none- numeric
## 11 2 -none- numeric
## 12 2 -none- numeric
## 13 2 -none- numeric
## 14 2 -none- numeric
## 15 2 -none- numeric
## 16 2 -none- numeric
## 17 2 -none- numeric
## 18 2 -none- numeric
## 19 2 -none- numeric
## 20 2 -none- numeric
## 21 2 -none- numeric
## 22 2 -none- numeric
## 23 2 -none- numeric
## 24 2 -none- numeric
## 25 2 -none- numeric
## 26 2 -none- numeric
## 27 2 -none- numeric
## 28 2 -none- numeric
## 29 2 -none- numeric
## 30 2 -none- numeric
## 31 2 -none- numeric
## 32 2 -none- numeric
## 33 2 -none- numeric
## 34 2 -none- numeric
## 35 2 -none- numeric
## 36 2 -none- numeric
## 37 2 -none- numeric
## 38 2 -none- numeric
## 39 2 -none- numeric
## 40 2 -none- numeric
## 41 2 -none- numeric
## 42 2 -none- numeric
## 43 2 -none- numeric
## 44 2 -none- numeric
## 45 2 -none- numeric
## 46 2 -none- numeric
## 47 2 -none- numeric
## 48 2 -none- numeric
## 49 2 -none- numeric
## 50 2 -none- numeric
fitList <- dlply(myd, .(Subject), glm, formula=Waited~log(LLmagn/SSmagn)+log(LLdelay+1)-1, family=binomial)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
julie<-as.data.frame(t(sapply(fitList,coef)))
julie$Subject<-as.numeric(row.names(julie))
summary(julie)
## log(LLmagn/SSmagn) log(LLdelay + 1) Subject
## Min. :-0.4077 Min. :-4.7275 Min. : 1.00
## 1st Qu.: 1.2147 1st Qu.:-3.1650 1st Qu.:13.25
## Median : 2.5886 Median :-2.1300 Median :25.50
## Mean : 2.9527 Mean :-2.0443 Mean :25.50
## 3rd Qu.: 3.5339 3rd Qu.:-0.9842 3rd Qu.:37.75
## Max. :28.9181 Max. : 0.2210 Max. :50.00
julie$betaM<-jill$betaM
julie$betaD<-jill$betaD
plot(julie[[1]]~jill$betaM, xlim=c(0,5), ylim=c(0,6), xlab=expression(beta["Magnitude"]), ylab=expression(paste("Fitted ", beta["Magnitude"])))
plot(function(x) x, 0,6, add=T, lty=2)
plot(julie[[2]]~ I(-jill$betaD), xlim=c(-4, .25), xlab=expression(beta["Delay"]), ylab=expression(paste("Fitted ", beta["Delay"])))
plot(function(x) x, -4,.25, add=T, lty=2)
plot(coef(myr)$Subject[[1]]~row.names(coef(myr)$Subject))
plot(coef(myr)$Subject[[2]]~as.numeric(row.names(coef(myr)$Subject))%%5)
# Testing ability to identify strength of an effect that was held constant
myr1.1<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 1), family=binomial)
myr1.2<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 2), family=binomial)
myr1.3<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 4), family=binomial)
myr1.4<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 8), family=binomial)
myr1.5<-glmer(Waited~log(LLmagn/SSmagn)+(log(LLmagn/SSmagn)|Subject), data=subset(myd, LLdelay == 16), family=binomial)
summary(myr1.1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + (log(LLmagn/SSmagn) | Subject)
## Data: subset(myd, LLdelay == 1)
##
## AIC BIC logLik deviance df.resid
## 920.1 944.7 -455.1 910.1 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4471 -0.5228 0.1026 0.4239 2.7325
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 0.2285 0.4781
## log(LLmagn/SSmagn) 2.2801 1.5100 0.28
## Number of obs: 1000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2185 0.1577 -7.725 1.12e-14 ***
## log(LLmagn/SSmagn) 2.3901 0.2920 8.185 2.71e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## lg(LLmg/SS) -0.324
summary(myr1.2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + (log(LLmagn/SSmagn) | Subject)
## Data: subset(myd, LLdelay == 2)
##
## AIC BIC logLik deviance df.resid
## 908.7 933.3 -449.4 898.7 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3333 -0.4405 -0.0588 0.3785 4.1378
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 2.019 1.421
## log(LLmagn/SSmagn) 1.313 1.146 0.33
## Number of obs: 1000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8600 0.2744 -6.778 1.21e-11 ***
## log(LLmagn/SSmagn) 2.0308 0.2383 8.522 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## lg(LLmg/SS) -0.206
summary(myr1.3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + (log(LLmagn/SSmagn) | Subject)
## Data: subset(myd, LLdelay == 4)
##
## AIC BIC logLik deviance df.resid
## 761.9 786.4 -375.9 751.9 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9583 -0.3444 -0.1047 0.2529 5.0112
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 4.358 2.088
## log(LLmagn/SSmagn) 3.127 1.768 0.10
## Number of obs: 1000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.3780 0.4462 -7.571 3.70e-14 ***
## log(LLmagn/SSmagn) 2.5795 0.3577 7.212 5.51e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## lg(LLmg/SS) -0.396
summary(myr1.4)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + (log(LLmagn/SSmagn) | Subject)
## Data: subset(myd, LLdelay == 8)
##
## AIC BIC logLik deviance df.resid
## 681.8 706.3 -335.9 671.8 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.4753 -0.2062 -0.0728 0.1714 4.8021
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 11.322 3.365
## log(LLmagn/SSmagn) 1.125 1.061 -0.08
## Number of obs: 1000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.5568 0.7402 -6.156 7.47e-10 ***
## log(LLmagn/SSmagn) 2.5384 0.3452 7.353 1.94e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## lg(LLmg/SS) -0.649
summary(myr1.5)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + (log(LLmagn/SSmagn) | Subject)
## Data: subset(myd, LLdelay == 16)
##
## AIC BIC logLik deviance df.resid
## 595.5 620.1 -292.8 585.5 995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1902 -0.1895 -0.0880 0.0935 5.4720
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 8.610 2.934
## log(LLmagn/SSmagn) 2.628 1.621 0.27
## Number of obs: 1000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.5283 0.7644 -5.924 3.14e-09 ***
## log(LLmagn/SSmagn) 1.5517 0.4670 3.323 0.000891 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## lg(LLmg/SS) -0.462
############################################################################################################
# Testing extension of above to make a group comparison (Subjects 1-25 vs. 26-50)
############################################################################################################
myd$Group<-ifelse(myd$Subject < 26, 1, 0)
# Multilevel Fit
myr2<-glmer(Waited~Group*log(LLmagn/SSmagn)+Group*log(LLdelay+1)-1-Group+(log(LLmagn/SSmagn)+log(LLdelay+1)-1|Subject), data=myd, family=binomial)
summary(myr2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ Group * log(LLmagn/SSmagn) + Group * log(LLdelay + 1) -
## 1 - Group + (log(LLmagn/SSmagn) + log(LLdelay + 1) - 1 | Subject)
## Data: myd
##
## AIC BIC logLik deviance df.resid
## 3328.5 3374.2 -1657.3 3314.5 4993
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.8625 -0.3317 -0.0390 0.2741 10.4762
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject log(LLmagn/SSmagn) 0.5186 0.7202
## log(LLdelay + 1) 1.6481 1.2838 0.18
## Number of obs: 5000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## log(LLmagn/SSmagn) 3.4547 0.2060 16.772 < 2e-16 ***
## log(LLdelay + 1) -1.8204 0.2715 -6.704 2.03e-11 ***
## Group:log(LLmagn/SSmagn) -2.2363 0.2661 -8.403 < 2e-16 ***
## Group:log(LLdelay + 1) -0.2194 0.3884 -0.565 0.572
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## l(LL/S l(LL+1 G:(LL/
## lg(LLdly+1) -0.068
## Grp:(LL/SS) -0.767 0.044
## Grp:l(LL+1) 0.041 -0.688 -0.042
xyplot(fitted(myr2)~LLdelay|Group, data=myd, type="a")
# Group comparison for individually derived weights generated in earlier analysis.
julie$Group<-ifelse(julie$Subject < 26, 1, 0)
colnames(julie)[1]<-"fit.betaM"
colnames(julie)[2]<-"fit.betaD"
julierM<-lm(fit.betaM~Group, data=julie)
julierD<-lm(fit.betaD~Group, data=julie)
summary(julierD)
##
## Call:
## lm(formula = fit.betaD ~ Group, data = julie)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.54914 -1.12641 -0.07801 1.07657 2.39938
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.9102 0.2866 -6.665 2.4e-08 ***
## Group -0.2682 0.4053 -0.662 0.511
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.433 on 48 degrees of freedom
## Multiple R-squared: 0.009037, Adjusted R-squared: -0.01161
## F-statistic: 0.4377 on 1 and 48 DF, p-value: 0.5114
summary(julierM)
##
## Call:
## lm(formula = fit.betaM ~ Group, data = julie)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.540 -1.332 -0.606 0.538 24.322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5962 0.7434 6.183 1.32e-07 ***
## Group -3.2870 1.0513 -3.127 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.717 on 48 degrees of freedom
## Multiple R-squared: 0.1692, Adjusted R-squared: 0.1519
## F-statistic: 9.776 on 1 and 48 DF, p-value: 0.003001
densityplot(~fit.betaM|Group, data=julie)
# Try again by omitting subject 50
julierM<-lm(fit.betaM~Group, data=subset(julie, Subject!=50))
summary(julierM)
##
## Call:
## lm(formula = fit.betaM ~ Group, data = subset(julie, Subject !=
## 50))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.7168 -0.7229 -0.1857 0.5805 2.1868
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.5828 0.2040 17.56 < 2e-16 ***
## Group -2.2736 0.2856 -7.96 2.93e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9995 on 47 degrees of freedom
## Multiple R-squared: 0.5741, Adjusted R-squared: 0.565
## F-statistic: 63.36 on 1 and 47 DF, p-value: 2.925e-10
densityplot(~fit.betaM|Group, data=julie)
############################################################################################################
# Trying a new data set with nonzero SS delays.
############################################################################################################
myd2<-read.csv("Simulated Choice data revised simplified nonzero.dat")
myr2<-glmer(Waited~log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1+(log(LLmagn/SSmagn)+log((LLdelay+1)/(SSdelay+1))-1|Subject), data=myd2, family=binomial)
summary(myr2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Waited ~ log(LLmagn/SSmagn) + log((LLdelay + 1)/(SSdelay + 1)) -
## 1 + (log(LLmagn/SSmagn) + log((LLdelay + 1)/(SSdelay + 1)) -
## 1 | Subject)
## Data: myd2
##
## AIC BIC logLik deviance df.resid
## 3550.3 3582.8 -1770.1 3540.3 4995
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -15.9597 -0.3347 0.0119 0.3329 11.1544
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject log(LLmagn/SSmagn) 1.766 1.329
## log((LLdelay + 1)/(SSdelay + 1)) 2.115 1.454 0.17
## Number of obs: 5000, groups: Subject, 50
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## log(LLmagn/SSmagn) 2.4943 0.2093 11.915 <2e-16 ***
## log((LLdelay + 1)/(SSdelay + 1)) -1.9384 0.2193 -8.839 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## l(LL/S
## l((LL+1)/+1 0.038
xyplot(fitted(myr2)~SSdelay, group=Subject, data=myd2, type="a", ylab="P(LL)")
