[R] mgcv gam bs=re questions
. For
example, consider s(SessIDX):factor(PID)1 0.008853821 above. Clearly
this is not what I want—it is not a measure of whether the five cases
differ significantly from each other in trend. I am guessing this random
effect concerns whether time points within PID1 differ significantly
from each other?
One might, of course, use gamm or gamm4 to address the random effects
question. I had hoped, however, to run all analyses within the same
program (gam) to avoid confounding differences in results with
differences in the estimation routines that different programs use.
Running it all in gam, however, encounters another problem, that it does
not appear that gam can use bs=”re” while fixing the trend to be linear.
For example:
M4 - gam(DVY ~ s(SessIDX, bs = re, k = 2, fx=TRUE) +
+factor(TX),
+data = SID5DVID1,
+family = quasipoisson(link=log), method=REML)
summary(M4,all.p=TRUE)
Error in b$smooth[[m]]$S[[1]] : subscript out of bounds
gam.vcomp(M4)
Error in if (x$smooth[[i]]$sp[j] 0) { : argument is of length zero
I have read extensively to try to understand all this (e.g., Wood 2013
“A simple test for random effects in regression models”) but have not
succeeded. Any help or advice would be appreciated very very much.
Will Shadish
--
William R. Shadish
Distinguished Professor
Founding Faculty
Mailing Address:
William R. Shadish
University of California
School of Social Sciences, Humanities and Arts
5200 North Lake Rd
Merced CA 95343
Physical/Delivery Address:
University of California Merced
ATTN: William Shadish
School of Social Sciences, Humanities and Arts
Facilities Services Building A
5200 North Lake Rd.
Merced, CA 95343
209-228-4372 voice
209-228-4007 fax (communal fax: be sure to include cover sheet)
wshad...@ucmerced.edu
http://faculty.ucmerced.edu/wshadish/index.htm
http://psychology.ucmerced.edu
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[R] metafor package
If your command looks something like this: forest(x,xlim=c(-6,7),cex=.8,order=prec, xlab=Effect Size, alim=c(-.5,4),mlab=RE Average for All Studies) In alim=c(-.5,4), increase the numbers. Make them higher than your highest confidence interval upper bound and lower than your lowest confidence interval lower bound. The bars may run into the text, but that may be fixed by stretching the output window in R. Will Shadish ** Date: Tue, 11 Feb 2014 22:11:20 + From: Nathan Pace n.l.p...@utah.edu To: r-help@R-project.org r-help@r-project.org Subject: [R] metafor package Message-ID: cf1fef68.2490a%n.l.p...@utah.edu Content-Type: text/plain Hi, I have a random effects meta analysis of a proportion (logit transformation) using rma.glmm. I have created a forest plot of the proportion (inverse logic transformation) using forest.rma. I have added the credibility interval. The forest plot is saved as a pdf. The dotted line and whiskers of the credibility interval are too faint. I need help on the argument(s) to widen the credibility interval dots and whiskers. I have looked at the forest.default function, but don?t see anything obvious to me. Nathan -- Nathan Pace, MD, MStat Department of Anesthesiology University of Utah 801.581.6393 n.l.p...@utah.edu -- William R. Shadish Distinguished Professor Founding Faculty Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4007 fax (communal fax: be sure to include cover sheet) wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] package mgcv - predict with bam: Error in X[ind, ] :, subscript out of bounds
Dear Simon, your note below says bs=re specifies a Gaussian random effect . I have been using bs = re for data modeled with Poisson and binomial distributions, or variants thereof (e.g., quasi-Poisson). Have I erred in assuming bs =re can be used to obtain random effects for such data? Will Shadish - Actually this is ok. mgcv exploits the duality between quadratically penalized smooths and Gaussian random effects to allow random effects to be specified this way. bs=re specifies a Gaussian random effect with corresponding model matrix given by model.matrix(~site-1). (More generally s(x,y,z,bs=re) specifies a gaussian random effect with model matrix given by model.matrix(~x:y:z-1), with obvious generalization to more or fewer variables). See mgcv help file ?random.effects for more. best, Simon -- William R. Shadish Distinguished Professor Founding Faculty Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4007 fax (communal fax: be sure to include cover sheet) wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] significance of random effect in mgcv gam
Thank you, Simon. On 12/4/2013 1:02 AM, Simon Wood wrote: Question. Am I correct that p = .126 above can be taken as the p-value associated with the random effect? - Yes. See http://biomet.oxfordjournals.org/content/100/4/1005.abstract for details of the approximate test used. On 03/12/13 20:09, William Shadish wrote: Dear R-helpers, I would like to test whether a random effect is significant when implemented with bs=re in mgcv gam. For example, if I run: M3b - gam(DVY ~ s(SessIDX, fTX, bs = re) + factor(TX), data = PCP, family = quasipoisson(link=log), method=REML) summary(M3b,all.p=TRUE) gam.vcomp(M3b) I obtain the the following output: summary(M3b,all.p=TRUE) Family: quasipoisson Link function: log Formula: DVY ~ s(SessIDX, fTX, bs = re) + factor(TX) Parametric coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 1.3282 0.2244 5.920 2.74e-07 *** factor(TX)1 -1.0546 0.7210 -1.463 0.15 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(SessIDX,fTX) 1.052 2 1.138 0.126 R-sq.(adj) = 0.388 Deviance explained = 39.5% REML score = 37.67 Scale est. = 1.4172n = 54 gam.vcomp(M3b) Standard deviations and 0.95 confidence intervals: std.dev lower upper s(SessIDX,fTX) 0.07842742 0.01095655 0.5613865 scale 1.19029872 0.97816911 1.4484316 Rank: 2/2 Question. Am I correct that p = .126 above can be taken as the p-value associated with the random effect? Thanks. Will Shadish -- William R. Shadish Distinguished Professor Founding Faculty Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4007 fax (communal fax: be sure to include cover sheet) wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] mgcv gam modeling trend variation over cases
Dear R-Helpers, I posted two days ago on testing significance of random effects in mgcv, but realize I did not make my overall purpose clear. I have a series of N short time series, where N might range from 3-10 and short means a median of 20 time points. The sample data below (PCP) has N = 4 cases with 9, 13, 16 and 16 observations over time respectively. The data set contains four variables: PID = case number, SessIDX = time (x axis), DVY = the outcome (y axis) and TX is treatment is applied = 1 or not = 0. My goal is to determine (a) is trend present in each case, (b) is it linear or nonlinear in each case, and (c) does trend vary significantly over cases (PID), the latter presumably to be measured with a random effect. I can do the first two (not shown here), but am not sure about the third. I am using generalized additive models, either mgcv gam or gamm4. For example, using mgcv gam my syntax is M1 - gam(DVY ~ s(SessIDX, bs = re) + factor(TX), data = PCP, family = quasipoisson(link=log), method=REML) summary(M1,all.p=TRUE) gam.vcomp(M1) Using gamm4, my syntax is PCP$fPID - factor(PCP$PID) M2 - gamm4(DVY ~ factor(TX) + s(SessIDX, by = factor(PID)), data = PCP, random =~ (1|fPID), family = poisson (link=log)) summary(M2$gam) summary(M2$mer) It is not clear to me whether either of these gives me what I want. In generalized linear mixed models, I am accustomed to the HLM approach (e.g., Raudenbush Bryk) where each case would have a trend coefficient, and the random effect would tell me if those four coefficients varied significantly. So that is what I am looking for, but adding the nonlinearity modeling of GAM. Is either of these formulations giving me what I want--a test of whether trend differs significantly over cases or not. Thanks for any help you can offer. I have worked very hard to solve this on my own, and just can't seem to do so. Will Shadish PCP PID SessIDX DVY TX 11 1 4 0 21 2 5 0 31 3 5 0 41 4 8 0 51 5 3 1 61 6 0 1 71 7 0 1 81 8 0 1 91 9 0 1 10 2 1 2 0 11 2 2 2 0 12 2 3 4 0 13 2 4 4 0 14 2 5 4 0 15 2 6 2 0 16 2 7 3 0 17 2 8 1 1 18 2 9 2 1 19 2 10 3 1 20 2 11 1 1 21 2 12 0 1 22 2 13 0 1 23 3 1 7 0 24 3 2 3 0 25 3 3 2 0 26 3 4 5 0 27 3 5 3 0 28 3 6 4 0 29 3 7 3 0 30 3 8 0 1 31 3 9 3 1 32 3 10 2 1 33 3 11 0 1 34 3 12 0 1 35 3 13 2 1 36 3 14 0 1 37 3 15 1 1 38 3 16 1 1 39 4 1 3 0 40 4 2 1 0 41 4 3 1 0 42 4 4 0 0 43 4 5 1 0 44 4 6 2 0 45 4 7 3 0 46 4 8 0 0 47 4 9 1 0 48 4 10 0 1 49 4 11 0 1 50 4 12 0 1 51 4 13 0 1 52 4 14 0 1 53 4 15 0 1 54 4 16 0 1 -- William R. Shadish Distinguished Professor Founding Faculty Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4007 fax (communal fax: be sure to include cover sheet) wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] significance of random effect in mgcv gam
Dear R-helpers, I would like to test whether a random effect is significant when implemented with bs=re in mgcv gam. For example, if I run: M3b - gam(DVY ~ s(SessIDX, fTX, bs = re) + factor(TX), data = PCP, family = quasipoisson(link=log), method=REML) summary(M3b,all.p=TRUE) gam.vcomp(M3b) I obtain the the following output: summary(M3b,all.p=TRUE) Family: quasipoisson Link function: log Formula: DVY ~ s(SessIDX, fTX, bs = re) + factor(TX) Parametric coefficients: Estimate Std. Error t value Pr(|t|) (Intercept) 1.3282 0.2244 5.920 2.74e-07 *** factor(TX)1 -1.0546 0.7210 -1.463 0.15 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Approximate significance of smooth terms: edf Ref.df F p-value s(SessIDX,fTX) 1.052 2 1.138 0.126 R-sq.(adj) = 0.388 Deviance explained = 39.5% REML score = 37.67 Scale est. = 1.4172n = 54 gam.vcomp(M3b) Standard deviations and 0.95 confidence intervals: std.dev lower upper s(SessIDX,fTX) 0.07842742 0.01095655 0.5613865 scale 1.19029872 0.97816911 1.4484316 Rank: 2/2 Question. Am I correct that p = .126 above can be taken as the p-value associated with the random effect? Thanks. Will Shadish -- William R. Shadish Distinguished Professor Founding Faculty Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4007 fax (communal fax: be sure to include cover sheet) wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
[R] two questions about xyplot
Dear R helpers,
I am generating three artificial short interrupted time series datasets
(single-case designs; call them Case 1, Case 2, Case 3) and then plotting
them in xyplot. I will put the entire code below so you can reproduce. I
have been unable to figure out how to do two things.
1. Each time series has 24 time points divided into four phases. Call them
phases A, B, C, D for convenience. I have used running() to compute the
means of the observations in each of these four parts; and I saved these as
objects called mn1 (for Case 1), mn2 (Case 2) and mn3 (Case 3). So mn1
contains the four means for A, B, C, D phases for Case 1, etc. I want to
insert these means into the xyplot in the appropriate place. For instance,
insert the first mean from mn1 into phase A of Case 1, the second mean into
phase B of Case 1, and so forth until insert the fourth mean from mn3 into
phase D of Case 3. Ideally, it would insert something like M = 49.02 or
Xbar = 49.02 into phase A for Case 1.
2. The xyplot code I use creates a line connecting the data points, and that
line is continuous over the entire graph. I would like to have the lines be
discontinuous between phases. Phase changes are indicated by panel.abline()
in the code, and occur at time points 6.5, 12.5, and 18.5. So, for example,
I would like a line connecting the datapoints from 1 to 6, then 7-12, then
13-18, then 19-24 (but not including 6-7, 12-13, and 18-19).
I appreciate any help you might be able to offer.
Will Shadish
Here is the code:
#
library(gtools)
library(lattice)
###g = .66
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
###change mean = to vary the effect size
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(1,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d1 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn1 - running(dvy, width=6, by=6)
#second dataset
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(2,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d2 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn2 - running(dvy, width=6, by=6)
#third dataset
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(3,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d3 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn3 - running(dvy, width=6, by=6)
#concatenate d1 d2 d3
d66 - rbind(d1, d2, d3)
d66df - as.data.frame(d66)
d66df$case - ordered(d66df$pid,
levels = c(1,2,3),
labels = c(Case 3, Case 2, Case 1))
p-xyplot(dvy ~ sessidx | case, data=d66df,
layout=c(1, 3), xlab= Sessions,
ylab = Number of Seconds,
type=l)
update(p, panel=function(...){
panel.xyplot(...)
panel.abline(v=6.5)
panel.abline(v=12.5)
panel.abline(v=18.5)
} )
--
William R. Shadish
Distinguished Professor
Founding Faculty
Mailing Address:
William R. Shadish
University of California
School of Social Sciences, Humanities and Arts
5200 North Lake Rd
Merced CA 95343
Physical/Delivery Address:
University of California Merced
ATTN: William Shadish
School of Social Sciences, Humanities and Arts
Facilities Services Building A
5200 North Lake Rd.
Merced, CA 95343
209-228-4372 voice
209-228-4007 fax (communal fax: be sure to include cover sheet)
wshad...@ucmerced.edu
http://faculty.ucmerced.edu/wshadish/index.htm
http://psychology.ucmerced.edu
__
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
Re: [R] two questions about xyplot
Dear Richard, your solution to the second question worked like a charm.
Thanks! So much to learn about this stuff, but at least it is fun.
On the first question, yes, I want a text to display the mean in each of
the 12 panels.
Will
On 9/23/2013 11:23 AM, Richard Kwock wrote:
Hi,
To answer your second question you can do something like this:
p-xyplot(dvy ~ sessidx | case, group = numph, data=d66df, col = c(1:4),
layout=c(1, 3), xlab= Sessions,
ylab = Number of Seconds,
type=l)
update(p, panel=function(...){
panel.xyplot(...)
panel.abline(v=6.5)
panel.abline(v=12.5)
panel.abline(v=18.5)
} )
By setting the group parameter in xyplot to be numph, xyplot will
plot different lines for each group of numph you have in each case.
For your first question, did you mean you want a text to display the
mean in each panel?
Richard
On Mon, Sep 23, 2013 at 10:41 AM, William Shadish wshad...@ucmerced.edu wrote:
Dear R helpers,
I am generating three artificial short interrupted time series datasets
(single-case designs; call them Case 1, Case 2, Case 3) and then plotting
them in xyplot. I will put the entire code below so you can reproduce. I
have been unable to figure out how to do two things.
1. Each time series has 24 time points divided into four phases. Call them
phases A, B, C, D for convenience. I have used running() to compute the
means of the observations in each of these four parts; and I saved these as
objects called mn1 (for Case 1), mn2 (Case 2) and mn3 (Case 3). So mn1
contains the four means for A, B, C, D phases for Case 1, etc. I want to
insert these means into the xyplot in the appropriate place. For instance,
insert the first mean from mn1 into phase A of Case 1, the second mean into
phase B of Case 1, and so forth until insert the fourth mean from mn3 into
phase D of Case 3. Ideally, it would insert something like M = 49.02 or
Xbar = 49.02 into phase A for Case 1.
2. The xyplot code I use creates a line connecting the data points, and that
line is continuous over the entire graph. I would like to have the lines be
discontinuous between phases. Phase changes are indicated by panel.abline()
in the code, and occur at time points 6.5, 12.5, and 18.5. So, for example,
I would like a line connecting the datapoints from 1 to 6, then 7-12, then
13-18, then 19-24 (but not including 6-7, 12-13, and 18-19).
I appreciate any help you might be able to offer.
Will Shadish
Here is the code:
#
library(gtools)
library(lattice)
###g = .66
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
###change mean = to vary the effect size
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(1,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d1 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn1 - running(dvy, width=6, by=6)
#second dataset
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(2,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d2 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn2 - running(dvy, width=6, by=6)
#third dataset
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(3,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d3 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn3 - running(dvy, width=6, by=6)
#concatenate d1 d2 d3
d66 - rbind(d1, d2, d3)
d66df - as.data.frame(d66)
d66df$case - ordered(d66df$pid,
levels = c(1,2,3),
labels = c(Case 3, Case 2, Case 1))
p-xyplot(dvy ~ sessidx | case, data=d66df,
layout=c(1, 3), xlab= Sessions,
ylab = Number of Seconds,
type=l)
update(p, panel=function(...){
panel.xyplot(...)
panel.abline(v=6.5)
panel.abline(v=12.5)
panel.abline(v=18.5)
} )
--
William R. Shadish
Distinguished Professor
Founding Faculty
Mailing Address:
William R. Shadish
University of California
School of Social Sciences, Humanities and Arts
5200 North Lake Rd
Merced CA 95343
Physical/Delivery Address:
University of California Merced
ATTN: William Shadish
Re: [R] two questions about xyplot
Richard,
This worked perfectly (adding # before update). Thank you so much for
your help. I've bought a couple of books on R Graphics so I can learn
this stuff better.
Will
On 9/23/2013 1:08 PM, Richard Kwock wrote:
Hi,
Getting text to show on the panel plots is a bit trickier, but doable.
# append to the dataset the mean for each group and line
d66df_mns - cbind(d66df, Means = c(rep(c(mn1, mn2, mn3), each = 6)))
# set the y_lim to extend a bit further above the graph to allow for
the means to be displayed
p-xyplot(dvy ~ sessidx | case, group = numph, data=d66df_mns, col = c(1:4),
layout=c(1, 3), xlab= Sessions,
ylab = Number of Seconds, ylim = c(min(d66df_mns$dvy), 110),
type=l)
# pass in the means as an argument to the panel function
update(p, panel=function(x, y, means = d66df_mns$Means, ... ){
# print(list(...))
# this will store the groups index by subscript into a variable
grps - list(...)$groups[list(...)$subscript]
unique_indices - !duplicated(grps)
# this will get the mean for each panel and for each line
mean_1 - (means[list(...)$subscript][unique_indices])
print(mean_1)
print(x[unique_indices])
panel.xyplot(x, y, ... )
panel.abline(v=6.5)
panel.abline(v=12.5)
panel.abline(v=18.5)
panel.abline(v=18.5)
# print the mean values here.
panel.text(x[unique_indices], 100, paste(M = , round(mean_1,
2)), adj = c(0,0))
} )
If you are working in lattice a lot, print(list(...)) is a handy
function that will show you what parameters values you are passing in
as ... in the panel function.
Hope that helps.
Richard
On Mon, Sep 23, 2013 at 11:23 AM, Richard Kwock richardkw...@gmail.com wrote:
Hi,
To answer your second question you can do something like this:
p-xyplot(dvy ~ sessidx | case, group = numph, data=d66df, col = c(1:4),
layout=c(1, 3), xlab= Sessions,
ylab = Number of Seconds,
type=l)
update(p, panel=function(...){
panel.xyplot(...)
panel.abline(v=6.5)
panel.abline(v=12.5)
panel.abline(v=18.5)
} )
By setting the group parameter in xyplot to be numph, xyplot will
plot different lines for each group of numph you have in each case.
For your first question, did you mean you want a text to display the
mean in each panel?
Richard
On Mon, Sep 23, 2013 at 10:41 AM, William Shadish wshad...@ucmerced.edu wrote:
Dear R helpers,
I am generating three artificial short interrupted time series datasets
(single-case designs; call them Case 1, Case 2, Case 3) and then plotting
them in xyplot. I will put the entire code below so you can reproduce. I
have been unable to figure out how to do two things.
1. Each time series has 24 time points divided into four phases. Call them
phases A, B, C, D for convenience. I have used running() to compute the
means of the observations in each of these four parts; and I saved these as
objects called mn1 (for Case 1), mn2 (Case 2) and mn3 (Case 3). So mn1
contains the four means for A, B, C, D phases for Case 1, etc. I want to
insert these means into the xyplot in the appropriate place. For instance,
insert the first mean from mn1 into phase A of Case 1, the second mean into
phase B of Case 1, and so forth until insert the fourth mean from mn3 into
phase D of Case 3. Ideally, it would insert something like M = 49.02 or
Xbar = 49.02 into phase A for Case 1.
2. The xyplot code I use creates a line connecting the data points, and that
line is continuous over the entire graph. I would like to have the lines be
discontinuous between phases. Phase changes are indicated by panel.abline()
in the code, and occur at time points 6.5, 12.5, and 18.5. So, for example,
I would like a line connecting the datapoints from 1 to 6, then 7-12, then
13-18, then 19-24 (but not including 6-7, 12-13, and 18-19).
I appreciate any help you might be able to offer.
Will Shadish
Here is the code:
#
library(gtools)
library(lattice)
###g = .66
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
###change mean = to vary the effect size
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(1,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir - rep(0,24)
sessidx - c(1:24)
m - rep(1,6)
n - rep(2,6)
o - rep(3,6)
p - rep(4,6)
numph - c(m,n,o,p)
phasebtm - c(b,c,b,c)
d1 - cbind(jid,sid,pid,dvid,desvar,dvdir,dvy,sessidx,numph,phasebtm)
mn1 - running(dvy, width=6, by=6)
#second dataset
z - rnorm(24, mean = 0, sd = 10)
w - rnorm(24, mean = 0, sd = 10)
tm - rnorm(6, mean = 10, sd = 10)
b - rep(0,6)
c - rep(1,6)
tmt - c(b,tm,b,tm)
for (t in 2:24) z[t] - 0.25 * z[t - 1] + w[t]
dvy - 50 + z + tmt
jid - rep(1,24)
sid - rep(1,24)
pid - rep(2,24)
dvid - rep(1,24)
desvar - rep(1,24)
dvdir
Re: [R] gamm in mgcv random effect significance
Gavin et al., Thanks so much for the help. Unfortunately, the command anova(g1$lme, g2$lme) gives Error in eval(expr, envir, enclos) : object 'fixed' not found and for bam (which is the one that can use a known ar1 term), the error is AR1 parameter rho unused with generalized model Apparently it cannot run for binomial distributions, and presumably also Poisson. I did find a Frequently Asked Questions for package mgcv page that said How can I compare gamm models? In the identity link normal errors case, then AIC and hypotheis testing based methods are fine. Otherwise it is best to work out a strategy based on the summary.gam So putting all this together, I take it that my binomial example will not support a direct model comparison to test the significance of the random effects. I'm guessing the best strategy based on the summary.gam is probably just to compare fit indices like Log Likelihoods. If anyone has any other suggestions, though, please do let me know. Thanks so much. Will Shadish On 6/7/2013 3:02 PM, Gavin Simpson wrote: On Fri, 2013-06-07 at 13:12 -0700, William Shadish wrote: Dear R-helpers, I'd like to understand how to test the statistical significance of a random effect in gamm. I am using gamm because I want to test a model with an AR(1) error structure, and it is my understanding neither gam nor gamm4 will do the latter. gamm4() can't yes and out of the box mgcv::gam can't either but see ?magic for an example of correlated errors and how the fits can be manipulated to take the AR(1) (or any structure really as far as I can tell) into account. You might like to look at mgcv::bam() which allows an known AR(1) term but do check that it does what you think; with a random effect spline I'm not at all certain that it will nest the AR(1) in the random effect level. snip / Consider, for example, two models, both with AR(1) but one allowing a random effect on xc: g1 - gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, correlation=corAR1()) g2 - gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random = list(xc=~1),correlation=corAR1()) Shouldn't you specify how the AR(1) is nested in the hierarchy here, i.e. AR(1) within xc? maybe I'm not following your data structure correctly. I include the output for g1 and g2 below, but the question is how to test the significance of the random effect on xc. I considered a test comparing the Log-Likelihoods, but have no idea what the degrees of freedom would be given that s(xc) is smoothed. I also tried: anova(g1$gam, g2$gam) gamm() fits via the lme() function of package nlme. To do what you want, you need the anova() method for objects of class lme, e.g. anova(g1$lme, g2$lme) Then I think you should check if the fits were done via REML and also be aware of the issue of testing wether a variance term is 0. that did not seem to return anything useful for this question. A related question is how to test the significance of adding a second random effect to a model that already has a random effect, such as: g3 - gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random = list(Case=~1, z=~1),correlation=corAR1()) g4 - gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random = list(Case=~1, z=~1, int=~1),correlation=corAR1()) Again, I think you need anova() on the $lme components. HTH G Any help would be appreciated. Thanks. Will Shadish g1 $lme Linear mixed-effects model fit by maximum likelihood Data: data Log-likelihood: -437.696 Fixed: fixed X(Intercept) Xz XintXs(xc)Fx1 0.6738466 -2.56883170.0137415 -0.1801294 Random effects: Formula: ~Xr - 1 | g Structure: pdIdnot Xr1 Xr2 Xr3 Xr4 Xr5 Xr6 Xr7 Xr8 Residual StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 1.693177 Correlation Structure: AR(1) Formula: ~1 | g Parameter estimate(s): Phi 0.3110725 Variance function: Structure: fixed weights Formula: ~invwt Number of Observations: 264 Number of Groups: 1 $gam Family: binomial Link function: logit Formula: y ~ s(xc) + z + int Estimated degrees of freedom: 1 total = 4 attr(,class) [1] gamm list g2 $lme Linear mixed-effects model fit by maximum likelihood Data: data Log-likelihood: -443.9495 Fixed: fixed X(Intercept) Xz XintXs(xc)Fx1 0.720018143 -2.562155820 0.003457463 -0.045821030 Random effects: Formula: ~Xr - 1 | g Structure: pdIdnot Xr1 Xr2 Xr3 Xr4 Xr5 Xr6 Xr7 Xr8 StdDev: 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 Formula: ~1 | xc %in% g (Intercept) Residual StdDev: 6.277279e-05 1.683007 Correlation Structure: AR(1) Formula: ~1
[R] gamm in mgcv random effect significance
Dear R-helpers, I'd like to understand how to test the statistical significance of a random effect in gamm. I am using gamm because I want to test a model with an AR(1) error structure, and it is my understanding neither gam nor gamm4 will do the latter. The data set includes nine short interrupted time series (single case designs in education, sometimes called N-of-1 trials in medicine) from one study. They report a proportion as outcome (y), computed from a behavior either observed or not out of 10 trials per time point. Hence I use binomial (I believe quasi-binomial is not available in gamm). Each of the nine series has an average of 30 observations give or take (total 264 observations), some under treatment (z) and some not. xc is centered session number, int is the z*xc interaction. Based on prior work, xc is also smoothed Consider, for example, two models, both with AR(1) but one allowing a random effect on xc: g1 - gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, correlation=corAR1()) g2 - gamm(y ~ s(xc) +z+ int,family=binomial, weights=trial, random = list(xc=~1),correlation=corAR1()) I include the output for g1 and g2 below, but the question is how to test the significance of the random effect on xc. I considered a test comparing the Log-Likelihoods, but have no idea what the degrees of freedom would be given that s(xc) is smoothed. I also tried: anova(g1$gam, g2$gam) that did not seem to return anything useful for this question. A related question is how to test the significance of adding a second random effect to a model that already has a random effect, such as: g3 - gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random = list(Case=~1, z=~1),correlation=corAR1()) g4 - gamm(y ~ xc +z+ s(int),family=binomial, weights=trial, random = list(Case=~1, z=~1, int=~1),correlation=corAR1()) Any help would be appreciated. Thanks. Will Shadish g1 $lme Linear mixed-effects model fit by maximum likelihood Data: data Log-likelihood: -437.696 Fixed: fixed X(Intercept) Xz XintXs(xc)Fx1 0.6738466 -2.56883170.0137415 -0.1801294 Random effects: Formula: ~Xr - 1 | g Structure: pdIdnot Xr1 Xr2 Xr3 Xr4 Xr5 Xr6 Xr7 Xr8 Residual StdDev: 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 0.0004377781 1.693177 Correlation Structure: AR(1) Formula: ~1 | g Parameter estimate(s): Phi 0.3110725 Variance function: Structure: fixed weights Formula: ~invwt Number of Observations: 264 Number of Groups: 1 $gam Family: binomial Link function: logit Formula: y ~ s(xc) + z + int Estimated degrees of freedom: 1 total = 4 attr(,class) [1] gamm list g2 $lme Linear mixed-effects model fit by maximum likelihood Data: data Log-likelihood: -443.9495 Fixed: fixed X(Intercept) Xz XintXs(xc)Fx1 0.720018143 -2.562155820 0.003457463 -0.045821030 Random effects: Formula: ~Xr - 1 | g Structure: pdIdnot Xr1 Xr2 Xr3 Xr4 Xr5 Xr6 Xr7 Xr8 StdDev: 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 7.056078e-06 Formula: ~1 | xc %in% g (Intercept) Residual StdDev: 6.277279e-05 1.683007 Correlation Structure: AR(1) Formula: ~1 | g/xc Parameter estimate(s): Phi 0.1809409 Variance function: Structure: fixed weights Formula: ~invwt Number of Observations: 264 Number of Groups: g xc %in% g 134 $gam Family: binomial Link function: logit Formula: y ~ s(xc) + z + int Estimated degrees of freedom: 1 total = 4 attr(,class) [1] gamm list -- William R. Shadish Distinguished Professor Founding Faculty Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4007 fax (communal fax: be sure to include cover sheet) wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] te( ) interactions and AIC model selection with GAM
Simon, could you clarify this paragraph. We are submitting to a psychology journal and I am certain they will be asking us about why it is still ok to use AIC for non-nested models. Thanks. Will Shadish On 8/2/2012 9:49 AM, Simon Wood wrote: For AIC model comparison, the usual advice applies that although the comparisons are a bit more reliable for nested models (since then some of the AIC approximations errors are the same in both cases, and cancel in the comparison), it is still ok to use AIC for non-nested models (i.e. it is *not* like most hypothesis testing where everything falls apart if you don't have nesting). best, Simon -- William R. Shadish Distinguished Professor Founding Faculty Chair, Psychological Sciences Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4390 fax wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] te( ) interactions and AIC model selection with GAM
Bert, thanks for the very useful link. I am duly chastised for not having searched thoroughly enough. Will On 8/2/2012 10:02 AM, Bert Gunter wrote: Well, geez! Without trying to upstage SImon, why not google it yourself?! http://en.wikipedia.org/wiki/Akaike_information_criterion (seemed pretty clear to me...) -- Bert On Thu, Aug 2, 2012 at 9:54 AM, Will Shadish wshad...@ucmerced.edu wrote: Simon, could you clarify this paragraph. We are submitting to a psychology journal and I am certain they will be asking us about why it is still ok to use AIC for non-nested models. Thanks. Will Shadish On 8/2/2012 9:49 AM, Simon Wood wrote: For AIC model comparison, the usual advice applies that although the comparisons are a bit more reliable for nested models (since then some of the AIC approximations errors are the same in both cases, and cancel in the comparison), it is still ok to use AIC for non-nested models (i.e. it is *not* like most hypothesis testing where everything falls apart if you don't have nesting). best, Simon -- William R. Shadish Distinguished Professor Founding Faculty Chair, Psychological Sciences Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4390 fax wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- William R. Shadish Distinguished Professor Founding Faculty Chair, Psychological Sciences Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4390 fax wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
Re: [R] curve comparison
I'm very much a novice with R, and not a statistician. But as my group has been working with generalized additive models (semi-parametric regression), we have followed Wood's advice about using the R anova function to do model comparison for different regressions. I would imagine at least some of yours might be nested, e.g., adding polynomial functions to some? If so, model comparisons like this might be justified. Will Shadish On 7/29/2012 10:38 PM, Luis Fernando García Hernández wrote: Dear R users, I have seven regression lines I´d like to compare, in order to find out if these are significatively different. The main problem is that these are curves, non normal, non homogeneous data, I´ve tried to linearize them but it has not worked. So I´d like to know if you know any command or source in R which explains how to perform this kind of comparison. Thanks in advance for your help! [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. -- William R. Shadish Distinguished Professor Founding Faculty Chair, Psychological Sciences Mailing Address: William R. Shadish University of California School of Social Sciences, Humanities and Arts 5200 North Lake Rd Merced CA 95343 Physical/Delivery Address: University of California Merced ATTN: William Shadish School of Social Sciences, Humanities and Arts Facilities Services Building A 5200 North Lake Rd. Merced, CA 95343 209-228-4372 voice 209-228-4390 fax wshad...@ucmerced.edu http://faculty.ucmerced.edu/wshadish/index.htm http://psychology.ucmerced.edu [[alternative HTML version deleted]] __ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
