[R-sig-ME] Should I use full models when using Powersim?

[Alday, Phillip]

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