[R-sig-ME] Should I use full models when using Powersim?
Ben Bolker
bbo|ker @end|ng |rom gm@||@com
Mon Feb 18 22:19:32 CET 2019
Agree with everything said here. A few super-low-tech partial
solutions ...
Others may disagree, but 5000 sims seems pretty extreme for estimating
power (which is after all never better than a crude estimate); if you
reduced that 10-fold you'd probably still get an adequate estimate of power.
It doesn't look like simr does parallel computation yet - if it did,
and if you have access to a reasonably powerful machine, you could
reduce the computation by another factor of <number of cores>. (You
could do this by brute force by manually running a bunch of jobs with
different random-number seeds, or using the parallel package's
parLapply(), then averaging the power curves you get ...)
On 2019-02-18 4:02 p.m., Alday, Phillip wrote:
> 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 using ugent.be
>>
>> ________________________________________
>> From: R-sig-mixed-models <r-sig-mixed-models-bounces using r-project.org> on behalf of r-sig-mixed-models-request using r-project.org <r-sig-mixed-models-request using r-project.org>
>> Sent: Saturday, February 16, 2019 12:00
>> To: r-sig-mixed-models using 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 using 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 using r-project.org
>>
>> You can reach the person managing the list at
>> r-sig-mixed-models-owner using 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 using popso.it>
>> To: "r-sig-mixed-models using r-project.org"
>> <r-sig-mixed-models using r-project.org>
>> Subject: [R-sig-ME] gls() funcion after R upgrade
>> Message-ID: <d0cd3fe9653e44dcbdb782bb8af4094a using 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 using popso.it<mailto:fabio.fanoni using 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 using utulsa.edu>
>> To: r-sig-mixed-models using r-project.org
>> Subject: [R-sig-ME] Help with: Multivariate models with binary
>> variables and variables with zipoisson distribution
>> Message-ID:
>> <CAFtRfJXA5jOSX0n343Zs_wcyXnC1Dy3iw06GV513Li7GxUx8nw using 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 using g.uky.edu>
>> To: r-sig-mixed-models <r-sig-mixed-models using r-project.org>
>> Subject: [R-sig-ME] Teaching Assistants for advanced multilevel
>> modeling at ICPSR
>> Message-ID:
>> <CAFW8BypZZqASQ7zdBCYFS3zBa04DPWtWebHjUVoF4cg0RXgPbg using 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 using 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 using r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
More information about the R-sig-mixed-models
mailing list