Should I use full models when using Powersim?
Dear Chi, Does your maximal model converge for your observed data? Are you able to detect all assumed "true" effects (i.e. are all of your effects of interest significant)? If the answer is no: then your data are too noisy to give you a reliable estimate of the effect and this can present itself as weird results in simulation-based power analyses. Basically, if your estimates are noisy, then your simulation will be noisy and your power estimates won't be reliable. This also holds for power analysis without simulation, but won't be as obvious: if you don't have a good estimate of your effect size, then your power analysis will be misleading. As far as computational time: use the model in your power analysis that you want to use in your final analysis. After all, you want to estimate how much power your final analysis has! There is a large amount of literature debating tradeoffs in Type I & II error when using different random-effects structures, and you haven't revealed anything else about your data and experimental design, so there's little specific advice I can give. However, in my experience, removing the interaction terms from the random effects generally speeds up the computation a lot while not really changing model fit. Since the interaction term is the effect you care about, I would reparameterize the model so that interaction is a main effect, e.g. data$ab <- interaction(data$a, data$b) fit <- glmer(B ~ a * ab + (1 + a + ab | Subject) .... If your original model take an hour or more to compute,then it's no surprise that 5000 simulations * 4 points on the power curve = 20 000 model fits take weeks! Best, Phillip PS: Don't name your dataframe "data"! There is a built-in function with that name in R and if you're not careful, you'll get all sorts of weird erros.
On 16/2/19 2:52 pm, Chi Zhang wrote:
Dear all,
I tried using powersim from R package simr to estimate the number of participants that I need for an experiment. I performed the simulation based on the data from my pilot study. The model I used is sketched below:
fit <- glmer(B ~ a+b+a:b
(1+a+b+a:b|Subject) +
(1+a+b+a:b|Item),
family = binomial(link="logit"),
data = data,
control = glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=50000),
tol = .0001))
in which Subject and Item mean the distinct id of subjects and items from the pilot study. I want to know test how the power of the interaction term (a:b) changes with the growth of the number of participants. The code I am using is:
fit2<- extend(fit, along="Subject", n = 84)
sim <- powerCurve(fit2, test = fcompare(~a+b), along = "Subject", breaks=c(48,60,72,84), nsim = 5000)
print(sim)
But the results of the simulation was rather bizarre. To begin with, the power of the interaction grew smaller when the number of participants increased from 72 to 84, which I believe is incompatible with the normal observation that the power increases with the number of participants. Second, I tried using the full random model to perform the simulation, but it is really slow (it took me weeks to get just one result). I was wondering if I can use a simpler random model to perform the simulation.
To reiterate my question: first, why my simulated power decreased with the increase of the number of participants? Is there something wrong with my code? Second, can I use a simpler random model for the simulation in order to save time? Thanks in advance!
Best,
Chi Zhang
------------------------
Chi Zhang
PhD Student
Department of Experimental Psychology
Ghent University
Henri Dunantlaan 2, B-9000 Gent, Belgium
Tel: +32 465386530
E-mail: chi.zhang at ugent.be
________________________________________
From: R-sig-mixed-models <r-sig-mixed-models-bounces at r-project.org> on behalf of r-sig-mixed-models-request at r-project.org <r-sig-mixed-models-request at r-project.org>
Sent: Saturday, February 16, 2019 12:00
To: r-sig-mixed-models at r-project.org
Subject: R-sig-mixed-models Digest, Vol 146, Issue 11
Send R-sig-mixed-models mailing list submissions to
r-sig-mixed-models at r-project.org
To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
or, via email, send a message with subject or body 'help' to
r-sig-mixed-models-request at r-project.org
You can reach the person managing the list at
r-sig-mixed-models-owner at r-project.org
When replying, please edit your Subject line so it is more specific
than "Re: Contents of R-sig-mixed-models digest..."
Today's Topics:
1. gls() funcion after R upgrade (Fabio Fanoni)
2. Help with: Multivariate models with binary variables and
variables with zipoisson distribution (Stacey Hannebaum)
3. Teaching Assistants for advanced multilevel modeling at ICPSR
(Poe, John)
----------------------------------------------------------------------
Message: 1
Date: Mon, 11 Feb 2019 08:40:28 +0000
From: Fabio Fanoni <fabio.fanoni at popso.it>
To: "r-sig-mixed-models at r-project.org"
<r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] gls() funcion after R upgrade
Message-ID: <d0cd3fe9653e44dcbdb782bb8af4094a at popso.it>
Content-Type: text/plain; charset="utf-8"
Dear all,
I have been using the package nlme for about a year through the software R
Recently, following an update of R (and the corresponding nlme package), in some situations the gls() function does not work.
The following lines of code
require(nlme)
db<-read.csv2("G:\\Dati.csv",as.is=TRUE)
f<-"DELTA_INVNORM_STD ~ 0+ ITOD + EU_CH"
gls(as.formula(f),correlation=corARMA(p=2,q=3,fixed=FALSE),data=db,weights=varExp(),method="ML")
with the old version of R returns the following output:
Generalized least squares fit by maximum likelihood
Model: as.formula(f)
Data: db
Log-likelihood: -36.13853
Coefficients:
ITOD EU_CH
-0.3138869 -0.5010393
Correlation Structure: ARMA(2,3)
Formula: ~1
Parameter estimate(s):
Phi1 Phi2 Theta1 Theta2 Theta3
0.50054107 -0.07357577 0.91847373 0.66706344 0.51689285
Variance function:
Structure: Power of variance covariate
Formula: ~fitted(.)
Parameter estimates:
power
-0.01893321
Degrees of freedom: 71 total; 69 residual
Residual standard error: 0.4077908
while with the new one version of R ends with the following error message
Error in `coef<-.corARMA`(`*tmp*`, value = value[parMap[, i]]) :
NA/NaN/Inf in foreign function call (arg 1)
here are the parameters of the old version of R
sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows >= 8 x64 (build 9200)
locale:
[1] LC_COLLATE=Italian_Italy.1252 LC_CTYPE=Italian_Italy.1252 LC_MONETARY=Italian_Italy.1252
[4] LC_NUMERIC=C LC_TIME=Italian_Italy.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] nlme_3.1-128
loaded via a namespace (and not attached):
[1] tools_3.3.2 grid_3.3.2 lattice_0.20-34
and the correspondents of the new one
sessionInfo()
R version 3.5.2 (2018-12-20)
Platform: i386-w64-mingw32/i386 (32-bit)
Running under: Windows >= 8 x64 (build 9200)
Matrix products: default
locale:
[1] LC_COLLATE=Italian_Italy.1252 LC_CTYPE=Italian_Italy.1252 LC_MONETARY=Italian_Italy.1252
[4] LC_NUMERIC=C LC_TIME=Italian_Italy.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] nlme_3.1-137
loaded via a namespace (and not attached):
[1] compiler_3.5.2 tools_3.5.2 grid_3.5.2 lattice_0.20-38
The data used are shown below
I wonder if someone could help me to solve this problem.
Best regards,
Fabio Fanoni
DELTA_INVNORM_STD;ITOD;EU_CH
-1,01888891391635;1,11744160507933;0,164921569515933
-1,7287582158078;1,33122211515022;0,192399322297369
-1,92084919508607;1,42620153798496;0,119394703266418
-1,07481132504896;1,40851069525448;-0,21451386962765
0,368325144516397;1,143582080184;-0,541418306765893
1,60232300735139;0,877753890833926;-0,540071359447086
1,79762928248747;0,585114281768911;-0,457099318942592
0,91995670533225;0,300270733467356;-0,292232806294851
-0,0262660321428139;0,20146037523289;-0,236469135573377
-1,02121152912805;0,0419782544708825;-0,345437279939175
-1,10672650470506;0,0574913198102401;-0,366584376685452
-1,39178681241528;0,33582796272111;-0,222730255141816
-1,1454072922613;0,119188584353956;-0,00950828000562807
-0,958672816825354;-0,0667970836192611;0,4570747351185
-1,32272038017045;-0,307712992787705;0,958678420940093
-0,907400017348963;-0,781894397724971;1,3426934842003
-1,24178782808971;-0,502123817789113;1,78543551198139
-0,727279530248652;-0,198729005124588;1,65424271463061
-0,663683090224516;-0,00292614483254617;1,29016249027665
-0,870662613154294;0,284152373447782;0,852269477692287
-0,740956260890321;0,212558271927232;0,359286261984589
-0,875175106507846;0,16921613723138;0,301636857743261
-0,159691021207123;0,344902273564796;0,406160071512065
0,485359221230746;0,405193753056144;0,545299871378597
0,91806657022351;0,51772409209882;0,664639525858513
0,94182403298868;0,624398008842597;0,721750147939642
0,772810225693767;0,669380547864261;0,746938090683181
0,454079892235428;0,829860577818214;0,875841077996076
0,216198540864624;0,85936757871922;1,06778126286636
-0,390013295840405;0,820301653758161;1,32356681853436
-1,08957286663887;0,735286971796754;1,51483352772482
-1,7181964072478;0,54188741408352;1,59807495122772
-2,14142427321141;0,3080114833751;1,30188093962788
-1,73631132303729;0,242276386254216;0,821154955431617
-1,14832298161027;-0,168559611137128;0,390400772122712
-0,199657255063557;-0,673990227176314;-0,425580724204736
0,637711021315926;-1,26311924357155;-0,777807797606142
1,13372564859946;-2,22697608948676;-1,02793616079627
1,3387228312938;-2,67326030450369;-1,25476231880355
1,13410628583832;-2,41425486942835;-0,769187327492251
1,20371911407773;-1,73317141026889;-0,4868668898436
1,04245681382429;-0,627677467262121;-0,517038541303494
1,19122635858486;0,318264079861633;-0,902669938921763
1,10644115835792;0,53399823108225;-1,58826681113873
0,657027862693417;0,566924775110004;-2,16826290930735
0,391546504792817;0,382793958142103;-2,37380727704873
-0,121513593301095;0,113508212335312;-2,28814134594747
0,090587410027589;-0,10456301293743;-1,92150191643467
0,293363258062962;-0,599038773090656;-1,5279235122598
0,490646403474033;-1,16804784867814;-1,10646327070376
0,663013157660536;-1,68022387383074;-0,281457211581005
0,48006697458688;-2,17535440458421;0,0140633243390067
0,33992207153128;-2,15751594621565;0,415319331867849
0,485131796838891;-1,91164287702336;0,74289724468591
0,645057199360022;-1,64072137764755;0,715150100824384
0,830719890823037;-1,20839230550334;0,886616664628225
1,00060688160202;-0,799350150825283;0,78775063121476
0,799516770753564;-0,474799224784852;0,620190215847712
0,884527141039385;-0,111321336688718;0,400772282236831
0,77592825475446;0,41873781137028;0,214623979323192
0,897761905934407;0,618216183184323;-0,371702774880583
1,23087946317063;0,86282014652723;-1,04275259746654
1,1730433788867;1,03592979384644;-1,50933562067236
1,18413315011977;0,843466359597336;-1,89254251877396
0,522306442129412;0,786432900675309;-1,199941514842
-0,423472231576842;0,534354713733984;-0,258289703958632
-1,02160314495135;0,445296924001273;0,367637345059225
-1,2194743130469;0,791437804698057;0,831122380946011
-0,682787491324478;0,735393577499657;0,636488297397258
-0,277058822726022;0,988210808618591;0,362382895620209
-0,0283252883031908;1,18581655608866;0,466771415869385
Fabio Fanoni
SVILUPPO MODELLI DI CREDITO
Piazza Garibaldi 16
23100 Sondrio (SO)
Banca Popolare di Sondrio
Tel. +390342528920
fabio.fanoni at popso.it<mailto:fabio.fanoni at popso.it>
AVVERTENZE LEGALI - I contenuti di questo messaggio, proveniente da un indirizzo di posta elettronica aziendale della Banca Popolare di Sondrio, e gli eventuali allegati possono essere letti e utilizzati, per esigenze lavorative, da chi opera alle dipendenze o per conto della stessa, o ? comunque cointeressato nell'inerente relazione d'affari. Le dichiarazioni, ivi contenute, non impegnano contrattualmente la Banca Popolare di Sondrio se non nei limiti di quanto eventualmente previsto in accordi opportunamente formalizzati. Se il messaggio ? stato ricevuto per errore ce ne scusiamo, pregando di segnalare ci? al mittente e poi di distruggerlo senza farne alcun uso poich? l'utilizzo senza averne diritto ? vietato dalla legge e potrebbe costituire reato.
DATI SOCIETARI - BANCA POPOLARE DI SONDRIO - Societ? cooperativa per azioni - Fondata nel 1871
Sede sociale e direzione generale: I - 23100 SONDRIO SO - piazza Garibaldi, 16
Indirizzo Internet: http:[doppiabarra]www[punto]popso[punto]it - E-mail: info[chiocciola]popso[punto]it
Iscritta al Registro delle Imprese di Sondrio al n. 00053810149, all'Albo delle Banche al n. 842, all'Albo delle Societ? Cooperative al n. A160536 Capogruppo del Gruppo bancario Banca Popolare di Sondrio, iscritto all'Albo dei Gruppi bancari al n. 5696.0 Aderente al Fondo Interbancario di Tutela dei Depositi e al Fondo Nazionale di Garanzia - Codice fiscale e partita IVA: 00053810149
Capitale sociale: EUR 1.360.157.331 - Riserve: EUR 947.325.264 (Dati approvati dall'Assemblea dei soci del 29 aprile 2017).
N.B. I "filtri antivirus" e "antispam" in uso su molti sistemi di posta elettronica possono talvolta ritardare o impedire, in tutto o in parte, il recapito dei messaggi. In tali casi, salvo verifica di avvenuta ricezione, pu? essere necessario modificare i contenuti o le modalit? d'invio.
[[alternative HTML version deleted]]
------------------------------
Message: 2
Date: Fri, 15 Feb 2019 14:25:09 -0600
From: Stacey Hannebaum <stacey-hannebaum at utulsa.edu>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] Help with: Multivariate models with binary
variables and variables with zipoisson distribution
Message-ID:
<CAFtRfJXA5jOSX0n343Zs_wcyXnC1Dy3iw06GV513Li7GxUx8nw at mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
Hi,
I am trying to build a multivariate model to test for associations between
individual variation in behavioral and other traits. Essentially, I want to
derive two calculations from the model: (1) the repeatability of each of
the three traits within individuals and (2) the correlation between each
pair of traits and between each trait and fitness and trappability (the
latter two have one value per individual).
I have been working on this for some time and have reached a standstill.
Here are the main issues I am facing. Any suggestions would be very
appreciated.
(1) How should a binary variable be modelled? For now, I have modelled it
as Gaussian (with the assumption that 0s and 1s are drawn from some
underlying Gaussian process) so that I can calculate correlations but I
don?t know if this is reasonable. Any thoughts?
(2) I tried to run the below model but received the following error:
Error in priorformat(if (NOpriorG) { :
V is the wrong dimension for some prior$G/prior$R elements
I am sure it has to do with modelling CHARGES as zipoisson - I think I need
to specify the structure for modelling the 0?s but I don?t know where to
start with this.
Below are the variables of interest in case this is useful to know:
logitTAPERETURN: logit-transformed latency to enter a nest once a novel
stimulus was applied / total possible time allowed for trial; trials were
terminated after 5 minutes whether or not bird entered nest; multiple
trials per individual.
CHARGES: number of times an individual attacked the novel stimulus (count
data with many zeros); multiple trials per individual. This could be
expressed as a rate (number of charges / latency to enter nest) but I
decided to model length of trial as a predictor variable using log(latency
to enter nest); trial repeated multiple times per subject
STATUSAFTER: whether the bird stayed in the nest or fled from the nest
(binary)
CAPTURES: number of times an individual was captured (count data with many
ones and no zeros); only one value per individual
RELFITNESS: number of fledglings for individual / mean number of fledglings
across all individuals; only one value per individual
prior_tr_c_sa_nc_f =
list(G=list(G1=list(V=diag(c(1,1,1,0.0001,0.0001),5,5),nu=1.002, fix = 4)),
R=list(V=diag(5),nu=1.002))
I fixed CAPTURES and RELFITNESS because there is not within-individual
variation in these variables.
mcmc_tr_c_sa_nc_f_us <- MCMCglmm(cbind(logitTAPERETURN, CHARGES,
STATUSAFTER, CAPTURES, RELFITNESS) ~
at.level(trait,1):(COLSIZE) +
at.level(trait,1):(TAPETRIALNUM) +
at.level(trait,1):(TAPENESTSTATUS) +
at.level(trait,2):(COLSIZE) +
at.level(trait,2):(stTAPERETURN) +
at.level(trait,2):(stTAPETEMP) +
at.level(trait,3):(stALARMTRIALNUM) +
at.level(trait,3):(stALARMTIME)
+
at.level(trait,3):(stALARMWINDSP) +
at.level(trait,3):(STATUSBEFORE) +
at.level(trait,4):(COLSIZE) +
at.level(trait,5):(COLSIZE) +
at.level(trait,5):(stCLUTCHSIZE) - 1,
random =~ us(trait):BAND,
rcov =~ us(trait):units,
family =
c("gaussian","zipoisson","gaussian","ztpoisson","gaussian"),
prior = prior_tr_c_sa_nc_f,
nitt=420000,
burn=20000,
thin=100,
verbose = TRUE,
data =
as.data.frame(MultivariateAnalysis20172018))
Thank you,
Stacey Hannebaum
Postdoctoral Research Associate
[[alternative HTML version deleted]]
------------------------------
Message: 3
Date: Fri, 15 Feb 2019 17:46:27 -0600
From: "Poe, John" <jdpo223 at g.uky.edu>
To: r-sig-mixed-models <r-sig-mixed-models at r-project.org>
Subject: [R-sig-ME] Teaching Assistants for advanced multilevel
modeling at ICPSR
Message-ID:
<CAFW8BypZZqASQ7zdBCYFS3zBa04DPWtWebHjUVoF4cg0RXgPbg at mail.gmail.com>
Content-Type: text/plain; charset="utf-8"
I'm looking for recommendations on potential teaching assistants for my
4-week advanced multilevel modeling course at the ICPSR summer program in
Ann Arbor, Michigan. If you know of advanced graduate students or postdocs
who might be qualified and interested please let me know off the listserv
or have them email me directly. It's a great program and Ann Arbor is an
absolutely lovely place to spend some time in the summer.
The class pulls from econometrics, biostatistics, and psychometrics and
covers both MLE and Bayesian multilevel modeling. I don't expect anyone to
be comfortable with everything so they'd need to be willing to go through
the material with me in an informal readings bootcamp before the class to
fill in gaps.
You can find last year's syllabus here:
http://www.johndavidpoe.com/teaching/multilevel-modeling-ii-advanced-topics/
[[alternative HTML version deleted]]
------------------------------
Subject: Digest Footer
_______________________________________________
R-sig-mixed-models mailing list
R-sig-mixed-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
------------------------------
End of R-sig-mixed-models Digest, Vol 146, Issue 11
***************************************************
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models