title: Cohen's d between groups over time question

text: I have a problem trying to properly phrase my results. I utilized: $$D =\frac{(\text{M}_{2e} –\text{M}_{2c}) – (\text{M}_{1e}-\text{M}_{1c})}{\text{SD}_{1\text{pooled}}}$$ for pooled SD I used: $$\text{SD}_{1\text{pooled}} = \sqrt{\left((n_1-1)s_1^2 + (n_2-1)s_2^2\right)\big /(n_1+n_2-2)}$$ My results for my manipulation check was $d = .08$. My results for my hypothesis testing was $d = .38$. Can someone please explain what this means? I know that Cohen's $d$ of .08 is below the .2 criteria of small effect, and the Cohen's $d$ of .38 is a small to moderate effect. Where I am confused is what this means for a between groups over time comparison... Thank you or any insight you may be able to provide me.


title: Multiple correlated random non-nested intercepts in R

text: I am trying to estimate a longitudinal model in R in which there are several random intercepts that are correlated with each other, and the data are non- nested. For example, consider a simple longitudinal model in which student test scores are regressed on three random intercepts, one random school effect, one random teacher effect, and one random teacher-school match effect. The data are non-nested because teachers frequently switch schools in the longitudinal data. I am primarily interested in the random effects parameters themselves, and so I want to specify all three random effects and allow them to be correlated with each other. There are no random slopes, only random intercepts. I tried:               model1 <- lmer(test_score ~ FEs + (1|schoolid) + (1|teacherid) + (1|matchid), data)      but this appears to force the random effects to be uncorrelated with each other. Is there a way, using `lmer` or any other R package to estimate a similar model in which the random intercepts are allowed to be correlated with each other (using non-nested data)? More specifically, the above model imposes the covariance structure: $$ Cov \left[ \begin{array} \\\ \theta_{s} \\\ \theta_{t} \\\ \theta_{m} \\\ \end{array} \right] = \left[ \begin{array} \\\ \sigma^{2}_{\theta_{s}}I_{s} & 0 & 0 \\\ 0& \sigma^{2}_{\theta_{t}}I_{t} & 0\\\ 0&0&\sigma^{2}_{\theta_{m}}I_{m}\\\ \end{array} \right] $$ What I would like to estimate instead is a random effects model with the following covariance structure: $$ Cov \left[ \begin{array} \\\ \theta_{s} \\\ \theta_{t} \\\ \theta_{m} \\\ \end{array} \right] = \left[ \begin{array} \\\ \sigma^{2}_{\theta_{s}}I_{s} & \sigma^{2}_{\theta_{st}}C_{st} & \sigma^{2}_{\theta_{sm}}C_{sm} \\\ \sigma^{2}_{\theta_{ts}}C_{ts} & \sigma^{2}_{\theta_{t}}I_{t} & \sigma^{2}_{\theta_{tm}}C_{tm}\\\ \sigma^{2}_{\theta_{ms}}C_{ms} &\sigma^{2}_{\theta_{mt}}C_{mt} &\sigma^{2}_{\theta_{m}}I_{m}\\\ \end{array} \right] $$ where $C_{st}$, for example, is a matrix of zeros and ones, with $C(i,j)=1$ where there is an observation of the pair $(school=i, teacher=j)$ in the data, and 0 if no such pair is observed.


title: What is the best way to enter data to analyze answers from a Likert scale?

text: I need to analyze data regarding an item based on a Likert scale. Field data was arranged in Excel using tables that show the frequency of answers among respondents. What is the best approach to enter data in the statistical software in order to analyze it? An example of what I have in Excel is attached. I have a list of 32 occupations (rows) that need to be analyzed among 8 different industrial sectors. Sample size = 100. ![enter image description here](http://i.stack.imgur.com/smDrb.png) ![enter image description here](http://i.stack.imgur.com/g3saq.jpg)


title: Hypothesis test of mean between groups with dependencies and pairing

text: I have two groups with small number of n (less than 10) where some observations are replicates for each other, some are paired difference observations, and some are completely independent. For example:               Group I        Group II     A <--paired--> B     C (unpaired)     D <--paired--> E     F <-| replicate of D      Is there any known way to test the hypothesis that the mean of Group I != mean of Group II?


title: Is there a regression method for fitting general N-order polynomials of two or more variables?

text: I would like to fit a general N-order polynomial of two or more variables. For example, consider (for reference) the general second-order polynomial in two variables $x$ and $y$: $$\beta_{00}+\beta_{01}y^1+\beta_{02}y^2+\beta_{10}x+\beta_{11}xy+\beta_{12}xy^2+\beta_{20}x^2+\beta_{21}x^2y+\beta_{22}x^2y^2.$$ Another example could be the general N-order polynomial in two variables $x$ and $y$: $$\sum_{n=0}^{N}\sum_{m=0}^{N}\beta_{nm}x^ny^m$$ I know that there is an extension of linear regression named Multiple Linear Regression, but is there a regression method for fitting general N-order polynomials of two or more variables?


title: jags missing data error

text: I am estimating a probit in jags on a time series cross sectional data set (country years) with 5 covariates, 4 of which are lagged, so the first observation for each country will be missing on those 4 variables. There is some scattered missingness elsewhere in the data, but nothing too bad. All of the variables are of numeric type, and missing values are coded as 'NA'. The model works fine when I omit all missing values, but when I don't I get the error               Compiling model graph     Resolving undeclared variables     Allocating nodes     Deleting model     Error in jags.model("model.bug", data.jags, n.chains = 4, n.adapt = 1000) :      RUNTIME ERROR:     Compilation error on line 4.     Unable to resolve node Xbeta[1]     This may be due to an undefined ancestor node or a directed cycle in the graph      Here is my data loading code               data.all <- read.csv("http://dl.dropbox.com/u/1253073/blowback-all.csv")          data.jags <- list('onset.b.lag' = data.all$onset.b.lag,                       'ongoing' = data.all$ongoing,                       'peace' = data.all$peace,                       'ln.gdppc.lag' = data.all$ln.gdppc.lag,                       'all.lag' = data.all$all.lag,                       'onset' = data.all$onset,                       'n' = nrow(data.all))      rjags call:               mod.all <- jags.model('model.bug', data.jags, n.chains = 4, n.adapt = 1000)      Here is my model code:               model{          for (i in 1:n) {     probit(p[i]) <- Xbeta[i]     onset[i] ~ dbern(p[i])     Xbeta[i] <- b.0 +                 b.onset.b.lag * onset.b.lag[i] +                 b.peace * peace[i] +         b.ln.gdppc.lag * ln.gdppc.lag[i] +         b.ongoing * ongoing[i] +         b.all.lag * all.lag[i]     }          b.0 ~ dnorm(0, .0001)     b.peace ~ dnorm(0, .0001)     b.onset.b.lag ~ dnorm(0, .0001)     b.ln.gdppc.lag ~ dnorm(0, .0001)     b.ongoing ~ dnorm(0, .0001)     b.all.lag ~ dnorm(0, .0001)      } My understanding is that when jags encounters an index where a variable is coded `NA` it takes a draw from the posterior. Any hints would be appreciated.


title: What is the probability of an occurrence within an interval

text: I am having an issue getting my head around this question as well. Gale-force winds occur on the average of 0.7 times every 54 hours. For an interval of 110 hours whats the probability r = 2? I believe I should solve it with Poisson Distribution (as it is in this weeks lesson, and question references using two decimal places for lambda). I would like to know the methodology - not necessarily the answer, but if the answer was supplied it might help me work back and forth through the problem. Thank you in advance.


title: Neural net cost function for Hyperbolic Tangent activation

text: In Andrew Ng's online machine learning course in the part about neural nets for classification the following convex cost function is given: $$\text{cost} = -y\log(h_0(x)) - (1-y)\log(1-h_0(x))$$ which is predicated on the output labels being either 1 or 0 (a sigmoid activation function). I would like to code a NN classifier using a hyperbolic tangent activation function, particularly that given in Le Cun 1998, shown below: $$1.7159\tanh((2/3)x)$$ which will take on values between >1 and < 0 in the limit. Obviously the above cost function would not be suitable for this activation function, so could anyone let me know what sort of similar cost function for this hyperbolic tangent activation I should be looking for? So far my web search has been fruitless.


title: How to calculate the p-value for a binomial test using pbinom?

text: 25 patients receive treatment, and they each get a difference score (After- Before). 18 of them get positive scores. I would like to test the null hypothesis that p=0.5. I know that `1-pbinom(17,25,0.5)` gives me the probability of 18 or more patients getting a positive score. My question is, to calculate the p-value, why do we have to do               2*(1-pbinom(17,25,0.5))      That is, why do we multiply by 2 here?


title: Estimates of regression coefficients are uncorrelated?

text: Consider a simple regression (normality not assumed): $Y_i = a + b X_i + e_i$ where $e_i$ is with mean 0 and standard deviation $\sigma$. Are the Least Square Estimates of $a$ and $b$ uncorrelated?


title: Seemingly unrelated regressions and difference in coefficients across groups

text: Does anyone know a good reference about how to use seemingly unrelated regressions to test for teh difference in coefficients across groups?


title: How to generate random auto correlated binary time series data?

text: How can I generate binary time series such that:   1. Average probability of observing 1 is specified (say 5%);   2. Conditional probability of observing 1 at time $t$ given the value at $t-1$ ( say 30% if $t-1$ value was 1)?


title: Filtering using ARMA model in R

text: I have two time-series, `x` and `y`. I would like to prewhiten `x` by fitting an ARMA(p,q) (or in my case ARMA(1,1)) process and then use the coefficients to filter `y`. This seems like a pretty standard thing to want to do. However, the `stats:::filter` function does only MA or AR filtering it looks like. What is the appropriate way to do this? Also, should one use the `arima` function in R to do this or are there other ways?


title: What is a good reference on simulating ordinal random variables?

text: I'm interested in any references for simulating ordinal random variables to be used in a Bayesian (multilevel) analysis.


title: Internal consistency reliability in item response theory models

text: How can internal consistency reliability of a test and of individual test items be quantified in Item Response Theory models? I know I can resort to Classical Test Theory, Cronbach's alpha, and other measures, but is there a way to characterize reliability within IRT?


title: Identify the parameters of the model $Y=\exp(\beta_0 + \beta_1 X + \beta_2 Z)+u_i$

text: I have the model $Y_i=\exp(\beta_0 + \beta_1 X_i +\beta_2 Z_i) + u_i$ where we assume $\mathbf{E}[u_i|X_i,Z_i]=0$ and $Var(X_i)>0,Var(Z_i)>0$, and I need to show that $\beta_0,\beta_1,\beta_2$ are identified. I know that it's possible to take a log transformation somehow to make it linear but I can't seem to figure out how to get rid of the error term. I'm completely fumbling around here in the dark but I figured the best thing would be to try to write it as a product rather than as a sum, so I said let $v_i=1+\frac{u_{i}}{\exp\left(\beta_{0}+\beta_{1}X_{i}+\beta_{2}Z_{i}\right)}$ so that we can write the model as $$Y_{i}=v_{i}\exp\left(\beta_{0}+\beta_{1}X_{i}+\beta_{2}Z_{i}\right) $$Now I can write $$\ln Y_{i}=\beta_{0}+\beta_{1}X_{i}+\beta_{2}Z_{i} +\ln v_i$$ Now if I could show that $\mathbf{E}[\ln v_i]$ was 0 I would be done. Now, I can see, roughly, that $\mathbf{E}[v_i]$ is maybe 1, if I do some tricky stuff with the law of iterated conditional expectations (I don't know exactly how I would do it but I think I can figure that out) but even if it is, it's not true in general that $\mathbf{E}[\ln v_i]=\ln \mathbf{E}[v_i]$, so that wouldn't actually help me anyway. Can anyone see how I should proceed? Has what I've done so far made sense? * * * **Edited to clarify** Apparently there is more than one way to say that a parameter is "identified". The definition that I am working with is this: a parameter is identified when, if we knew the distribution of observed variables, we could uniquely specify the parameter.


title: How can I estimate unique occurrence counts from a random sampling of data?

text: Let's say I have a large set of $S$ values which sometimes repeat. I wish to estimate the total number of _unique_ values in the large set. If I take a random sample of $T$ values, and determine that it contains $T_u$ unique values, can I use this to estimate the number of unique values in the large set?


title: Probability inequalities

text: I am looking for some probability inequalities for sums of unbounded random variables. I would really appreciate it if anyone can provide me some thoughts. My problem is to find an exponential upper bound over the probability that the sum of unbounded i.i.d. random variables, which are in fact the multiplication of two i.i.d. Gaussian, exceeds some certain value, i.e., $\mathrm{Pr}[ X \geq \epsilon\sigma^2 N] \leq \exp(?)$, where $X = \sum_{i=1}^{N} w_iv_i$, $w_i$ and $v_i$ are generated i.i.d. from $\mathcal{N}(0, \sigma)$. I tried to use the Chernoff bound using moment generating function (MGF), the derived bound is given by: $\begin{eqnarray} \mathrm{Pr}[ X \geq \epsilon\sigma^2 N] &\leq& \min\limits_s \exp(-s\epsilon\sigma^2 N)g_X(s) \\\ &=& \exp\left(-\frac{N}{2}\left(\sqrt{1+4\epsilon^2} -1 + \log(\sqrt{1+4\epsilon^2}-1) - \log(2\epsilon^2)\right)\right) \end{eqnarray}$ where $g_X(s) = \left(\frac{1}{1-\sigma^4 s^2}\right)^{\frac{N}{2}}$ is the MGF of $X$. But the bound is not so tight. The main issue in my problem is that the random variables are unbounded, and unfortunately I can not use the bound of Hoeffding inequality. I will be to happy if you help me find some tight exponential bound . Thanks in advance


title: How can I determine which columns best correlate with target column in a SQL Server Table?

text: I have a table with a couple hundred columns and a few hundred thousand rows and want to determine which combination of 5-10 columns best correlate with a target column. I have SQL Server and Excel but no budget beyond that. My current plan is:   1. Pull the table into Excel   2. Use the correl() function on each column   3. Combine subsets of the best correlated columns to see if their combined correlation coefficient improves. (i.e. new ab merge column =A1*1000+B1) Is there a better way to do this?


title: How do I calculate and compare incidence of infection between age-groups

text: We work on cervical cancer detection and we started accumulating data in 2006 by performing HPV testing on all our samples. We have accumulated a data set with almost 400,000 women. Data enters the database on a daily basis. Most women come for a yearly smear so most women are more than one time in the database (2-13x). Most women are HPV negative (appr 85%) but some become HPV positive during this observation period. I would be very interested in looking at the rate of infection (incidence) of HPV16. I suppose I will have to look at the women with at least 2 samples and restricting to women going from a negative sample to a positive sample. I have a few questions:   1. How do I start? Do I use a survival analysis and how do I do this specifically in Stata or Prism5?   2. Do I censor some women and how (these women who are HPV16+ on the first sample, we don't know when they acquired this)?   3. What do I do with all the HPV negative women, do they stay in the analysis?   4. How would I compare the incidence rate according to age? I would specifically like to know if young women (<30 years) have a different pattern when compared to women of >30 years.


title: Knot selection for cubic regression splines

text: I was wondering if anybody had experience in how to set the knot points when using cubic regression splines. Some background: I have a response and predictor variable, and I want to determine the trend relationship between the two. To see what it looks like without making too many assumptions, I've fit a smoothing spline curve using the `gam` function in R. The trend is obviously not linear, but otherwise well-behaved: smooth, and not too wiggly. I'd now like to model this trend using a simple, cubic regression spline (there are various practical issues with using the `gam` fit, or I'd just use that). Of course, using a regression spline requires the knots to be specified in advance. It's not too hard to do that with linear splines: I'd insert a knot where the slope of the smooth fit changes substantially, eg around local minima/maxima. However, cubic splines appear to be a more complicated story. Any guidance on where I should put the knots would be much appreciated.


title: Statistical significance when A/B test has multiple values

text: I'm sure this is a pretty standard statistics question, but I'm no expert... I'm running an A/B test on my website to see if a change results in users adding more content. So there are 2 basic things I'm looking at; the # of users adding at least 1 piece of content and the total # of pieces of content added by all users. I really care much more about the total # of pieces of content added by all users. I'll make the change permanent if I know it's at least not worse than the existing site. So I need to know how many samples (users logging in) I need to have a 95% confidence level. Normally I can use one of the many web A/B test calculators that use chi-square and similar tests to figure out if my test is statistically significant or to figure out the sample size I need. In the first case of seeing how many users added content, I can do this. But to see the total pieces of content added among all users, I can't use those tests as their isn't a "conversion" event. So what's the best way to see what sample size I need to be statistically significant to a 95% confidence level? And how can I see if my test shows whether there was a difference? Again I just want to make sure the new change isn't worse (or isn't "much" worse, do I need to and how should I define "much"?).


title: Methods for checking if the number of comparisons is sufficient

text: This question is motivated by a issue regarding network motifs. To determine if a (connected, induced) subgraph $H$ occurs with significantly high frequency in an input network $G$, we generate an ensemble of comparison networks similar to $G$ and count the number of occurences of $H$ in them. Thus we obtain the number of copies of $H$ in $G$ (call this frequency $f$), and the number of copies of $H$ in the comparison networks (call these frequencies $x=(x_1,x_2,...,x_n)$). A trade-off arises: We want $n$ to be large in order to give a better comparison, and $n$ to be small to be able to actually perform the computation. I'm concerned, however, that sometimes $n$ is chosen too small. **Question** : How can we determine if $n$ is too small? I.e. how can we determine if we have sufficiently many comparisons? My feeling is to compare $f$ with the mean of datasets $x$ and `c(x,ceiling(mean(x)))` [in R notation]. I.e. in the second case, we add an artificial point to the dataset, which is close to the mean of $x$, integral and non-zero. The idea is, if we instead chose to use $n+14 comparison networks, and the new comparison graph turned out to have `ceiling(mean(x))` copies of $H$, would we have come to a different conclusion? If so, this would indicate that $n$ is too small.


title: Sampling from more than one population

text: In a practice exam I encountered a question about the standard error of a sample (calculating the mean order size of a random sample of 100 customers of a paper firm). A researcher (the subject of this question) advised against increasing the sample size because -although it would decrease the standard error- it also would increase the risk of _sampling from more than one population_. According to the provided answer the researcher was correct. While I understand what problems may be caused by sampling from more than one population, I fail to see how this could be taken into account in this case. I mean, without knowing the total population of clients and its stratification (not given in the question), there's nothing to be said about this risk, and in fact any sample size (apart from N=1) could contain members from different populations (if there are any). In that case, N=100 could be seen as too big as well. So, why would a sample of N=100 be free from that risk (or at least have a low risk) and a sample size of N>100 would all of the sudden be so much more risky in this respect that it might be advisable not to increase the sample size? And, to turn this question around, what would be the standard operating procedure to protect oneself from this multi-population risk?


title: completely randomized VS matching - housing market correspondence test

text: I work in experimental economics, labor and housing market. I am trying to answer the following question: why do researchers use matched pairs of applications in labor market correspondence tests (typically in hiring context), whereas they use a completely randomized approach in the housing market correspondence tests (typically rental housing market)? I thought a possible explanation might be the following:   1. The completely randomized approach reduces the potential bias (i.e. similar probabilities for discriminators and non-discriminators to end up in the groups I compare).    2. The completely randomized approach also allows utilizing a larger sample, avoiding over/under representation of discriminators or non-discriminators.    3. The randomization process also shall be preferred to the matching approach (usually used in correspondence tests in the labour market), since the researcher might not be fully aware of which variables are relevant for the matching.    4. Finally, with the matching approach the detection probability will be limited. What do you think? Is the answer complete? Anything to be corrected/added?


title: Software (or webapps) for teaching kids statistics or probability?

text: I would like (in the distant future) to teach statistics to kids. For that matter, I'd be happy to know of software (obviously I am tending towards FOSS), or webapps, that are helpful in explaining statistical/probabilistic ideas to kids (or adults for that matter). This can be used either by the instructor, the kids, or both. Suggested format of the answer: Software name, what it help teach, who should use it, link.


title: How to interpret regression estimates

text: I am studying the effect a certain chemical and age on an outcome. As variables, `age` is a factor of with levels 1, 2, and 3 and `covar1` is continuous. After fitting the model               mod1 <- glm(out1 ~ covar1*age, family=poisson, data=df)      The estimates are               Parametric coefficients:                     Estimate Std. Error t value      Pr(>|t|)         (Intercept)    -0.6236832  0.1037201  -6.013 0.00000000185 ***     covar1          0.0204307  0.0101802   2.007        0.0448 *       age2            1.2523312  0.0649721  19.275       < 2e-16 ***     age3            0.8373618  0.0673236  12.438       < 2e-16 ***     covar1:age2    -0.0145580  0.0098667  -1.475        0.1401         covar1:age3     0.0109891  0.0101829   1.079        0.2805         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      How do I interpret the effects of `covar1` and `age` on my outcome?


title: Fitting a random walk

text: Is there such thing as a procedure to 'fit' an random walk? I have done analysis on the increments of my data and have fitted several distributions and performed goodness-of-fit tests. With this done, is that the random walk model 'fitted' and then I should test the validity of the model? Something that looks interesting is the fit of the t-distribution with changing df depending on the size of sample (daily change over a month (30obvs), daily change over a year (360 obvs) etc.) the t-distribution fits well for all but the df changes. Is there a way to extrapolate this into one model for future predicitions?


title: Making two equal groups from exponentially distributed data

text: I'm running an experiment where we are taking 1500 of a certain type of product and changing the price of half of them to see if this influences sales. There are 5 categories on which I will be comparing the data that are exponentially distributed (i.e. sales volume over time for different regions). I need to put the selected products into Groups A and Group B where Group A will be the products with the price change. I also need to decide how long to run this test to make it valid. A few questions that arise...   1. Before starting with any statistics... Is it proper to just take all the products with an odd index number and place those in Group A and place all the even indexed products in Group B? My issue with doing this is that the two groups will not be similar with respect to average sales volume for the different regions, but does that really matter? It seems to me that even if the values were different I'm just testing to see at what rate Group A or Group B changed; the similarity for the two starting groups is meaningless.   2. However, if it is important that the two groups are similar before starting the test, is it possible to "finagle" the data into equal groups with me selecting the products individually to go into Group A or Group B? Someone mentioned that this is improper since one you select items to get awesome p-values for a t-test you "lose independence". I don't really know what that means though.   3. I've been looking on the boards and saw some stuff on the "equivalence" package in R, however, even after using that a bit it mentions in the documentation that your data must be normal. How would I tweak this for exponential data?   4. Let's say I figure out a method to divide the two groups. How long should I have the experiment run? Any help would be much appreciated.


title: GLM with quasipoisson nested distribution

text: I want to make a GLM with a quasipoisson response distribution. My response variable is inflorescences by plant and my explanatory variable is coverage. For inflorescence counts, I used the average from four subsamples counted from each plant. However, somebody recommended that I use the original data. In this way I have four data for inflorescence for each plant (like nested designs). How can I write the script in order to consider the subsamples?               glm1 <- Glm(Inflorescences(4 data by plant)~coverage, quasiposisson)


title: How to interpret interaction between 2 independent variables with ANCOVA?

text: I'm measuring the influence of the use of male gaze and nudity in magazines to see if girls self-objectify or not. Now, I want to determine the influence of age.               Independent variables:   male gaze (yes / no)                               clothes (few clothes / lots of clothes)     Dependent variable:      Self-objectification score     Covariate:               Age      In the experiment, the participants were assigned to 4 conditions (`male gaze - few clothes`; `male gaze - lots of clothes`; `no male gaze - lots of clothes`; `no male gaze - few clothes`). In my **ANCOVA** , I only want to determine the main effects and interaction effect between `male gaze` and `clothing`. This was my output:   * no significant effect for `male gaze`, neither for `clothing` or `age`   * BUT significant effect for `male gaze * clothing` (interaction effect) So, apparently, `age` had no significant influence on self-objectification, but the interaction effect between `male gaze` and `clothing` was significant (.010). What does this mean? That the participants who were assigned to this condition ALWAYS objectify? I have difficulties interpreting this, especially since the **ANOVA** showed that there was no significant interaction effect between `male gaze` and `clothing`...


title: How are optional stopping rules based on e.g. sample confidence (width of confidence interval) biased?

text: Inspired by this: http://pss.sagepub.com/content/22/11/1359 In the context of open-ended data collection where the necessary sample size cannot be properly estimated, for the purpose of a frequentists test; I understand that a stopping condition based on the main outcome is circular. For example, if I stop sampling once my p value happens to be below .05, my p value is biased (so much as to be mostly worthless). However, say I choose another stopping rule, such as the width of my 95% confidence interval (without regard for other aspects of the test, such as e.g. if the CI includes 0), am I introducing any bias (but for, of course, CI width and related statistics)? As far as I understand it, this is not a problem in a Bayesian analysis, but I am wondering about the options for conditional stopping precluding frequentist tests.


title: Cumulative Incidence Function - Competing Risks Discrete Time Survival Analysis

text: As I'm doing a competing risks survival analysis, I need to construct a cumulative incidence graph rather than a Kaplan-Meier survivor curve. However, in Stata, the command to construct this graph is only available for continuous data so has anyone had experience in constructing a CIF using discrete time data and how did you do it?


title: Alternatives to the multinomial logit model

text: I am trying to estimate a model of occupational choice with three choices. Are there any alternatives to using the multinomial logistic regression when handling such unordered categorical outcomes? When dealing with binary dependent variables there seems to be several choices such as the LPM model as well as the binary probit and logit model. When dealing with unordered categorical variables the literature however keeps recommending the multinomial logit model without comparing it to alternatives.


title: Using the self-organizing map for sequences of categorical data

text: I have a number of vectors of categorical data (ex. {'re','ty','cf', ...} ) and I want to perform an unsupervised learning on them. I came across `self- organizing map` technique and I am trying to see if I can apply this method on my data. The vectors that I have are not similar in size and I was thinking if the `self-organizing map` can be applied on my data. Is there any version of this method that can do so?


title: What do NORB and CIFAR stand for?

text: The MNIST dataset is a standard benchmark data set of digit images. MNIST stands for 'Mixed National Institute of Standards and Technology'. The NORB dataset is a commonly used dataset of binocular images used in object recognition. What does the abbreviation NORB stand for? The CIFAR dataset is a commonly used dataset of RGB images used in object recognition. What does the abbreviation CIFAR stand for?


title: Logistic regression when some variables are measured just for a subset of samples

text: I have been asked an interesting modelling question by a clinician. The outcome he is trying to predict whether a particular type of cancer will occur. One of the covariates is an indicator variable if the patient previously had another kind of cancer. Other covariates are cancer diameter and margin size of this prior cancer. The margin size is how many millimetres the surgeon cut around the cancer when it was removed. For patients that have not had the prior cancer, these values are obviously missing. I believe that making diameter be 0 millimetres makes sense, but there can be no sensible value for margin size. What model formula should I use? What theory should I learn?


title: What is statistical Bias?

text: I am asking this question in context to section 4.1 in this paper: security control methods for statistical Database (http://www.utdallas.edu/~muratk/courses/privacy08f_files/stat_database_sec.pdf) As I am not statistician or Mathematician so I had hard time understanding it. If someone could explain what is statistical bias or the bias problem in section 4.1 of this paper, that would be great :). Thanks! :)


title: Can I create a composite scale with Likert and non-Likert scale data?

text: I would like to create a composite scale for different attitudes by calculating the mean score from different items on a questionnaire. Most of the items use Likert scales however some of the items on the questionnaire just use yes/no answers. Is it still possible to combine these scores to create a composite scale? How would I do this? e.g. Question 1 Likert scale 1-5 Question 9 Yes/No


title: Market mix modeling question

text: Suppose I want to model the sales of Coke sold at WalMart for a particular state in the U.S. Suppose I only have access to a year's worth of weekly sales data for that particular store. I want to model the sales of Coke. My question pertains to the modeling of catastrophic events or rare events. For example, suppose I want to model the effect of Hurricane Sandy on the sales of Coke for the particular WalMart. Also, I would like to model the effect of President Obama's re-election on the sales of Walmart. Intuitively, I feel I cannot model these directly. For modeling the influence of Hurricane Sandy, would it make sense to partition the data into two subsets where one subset includes the influence of Hurricane Sandy and the subset does not? Would an approval rating regarding the president's handling of Hurricane Sandy be a good variable that encapsulates the influence of Hurricane Sandy on the sales of Coke at that particular Walmart? As a side note, are there any good resources/books about market mix modeling?


title: SAS and T-Tests - Degrees of Freedom

text: Why is it that in a normal t-test we have n-2 degrees of freedom while in the SAS output, and consequently the confidence interval constructed from this output, we end up using our 'error' degrees of freedom. As an example say we have 1 DF for our linear model and 14 DF 'error'. Then there are 15 DF in total. However, we take t with 14 degrees of freedom = 15 - 1 rather than 15 - 2. This is n-1 to my knowledge, not n-2. Moreover, in a regular t-test we are taking n-2 since we have two degrees removed for b0 and b1 as parameters. Why is the linear model only given 1 degree of freedom in SAS?


title: Significant difference between difference scores of 2 groups on two measures?

text: I am investigating a discrepancy between male and female self reports of sexual experiences. The original survey consists of a female version (asking about victimization) and a male version (asking about perpetration). Typically, when given the original version, female's reported rates of victimization are about 2/3s higher than male rates of perpetration. I have modified the original survey (both male and female versions) in order to determine whether the wording of the modified version will have an impact on the female/victim--male/perpetrator discrepancy. One of my hypotheses is that the modified version will produce a narrower discrepancy between female reports of victimization, and male reports of perpetration. I need to figure out what test (or series of tests) I can use to determine if there is a significant difference between the discrepancy rate of original survey, and the discrepancy rate of my modified version. Additional info:   * males and females are not matched, and I have different sample sizes of males and females   * each subject was administered both versions (original and modified) of the survey, according to gender   * subjects answered the original survey first, and then were given the modified survey   * my data will be nominal -- e.g. "Yes" I've had this experience, or "no" I haven't had this experience.


title: Scoring predictions of an ordinal variable

text: I read about using scoring rules to evaluate the performance of predictive models. In the Wikipedia article about the Brier score, it is stated: > The Brier score is appropriate for binary and categorical outcomes that can > be structured as true or false, but is inappropriate for ordinal variables > which can take on three or more values (this is because the Brier score > assumes that all possible outcomes are equivalently "distant" from one > another). What are some examples of scoring rules used for evaluating predictions of ordinal variables? For example, if the set of possible outcomes is $\\{0, 1, 2, 3, 4, 5, 6, 7, 8, 9\\}$ and a predictive model yields $p_6 = 0.3$, $p_7 = 0.5$, $p_8 = 0.2$ (where $p_i$ is the believed probability that the outcome will be $i$), this predictive model should receive a lower score if the outcome is $3$ than if it is $4$.


title: Fallacy in p-value definition

text: I have two more questions. Why Statsoft's textbook says that > the p-value represents the probability of error that is involved in > accepting our observed result as valid whereas Wikipedia gives a slightly different definition: > p-value is the probability of obtaining a test statistic at least as extreme > as the one that was actually observed, assuming that the null hypothesis is > true Wikipedia warns us against confusing Fisher's p-value criteria with Pearson's Type I error (probability of taking a wrong decision). I read that statsoft's definition falls into the Type I fallacy. Am I wrong? Can you explain how these definitions related?


title: Meaning of R-Square in this Context

text: Good evening, I am attempting the following question where I was able to successfully find the R-square (the coefficient of determination) at 0.8945. While I am able to identify this, my weak foundation in regression leave me puzzled to what this actually means. Appreciate some advice please. > An accountant for a large department store would like to develop a model to > predict the amount of time it takes to process invoices. Data are collected > from the past 30 working days. The number of invoices processed and the > completion time (in hours) are recorded. All the computer outputs are given > in the end of this exam paper. You may directly use the computer outputs to > answer the following questions if necessary. ![enter image description here](http://i.stack.imgur.com/Nh67g.jpg) ![enter image description here](http://i.stack.imgur.com/pIyJV.jpg)


title: boosting with ada package. How do I take the most probable answer from predict?

text: I trained classifier using ada. Now I executed:               predict(adaDol,newdata=cords())      and received response:                [1] FALSE TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE      [25] TRUE  TRUE  FALSE TRUE  TRUE  TRUE  TRUE  TRUE  FALSE TRUE  TRUE  TRUE  TRUE  FALSE FALSE FALSE FALSE TRUE  FALSE TRUE  FALSE TRUE  FALSE TRUE      [49] FALSE FALSE     Levels: FALSE TRUE      and some probabilites               predict(adaDol,newdata=cords(), type="probs")          $probs                [,1]      [,2]      [1,] 0.5037949 0.4962051      [2,] 0.4857610 0.5142390      [3,] 0.4654769 0.5345231      [4,] 0.4711432 0.5288568      [5,] 0.4714576 0.5285424      [6,] 0.4902077 0.5097923      [7,] 0.4813705 0.5186295      [8,] 0.4772360 0.5227640      [9,] 0.4595959 0.5404041     [10,] 0.4635512 0.5364488     [11,] 0.4828989 0.5171011     [12,] 0.4602214 0.5397786     [13,] 0.4857656 0.5142344     [14,] 0.4886816 0.5113184     [15,] 0.4677266 0.5322734     [16,] 0.4805374 0.5194626     [17,] 0.4889051 0.5110949     [18,] 0.4836400 0.5163600     [19,] 0.4864160 0.5135840     [20,] 0.4925470 0.5074530     [21,] 0.4854747 0.5145253     [22,] 0.4876254 0.5123746     ...      How do I calculate the most probable outcome the way ada does it while training?


title: Does anyone know "nice" ways of plotting confidence intervals for use in presentations?

text: I am now plotting some bars and pictures for use in presentation. Really interested in different ways of showing confidence intervals on the plots.


title: Comparing the Goodness of fit of two polytomous set of items not nested

text: I'm working on two identical set of items that only differ in their frequency- of-occurence dimension (in the first is max.7 while in the second is max.4). They are both fit with the Generalised partial credit model. I would like to understand if the scale has a better fit with 7 categories or with 4. Is there a test I could use?


title: finding position by regression fit

text: I try to estimate an underfloor pipe position based on a temperature distribution on a surface of the floor. I know there is the pipe emitting a hot temperature all way along. I can measure a distribution of this temperature in a rectangular room. I know that the pipe can be in two shapes only as a straight line or a bended section – second order polynomial. When I plot the data it shows that the pipe comes into the room as the straight section, after while it bends and it becomes the straight section again The coordinates for each sampled temperature is as a Cartesian coordinates (x,y) in the room. And in addition I have the temperature in the sample point, hence I have points: $s_1(x_1,y_1,t_1), s_2(x_2,y_2,t_2),… s_n(x_n,y_n,t_n)$ I have many sample points in this grid and can measure more if necessary. I am thinking to apply a regression fit to this data, but can’t find appropriate information how to apply it to such problem. I want to finish with three equations for three section each, as 1st and 2nd order polynomials. I have two problems: I am not sure how to formulate it as a regression to prove I am correct (I want to do it in R later on) I don’t know how to separate three sections by an algorithm instead of imprecisely read it from the plot. Can you point me to solution or some additional information about problem like that?


title: A priori power analysis for generalized linear mixed-effects model

text: Let's say I have a data set (with a binary DV) which I used a generalized linear mixed-effects model on. I am currently trying to run a simulation-based power analysis to determine the sample size I should try to collect for a replication attempt. I know about the simulate function for fitted lme4 model objects, which will create simulated responses that follow the size and structure of the actual dataset. My question is, **how can I use this fitted model object to run an a-priori power analysis for my replication?**. In other words, is there an easy way to extract the necessary information from the simulate function, vary the sample size, then to simulate responses? EDIT: Here is a simple example with ten subjects that should give me a good starting place. The data here are totally fake so the effects are not close to what I have in the pilot data. The two variables of interest are age in months (between subjects) and whether the stimulus is "similar" or "dissimilar" (within subjects). Each individual is exposed to three items within each similarity group, for a total of six binary responses for each subject. As seen in the model, I am interested in a possible interaction between age and similarity and include a random slope of similarity to do this. Let me know if I need to modify anything/provide more information to make things easier.               set.seed(28)          fakeData <- data.frame(ID = rep(1:10, each = 6),                            months = rep(round(runif(10, min = 60, max = 80)), each = 18),                            similarity = rep(c("Sim", "Dissim"), each = 3),                            response = sample(c(0, 1), 60, replace = TRUE))          fakeModel <- glmer(response ~ months * similarity + (similarity | ID), family = binomial, data = fakeData)     summary(fakeModel)


title: help with expectation propagation: updating factors

text: I have been trying unsuccessfully for a while to understand and setup an inference problem and use EP with it. I have tried to extract a simple case to highlight he issue. I hope someone will be kind enough to offer some insight into this. I have a graphical model as shown in the image below: ![A simple PGM](http://i.stack.imgur.com/vyHQU.png) I am trying to estimate a parameter set $w$. Note that $y$, $x$ and $\phi$ are observed. Now, this can be written as a factor graph as follows: ![Corresponding factor graph](http://i.stack.imgur.com/8WlRD.png) Now, the parameter $w$ is modelled as a multivariate distribution with a mean vector specified by $m$ and has a covariance structure $\Sigma$. The idea is that these parameters need to be spatially regularised and this is enforced by this covariance structure. Now, the likelihood can be broken down as product of individual likelihoods (I know…the factor graph should show multiple factors one for each likelihood). Anyway, the likelihood can be written as: $$ P(y|w,x,\phi) = \prod_n^{N} p(y_n|w_n,x_n,\phi) $$ and an individual likelihood can be written as: $$ p(y_n|w_n,x_n,\phi) = \frac{\phi}{2\pi} \exp^{-0.5 (y_n - t(x_n,w_n))\phi(y_n - t(x_n,w_n))} $$ So, this is minimising the sum of square difference between the observed $y_{i}$ and the transformed $x_{i}$ according to the estimated parameter $w$ and some transformation function $t$. Now, I have been watching Tom Minka's lecture and doing reading on EP. However, according to the lectures I will be refining the factors in turn. So, imagine that I am refining the ith likelihood term, so I will have a product as follows: From all the likelihood terms, I will have a term $q^{\backslash i}(y|w^{\backslash i}, x^{\backslash i}, \phi)$. So this is a product of all the approximated likelihood terms (except the ith term). Now, this has to be multiplied with the approximate prior term $q(w)$. Now $q(w)$ is a multivariate distribution with the same form as $P(w)$, a MV normal distribution. So, I am not sure how to do this step of $q(w)q^{\backslash i}(y|w^{\backslash i}, x^{\backslash i}, \phi)$. I have a feeling that I am missing a trick where the factor graph perhaps need to be refined further but am totally at a loss on how to do it. I would really appreciate any tips/guidance you can give me. I have been really banging my head against this for weeks.


title: How to reconstruct the original data when using manifold learning?

text: I'm using Isomap to reduce the dimensionality of my data. Isomap use geodesic distance rather than Euclidean distance to perform a MDS. Now I want to reconstruct my original data with a lower dimensional coordinate. How can I do that? When the lower dimensional coordinate lies on a data point, things could be easy. Just perform a certain kind of trace back. But what if it lies out of data points? ** I intended to put the tag "manifold-learning", but I cannot find it. Really?


title: What is the formula for lognormal hazard?

text: I'm plotting a bunch of survivor and hazard curves. The lognormal survivor function is: $S(t)=1-\Phi(\frac{log(t)-\mu}{\sigma}) $ Where $\mu$ is the scalar parameter. From a website (http://www.itl.nist.gov/div898/handbook/eda/section3/eda3669.htm), the null model gives a hazard function of: $h(t,σ)=\frac{\frac{1}{tσ}ϕ(\frac{\ln t}{\sigma})}{Φ(\frac{−\ln t}{\sigma})}$ where ϕ is the probability density function of the normal distribution and Φ is the cumulative distribution function of the normal distribution. Where does the scale parameter come in?


title: Observed information matrix is a consistent estimator of the expected information matrix?

text: I am trying to prove that the observed information matrix evaluated at the weakly consistent maximum likelihood estimator (MLE), is a weakly consistent estimator of the expected information matrix. This is a widely quoted result but nobody gives a reference or a proof (I have exhausted I think the first 20 pages of google results and my stats textbooks)! Using a weakly consistent sequence of MLEs I can use the weak law of large numbers (WLLN) and the continuous mapping theorem to get the result I want. However I believe the continuous mapping theorem cannot be used. Instead I think the uniform law of large numbers (ULLN) needs to be used. Does anybody know of a reference that has a proof of this? I have an attempt at the ULLN but omit it for now for brevity. I apologise for the length of this question but notation has to be introduced. The notation is as folows (my proof is at the end). Assume we have an iid sample of random variables $\\{Y_1,\ldots,Y_N\\}$ with densities $f(\tilde{Y}|\theta)$, where $\theta\in\Theta\subseteq\mathbb{R}^{k}$ (here $\tilde{Y}$ is a just a general random variable with the same density as any one of the members of the sample). The vector $Y=(Y_1,\ldots,Y_N)^{T}$ is the vector of all the sample vectors where $Y_{i}\in\mathbb{R}^{n}$ for all $i=1,\ldots,N$. The true parameter value of the densities is $\theta_{0}$, and $\hat{\theta}_{N}(Y)$ is the weakly consistent maximum likelihood estimator (MLE) of $\theta_{0}$. Subject to regularity conditions the Fisher Information matrix can be written as $$I(\theta)=-E_\theta \left[H_{\theta}(\log f(\tilde{Y}|\theta)\right]$$ where ${H}_{\theta}$ is the Hessian matrix. The sample equivalent is $$I_N(\theta)=\sum_{i=1}^N I_{y_i}(\theta),$$ where $I_{y_i}=-E_\theta \left[H_{\theta}(\log f(Y_{i}|\theta)\right]$. The observed information matrix is; $J(\theta) = -H_\theta(\log f(y|\theta)$, (some people demand the matrix is evaluated at $\hat{\theta}$ but some don't). The sample observed information matrix is; $J_N(\theta)=\sum_{i=1}^N J_{y_i}(\theta)$ where $J_{y_i}(\theta)=-H_\theta(\log f(y_{i}|\theta)$. I can prove convergence in probability of the estimator $N^{-1}J_N(\theta)$ to $I(\theta)$, but not of $N^{-1}J_{N}(\hat{\theta}_N(Y))$ to $I(\theta_{0})$. Here is my proof so far; Now $(J_{N}(\theta))_{rs}=-\sum_{i=1}^N (H_\theta(\log f(Y_i|\theta))_{rs}$ is element $(r,s)$ of $J_N(\theta)$, for any $r,s=1,\ldots,k$. If the sample is iid, then by the weak law of large numbers (WLLN), the average of these summands converges in probability to $-E_{\theta}[(H_\theta(\log f(Y_{1}|\theta))_{rs}]=(I_{Y_1}(\theta))_{rs}=(I(\theta))_{rs}$. Thus $N^{-1}(J_N(\theta))_{rs}\overset{P}{\rightarrow}(I(\theta))_{rs}$ for all $r,s=1,\ldots,k$, and so $N^{-1}J_N(\theta)\overset{P}{\rightarrow}I(\theta)$. Unfortunately we cannot simply conclude $N^{-1}J_{N}(\hat{\theta}_N(Y))\overset{P}{\rightarrow}I(\theta_0)$ by using the continuous mapping theorem since $N^{-1}J_{N}(\cdot)$ is not the same function as $I(\cdot)$. Any help on this would be greatly appreciated.


title: How to determine reliability of a change in an arbitary metric

text: I have been using a/b testing in order to optimise changes I make to my websites. Previously, I have set a single conversion goal and have ended up with a conversion rate for each group. I then used R (with the `power.prop.test` function) to calculate the power and significance of the change in conversion rate to decide when it is safe to stop the test and be confident that the values have meaning. I am now able to get information on a wide range of metrics about the site for each group. However I am not sure how to do the same test for an arbitrary metric - e.g., revenue. Could someone point me in the right direction? I am happy to use R so just pointing me at the right set of functions would be great, but I do like to try and understand what I am calculating so any explanations or links to useful site would be even better.


title: Forecasting & social influence data/experiment - Seeking research strategies

text: In my experiment, individuals assign probabilities to the likelihood of future events, and update their forecasts as frequently as they like. Most questions stay open (receiving new forecasts) for months. Forecasters are placed in different experimental conditions. Some forecast in isolation, some can observe forecasts made by others, and a final group can not only see other's forecasts, but is assigned to one of twenty-five teams of 15-20 people each, and encouraged to share information with their team. One of my goals is to characterize social influence. e.g. "If an average person raises their forecast 10%, what effect does it have on their teammates subsequent forecasts?" But their are many other things to consider: if the forecast change is accompanied by a verbal justification, how does that change its influence? If the forecast change moves the forecaster closer to the team mean, does that (as I would expect) imply it will have less influence on subsequent teammates' forecasts than if it moves the forecaster further from the team mean. I think it would be interesting to compare forecaster influence with forecaster accuracy. I assume they are positively correlated, but how strongly. Can we distinguish good/bad influences in terms of their effect on accuracy? Accuracy can be measured easily enough using something like the brier scoring rule. But estimates of social influence will inevitably rest on some pretty heroic assumptions. How can I distinguish between a team reacting to a previous forecast, and a team reacting to some piece of news which also affected the previous forecast? My questions for the crowd are as follows: 1) Given that I have this interest in influence, how would should I go about estimating it? Might I take advantage of random assignment to teams and/or to forecasting conditions without teams at all? 2) What other sorts of questions might I try to answer with this data?


title: Explaining augmented Dickey-Fuller regression output

text: I have monthly returns data going back to 1991 and I'm trying to work out if the data has a tendency to mean revert over time. In order to work this out I've used the Augmented Dickey Fuller test on GRETL - an example of my output is here: ![output](http://i.stack.imgur.com/ifOUO.jpg) I've also run checks on regression models without a trend, and without constant (this form had the lowest p value). I'm now struggling to understand what this output means. I know that the stars refer to how significant the coefficients are, but what does everything else mean. I'm guessing the constant refers to the mean it trends back to? How does the time term fit into this? What impact do the lags have? How could I plot this regression output on the graph of monthly returns? Sorry, I know this is all pretty basic stuff but I'm really struggling.


title: Differencing an i.i.d. time series

text: Why differencing a normal, i.i.d. time series $X$, generates a negatively correlated series $Y$?               > Acf(rnorm(5000, 0, 40))     > Acf(diff(rnorm(5000, 0, 40)))      ![enter image description here](http://i.stack.imgur.com/c1FXP.png) I stumbled on this problem while looking at a time-series that I wanted to predict using an ARIMA model. The original series appeared to be non- stationary by looking at the plot, so I decided to apply a difference at lag 1 and check the acf/pacf. This seemed to indicate an MA(1) model but obviously there is something I'm missing. ![enter image description here](http://i.stack.imgur.com/ERY6A.png)


title: Problem understanding the logistic regression link function

text: I am trying to learn the logistic regression model. I came to know that there is no linear relationship between predictor variables and response variables since response variables are binary (dichotomous). The link function used for logistic regression is logit which is given by $$ \log \frac {p}{1 - p} = \beta X $$ This tells that the log odds is a linear function of input features. Can anyone give me the mathematical interpretation of how the above relation becomes linear i.e. how logistic regression assumes that the log odds are linear function of input features? Since I am poor at statistics, I can't understand complex mathematical answer.


title: How to present a empirical study when using econometric models?

text: I've got a (probably easy) question in how to handle empirical studies, when there are a lot of effects involved. I have a whole bunch of variables and I'd like to analyze just a few of them. But the problem is, that the model is wrong... So standard errors and coeficients itself are probably biased and so t-statistics are as well. So general: Everything is wrong. Pretty frustrating task to find out how to handle this problem, when it is not possible to say what coeficeints have a clear influence on $y$. What would you do in this case? It's possible to compact the coeficients, but the problem is still present. Do you have some experiences how to handle this problem? Or is there anyone who knows a good paper where it's been discussed? Fyi: I'm going to do Cross Validation afterwards, to compare models... But it's still required to make an analysis of the estimation before looking which model is good. And I'm bounded to do OLS, before looking for better models. The general question is: How are other studies dealing with biased std. errors or coefficients? Please help :( Edit: I'm analyzing a whole bunch of effects on wage. Thererefor I've a lot of effects. I know that heteroskedasticity occurs, wage is skewed and the sample size is relatively small. I'm not interested in changing the model, since I've to do OLS without transforming variables or s.th. Just the regular OLS. Unfortunately I don't really know how to interpret all the effects, when I can't get rid of the non significant ones, because significance isn't good defined (because of bias). Is there a theory that says in general:"Although bias occurs you can assume that the effects with high significances are more clear different from 0, than other effects that are not significant?"


title: If I count cells should I consider them as Poisson distributed?

text: I calculate (with flow cytometry) pecentage of lymphocytes with a specific receptor (Lph*) as a ratio to general number of lymphocytes (Lph). Should I consider them (Lph*) as Poisson distributed?   (My data set is here.)


title: The number 2 in Deviance

text: What is the additional coefficient 2 in AIC, or other deviance criteria? In the books or on the internet I keep reading "historical reasons". What the heck are those historical reasons? Is there a merit for the 2 keep appearing? http://en.wikipedia.org/wiki/Deviance_(statistics) http://en.wikipedia.org/wiki/Bayesian_information_criterion Clearly the correct deviance for $AIC$ is $LogLikelihood - k$. Not only that it is negated it is also scaled by 2 in presentations: $AIC=-2LogLikelihood + 2k$ But then to fix it in relative likelihood we divide by 2. $Relative Likelihood = e^{(AIC_{min}−AIC_i)/2}$ Thanks


title: Categorical variable as response in poisson regression

text: I have data on damages on flowers from different treatments. The damage was originally count data (number of damages per flower) but the person collecting the data categorized the data in four levels - no damage (0) to a high damage level (3). I know I could use logistic regression to model the odds but I wonder if it is possible to use poisson regression instead and consider the observed level of damage as a count variable? If possible, would the intepretation of the dichotomized treatment coefficients be the difference in the logs of expected counts of damage levels for treatment compared to baseline while holding the other variables constant in the model? The weird intepretation makes me think this is not at all possible.


title: Why there are two different results?

text: ![enter image description here](http://i.stack.imgur.com/okkhr.png)I found my four independent variables have a high positive association with my depepednet variable. However, except one independent variable, P value of all other indepdent variables are not significant in my regressino analysis. My hypothesis is that each of my independent variables separately has a postive relationship with my dependent variable, so I have four hypothese. I was wondering if I should reject three of my hypothese under such situation? Sorry, forgot to say, There is no collinearity problem with my independent varibales Thank you. More information:Please see attached picture.![enter image description here](http://i.stack.imgur.com/Z0woI.png) It is my research model. After bivariate correlation analysis and partial correlation, all independent variables are positively associated with dependent variable. A hierarchical regression was used with control variable.I first put control variable in first block of SPSS, and then add the 4 other independent variables in second block. The results are in attached docuement. Three independent variables are not significant.... Correlation table: ![enter image description here](http://i.stack.imgur.com/Gfnol.png)


title: Fitting a generalized least squares model with correlated data; use ML or REML?

text: Reading the Linear Mixed Model (LMM) literature I am aware that fitting a model using REML provides better estimates of variance parameters than fitting via ML. However, we should not compare nested models fitted with REML that have different fixed effects. Recently, I have been fitting some models using GLS via the `gls()` function in the **nlme** package for R. The default fitting method for that function is REML. Do the same principals of REML vs ML for LMM also apply to GLS? Specifically, I am fitting a model with and without a linear trend with a correlation structure in the residuals:               m1 <- gls(Response ~ Time, data = foo, correlation = corAR1(form ~ Time))     m0 <- gls(Response ~ 1, data = foo, correlation = corAR1(form ~ Time))      In the above, I should fit the models using ML as they have different fixed effects. Is this correct? Secondly, consider two GLS models that differ in the correlation structure:               m1 <- gls(Response ~ Time, data = foo, correlation = corARMA(form ~ Time, p = 1))     m2 <- gls(Response ~ Time, data = foo, correlation = corARMA(form ~ Time, p = 2))      What fitting method should ideally be used here? REML or ML? Here my intuition would say fit via REML as we are estimating (co)variance parameters. Is my intuition correct or have I got this all mixed up?


title: Mean comparison for statistical significance

text: I was told to come here by a friend to ask a stats question. I took Social Stats and Methods of Social Research about 4 years ago and to be honest, I don't remember much about it. Now I'm writing a research paper and I can't find any basic tutorials to assist with a problem that would probably be pretty simple to a Stats major. So I'm running tests on the lunar effect, and I guess the qualitative details aren't important. So In "Group F" I have a sample size of 44 Days, and when measuring incidence within the days, the total was Σ=1593, which generates a mean of 36.20 per day on average. In "Group J" I have a sample size of 27 Days, and when measuring incidence within the days, the total was Σ=997, which generates a mean of 36.93 per day on average. I'm trying to see if there is a statistically significant difference between the means comparing Group F and Group J. I guess CL=95% or α=0.05 (I'm not sure) and I don't remember anything about standard deviation. So my question is, what is my step by step process to figure this out (becuase i have to do it 9 more times) and how do I find out my standard deviation and CL/α? Edit) Greg Snow, I'll answer the questions first. Yes, I have each individual number for each day. So tell me if this is how you find the standard deviation. 1) Subtract the mean from each numerical result from each day. 2) Square all of those values. 3) Add the results 4) Take the square root of that total... That's correct, right? There where certainly no 0's. lowest number I ever saw was something in the high 20's (I don't have the numbers in front of me). Highest was around low 70's (those maxes an mins are considering both groups)... No. Each day is independent and isn't dependent on days before or after. The days don't overlap. The number of days is more related to their frequency of occurrence (full moon nights vs. non-full moon nights). I guess you could call full moon days predetermined lol. But the non-full moon days are the 5th, 10th, 15th, 20th, and 25th day after the full moon (using a 3 full moon day method). That should answer your questions EDIT) Ok I see. thank you for your help! However, I have a question. What's the difference between dividing the variance by N or N-1? I know it's Population vs. Sample, but how do you know which to use?


title: Chi-squared versus logistic regression

text: I used `summary.formula` from `Hmisc` with continuous `Age` and binary outcome `O` with `test=TRUE`. This returned a p-value for `Age` predicting `O` (if I understand this correctly). I then ran a `glm` using `Age` and `O` (univariate logistic regression), which returned a different p-value. I thought that the p-value should be the same? **Edit:** `library(Hmisc)` `summary.formula(Outcome ~ cut2(Age, seq(15,75,10)), method="reverse", test=TRUE)` p-value = 0.8 `x = glm(Outcome ~ Age, family=binomial(link="logit"))` p-value = 0.4 (I'm getting quite confused: the example above which gives the 0.8 p-value is binning the ages, but I was told that the p-value it returns is for all the ages. However, if I run the summary without the binning, the p-value is 0.5.)


title: Comparison of 2 distributions

text: I have 2 distributions, 1 is, as far as im aware, normally distributed. Distribution 1 is the control group. 1st distribution Mean = 0.000002757; Median = 0; StDev = 0.00119307; Number of data points = 91601 2nd distribution Mean = 0.000058; Median = 0.000125; StDev = 0.001646243; Number of data points = 94045 **(Question 1)** Im assuming i can do a hypothesis test of the difference of the means to determine if the 2nd distribution is significantly different to the first?               1) H0 = No significant change.          2) Pick 5% significance level.          3) Z-table gives 1.65.          4) Calc Stdev of distribution of the difference of the means = 6.66*10^-6          5) Find threshold value =>(6.6610^-6)*1.65 = 1.098*10^-5      Since the difference of the difference of the means > 1.098*10^-5 i can reject the null hypothesis? **(Question 2)** Does it matter that the sample sizes are slightly different, regarding the calc for step 4? **(Question 3)** If my 1st distribution was not normal, could i still perform the same calculation to determine if the 2nd distribution was significantly different?


title: Is ARMA(0,0) equivalent to white noise?

text: If the EACF of my TS suggests ARMA(0,0) and the Box-Ljung test does not suggest my TS has correlation, can I conclude that my TS is white noise or merely that there is no reason to suspect that it is not TS? Conversely, if my TS is white noise, will I necessarily get ARMA(0,0) from EACF, ACF, PACF, etc?


title: Simple exponential smoothing

text: I simulated a time series using expressions (3.10a), (3.10b) from (Hyndman et al., 2008). Next, I'd like to use a simple exponential smoothing method to forecast for the next period. For a given initial point $\hat{y}_1$ I obtained the smoothing parameter $0<\alpha<2$, etc. My question is as follows: am I going to predict the next value of the local level $l_t$ or the next value of the observed time series $y_t$? Regretfully, I still do not understand how one can predict a quantity which includes a random component, and if I got it right, time series without trend consists of level and such a component.


title: Parameter estimation with generalized linear models

text: By default when we use a `glm` function in R, it uses the iteratively reweighted least squares (IWLS) method to find the maximum likelihood estimation the parameters. Now I have two questions.   1. Does IWLS estimations guarantee the global maximum of the likelihood function? Based on the last slide in this presentation, I think it does not! I just wanted to make sure of that.   2. Can we say that the reason for question 1 above is because of the fact that almost all the numerical optimization methods may stuck at a local maximum rather than a global maximum?


title: How to justify the error term in ANOVA?

text: A probably very basic question about multi-factorial ANOVA. Assume a two-way design where we test both main effects A, B, and the interaction A:B. When testing the main effect for A with type I SS, the effect SS is calculated as the difference $RSS(1) - RSS(A)$, where $RSS(1)$ is the residual error sum of squares for the model with just the intercept, and $RSS(A)$ the RSS for the model with factor A added. My question concerns the choice for the error term: **How do you justify that the error term for this test is typically calculated from the RSS of the full model A + B + A:B that includes both main effects and the interaction?** $$ F_{A} = \frac{(RSS_{1} - RSS_{A}) / (df_{RSS 1} - df_{RSS A})}{RSS_{A+B+A:B} / df_{RSS A+B+A:B}} $$ ... as opposed to taking the error term from the unrestricted model from the actual comparison (RSS from just the main effect A in the above case): $$ F_{A} = \frac{(RSS_{1} - RSS_{A}) / (df_{RSS 1} - df_{RSS A})}{RSS_{A} / df_{RSS A}} $$ This makes a difference, as the error term from the full model is probably often (not always) smaller than the error term from the unrestricted model in the comparison. It seems that the choice for the error term is somewhat arbitrary, creating room for desired p-value changes just by adding/removing factors that aren't really of interest, but change the error term anyway. In the following example, the F-value for A changes considerably depending on the choice for the full model, even though the actual comparison for the effect SS stays the same.               > DV  <- c(41,43,50, 51,43,53,54,46, 45,55,56,60,58,62,62,     +          56,47,45,46,49, 58,54,49,61,52,62, 59,55,68,63,     +          43,56,48,46,47, 59,46,58,54, 55,69,63,56,62,67)          > IV1 <- factor(rep(1:3, c(3+5+7, 5+6+4, 5+4+6)))     > IV2 <- factor(rep(rep(1:3, 3), c(3,5,7, 5,6,4, 5,4,6)))     > anova(lm(DV ~ IV1))                           # full model = unrestricted model (just A)               Df  Sum Sq Mean Sq F value Pr(>F)     IV1        2  101.11  50.556  0.9342 0.4009     Residuals 42 2272.80  54.114          > anova(lm(DV ~ IV1 + IV2))                     # full model = A+B               Df  Sum Sq Mean Sq F value   Pr(>F)         IV1        2  101.11   50.56  1.9833   0.1509         IV2        2 1253.19  626.59 24.5817 1.09e-07 ***     Residuals 40 1019.61   25.49                               > anova(lm(DV ~ IV1 + IV2 + IV1:IV2))           # full model = A+B+A:B               Df  Sum Sq Mean Sq F value    Pr(>F)         IV1        2  101.11   50.56  1.8102    0.1782         IV2        2 1253.19  626.59 22.4357 4.711e-07 ***     IV1:IV2    4   14.19    3.55  0.1270    0.9717         Residuals 36 1005.42   27.93      The same question applies to type II SS, and in general to a general linear hypothesis, i.e., to a model comparison between a restricted and an unrestricted model within a full model. (For type III SS, the unrestricted model is always the full model, so the question doesn't arise there)


title: If F = F1+F2 is some measured random field, how to find out if F2 has been correctly estimated and subtracted?

text: I have a random field, F, defined on the surface of a sphere pixellized into about a million points. I know that F must consist of two components - F1 and F2, F1 being a true signal and F2 being some contaminant. They cannot be measured independently but it is known that they should be completely uncorrelated to each other. F2 is estimated using some theoretical method and subtracted from F to leave behind F1. My question is: how can I check whether there is any residual contaminant left behind in F1 ? What would be a good statistical to do this check ? Any ideas would be welcome.


title: What's the appropriate spss test?

text: ### Context I have 25 men and 25 women as participants, and they did exactly the same thing: Each of them heard an attractive dialogue and they had to choose between a photo of a woman in red and a woman in green. Then each of them heard an unattractive dialogue and they chose again between the red and the green shirt. My hypothesis is that men are attracted to women in red in contrast to women. So, I 've got gender (2levels:0=male, 1=female), attraction (1=yes, 0=no) and trials/colour (0=green, 1=red) I am interested in showing that the colour (red) predicts attraction as far as men are concerned and that there are gender differences! Men are much more attracted to it. ### Questions   * Should I use a test which combine all of them or a test that combine first the attraction and colour for men, then the attraction and colour for women and then the gender differences? I've been told to run a pearson's correlation (one for men and one for women) but I think that it demands interval data.I was also told to use the chi- square test but it's not for participants who parrticipated in the same experimental conditions.   * What about repeated measures logistic regression? If so, can you give me some advice how to process my data? Any suggestions?  P.S. I really appreciate everyone who answered my previous posts here and here, they helped me to go one step further!


title: Forecasting irregular pattern and volume

text: I have tried a number of different models to forecast the time series shown below, but so far I haven't found any models that satisfy me. I am looking for ideas for a suitable model. ![enter image description here](http://i.stack.imgur.com/KSvjl.png) The objective of the forecasting is two-fold:   1. predicting the interval and thereby determining when the next peaks will happen.    2. What will the volume (size of peak) be on the next occurrences.  The time series data behind the chart are below:               sampleTS <- structure(c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2300L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 2700L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2750L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 2500L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 2600L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1700L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 2356L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      1967L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 3130L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L,      0L, 0L, 0L), .Dim = c(297L, 1L), .Dimnames = list(NULL, "V1"), .Tsp = c(15340,      15340.8109589041, 365), class = "ts")


title: How to test the statistical significance of AUC?

text: (Step 1) Using my predictive model, I predicted 1000 scores for my sample dataset. (Step 2) I then calculate the random score using the same method for a randomized dataset. I firstly fit the distribution of the random score. (Step 3) For each of my predictive score (1000 scores, in step 1), I calculated the p-values of getting a score larger than my predictive score for my sample sample dataset. Thus, 1000 p-values for my sample dataset are obtained. (Step 4) As the real classification is known, by looking at the enrichment of true positive, I found when filter the sample dataset by p-values < 0.05. The best true positive values is obtained, which represents about 150 data from my sample dataset. I then, want to test the predictive power of my model by doing AUC of the ROC (sensitivity vs 1- specificity plot). However, I am facing a problem now, should I include all 1000 data for ROC plot to get the AUC, or should I only include those 150 data (p < 0.05) for my AUC analysis? When I said p < 0.05, the p-value to obtain a score higher than my predictive score by random. In general speaking, does it means that 50% of my data are obtained by chance? * * * ## Edit Thanks for the comments from @AlefSin, @steffen and @Frank Harrell. For easier to discuss, I have prepare a sample dataset (x) as follows:   * My model predicted score (assume it is normal distributed with mean=1, sd=1)   * random set (assume also has mean=1, sd=1)   * Probability for each predicted score   * class prediction are listed as below, as listed in four columns * * *               x <- data.frame (predict_score=c(rnorm(50,m=1, sd=1)))     x$random <- rnorm(50, m=1, sd=1)     x$probability <- pnorm(x$predict_score, m=mean(x$random),sd=sd(x$random))     x$class <- c(1,1,1,1,2,1,2,1,2,2,1,1,1,1,2,1,2,1,2,2,1,1,1,1,2,1,2,1,2,2,1,1,2,1,2,1,2,2,1,1,1,1,1,1,2,2,2,2,1,1)      I then did AUC for all data as follows for all data points:               library(caTools)     colAUC(x$predict_score, x$class, plotROC=T, alg=c("Wilcoxon","ROC"))      result:                      [,1]     1 vs. 2  0.6      Let's said if the enrichment of true positive is higher (Your runs may be differnt from mine, as rnorm give different results everytime) when I filter the dataset by p < 0.5, I did the AUC for a subset of the data as folows:               b <- subset(x, x$probability < 0.5)     colAUC(b$predict_score, b$class, plotROC=T, alg=c("Wilcoxon","ROC"))      result:                            [,1]     1 vs. 2 0.7401961      My question is: When I do AUC analysis, is it a must to do the analysis with the whole dataset, or should we do filter the dataset first based on enrichment of true positive or what ever criteria before doing AUC?


title: How to test whether a regression coefficient is moderated by a grouping variable?

text: I have a regression done on two groups of the sample based on a moderating variable (say gender). I'm doing a simple test for the moderating effect by checking whether the significance of the regression is lost on one set while remains in the other. Q1: The above method is valid, isn't it? Q2: The level of confidence of my research is set at 95%. For one group, the regression is significant at .000. For the other, it is significant at 0.038 So, I believe I have to accept both regressions as significant and that there's no moderating effect. By accepting the regression is significant at while it's proven to be not at 0.01 am I causing a Type I error (accepting the falsy argument)?


title: Interpretation of the Granger causality test

text: I would like to clarify how the Granger causality can/should be used in practice, and how to interpret the statistical significance given by the test. Also, I would like to fill this table with things like "we don't know" or if we know something, what do we know (It will for sure not be causality, but maybe something else?). $\begin{matrix} & X granger cause Y significant & X granger cause Y non significant\\\ Y granger cause X significant & ... & ...\\\ Y granger cause X non significant & ... & ... \end{matrix} $


title: Robust OLS verus ML with sandwich estimator

text: If you compare the standard errors of the OLS coefficients with the White correction, versus the ML estimates with the variance estimated with the sandwich estimator, which standard errors do you expect to be the biggest? And is it correct that both are robust to heteroscedasticity?


title: How to choose training and test sets

text: I would like to propose a single model (decision tree), that is very variable, and validate it. I have choosen parameters after I had obtained good quality measures with a cross-validation. I could build the model on the whole data set and show cross-validated measures. But I can't get a special graph (called Reliability Plot) specific for that model. I should split my data set in training and test sets to obtain that specific graph. The model builded on the training set is different from the optimezed on the whole dataset. Could I choose my training set (50% of the total) to obtain the same model as the whole data set builded one? There is something unwise or wrong in this method? Thanks


title: Question on Exponential Distribution

text: I have a question on the exponential distribution: **Question** > The waiting time to receive food after placing an order at the local > sandwich shop follows an exponential distribution with a mean of 70 seconds. > 80 percent of the patrons wait less than how many minutes to receive their > food? **Attempt** > My attempt was to assume $\lambda$ as $7/6$ minute per order. However, I am > clueless about how to move on. I attempted to use the equation $F(x) = 1 - > e^{-\lambda x}$ with result to be $0.61$, which differs greatly from the > model answer $1.878$ minutes. Appreciate some guidiance as to what I am not doing right please. Is the equation used even correct and if $x$ should be $0.8$?


title: Dependent t-test - multiple readings

text: Let's suppose that I'd like to apply dependent t-test to assess if there is a significant difference in the value of the dependent variable between a group of subjects $X_{1}$ and a related group of subjects $X_{2}$. The same subjects are present in both $X_{1}$ and $X_{2}$. Now, the twist is that we have multiple readings of the dependent variable for each subject in each group. Can dependent t-test be applied in this setting? (for instance, by considering the mean of the dependent variable for each subject) Or should I use a different kind of test?


title: What is the meaning of a gold standard?

text: While reading a few papers, I came across the term "gold set" or "gold standard". What I don't understand is what makes a dataset gold standard? Peer acceptance, citation count and if its the liberty of the researcher and the relevance to problem he is attacking?


title: Multi level mixed-effects model in glmer

text: I have data for a set of users (UID), over numerous dates (date), and over numerous sessions (these are just small chunks of time that a user was active) within each date. For each session I have a set of features (Xi) calculated in which I would like to estimate the gender of the user which is a binary variable (Y). Note that each user does not have the same set of dates, or the same set of sessions within each date. An example of the dataset is as follows:                      UID  Date  Session  X1   X2   Y          1  AAA  1/1   2        3.4  2.2  0         2  AAA  1/1   4        2.5  4.2  0         3  AAA  1/1   5        1.8  4.8  0         4  AAA  1/1   8        6.5  2.9  0         5  AAA  1/2   1        3.5  8.0  0         6  AAA  1/2   2        7.3  5.1  0         7  AAA  1/2   5        1.7  3.5  0         8  BBB  1/1   1        3.7  1.3  1         9  BBB  1/1   2        8.2  9.4  1         10 BBB  1/4   4        2.7  7.3  1         11 BBB  1/4   7        5.2  6.5  1      I am trying to fit this dataset to a generalized mixed-model in glmer but I am stuck on writing the model. Below are a couple of models I am trying but have no idea which (if any) are correct. Notes (assumptions which may be wrong): 1) I believe that since session is the highest level of grouping I do not actually need to write that into any group. 2) I believe UID is a fixed effect since each UID has an underlying true gender and therefore is fixed. I believe data and session are random effects since these are actually changing with time. Models:                   m1 <- glmer(Y ~ X1 + X2 + (session|UID/Date))         m2 <- glmer(Y ~ X1 + X2 + (session|Date) + (1/UID))         m3 <- glmer(Y ~ X1 + X2 + (session|Date) + (Data/UID))      I think m1 may be correct but I am not sure. Thank you for any help that you can provide.


title: How to plot multiple loess curves in a xyplot panel data?

text: Imagine the following simplified longitudinal data with 4 individuals and 3 variables (y1, y2, y3) measured at 50 occasions :               dat<-data.frame(subject=as.factor(rep(1:4, each=50)),time=rep(1:50,4),y1=rnorm(200,50,20),y2=rnorm(200,70,20),y3=rnorm(200,20,20))      It is easy to plot at the same time the raw time series with this code :               xyplot(y1+y2+y3~time|subject,dat,type='l')      which result is : ![enter image description here](http://i.stack.imgur.com/7aWcy.jpg) I would like to plot exactly the same, but with loess curves instead of lines. And here is my problem! The following code :               xyplot(y1+y2+y3~time|subject,dat,     panel=function(x,y)         {panel.loess(x,y,span=0.2)     })      doesn't work because it considers y1+y2+y3 as a single variable, thus plotting only one loess curve by panel. How is it possible to plot each loess curve in each panel? PS. I don't want to use `type="smooth"` because I need to specify the span of the loess.


title: How to specify Bayesian mixed effects model in BUGS

text: I posted this earlier in the week then retracted the question when I found a good source, not wanting to waste people's time. I haven't made much progress I'm afraid. In trying to be a good citizen here I will make the problem as clear as possible. I suspect there will be few takers. I have a dataframe in R I wish to analyse in BUGS or R. It is in long format. It consists of multiple observations on 120 individuals, with a total of 885 rows. I am examining the occurrence of a categorical outcome - but that's not really relevant here. The question is about something deeper. The model I have been using up to here is                 mymodel<-gee(Category ~ Predictor 1 + Predictor 2..family=binomial(link="logit"),       data=mydata,         id=Person)      with a marginal model essentially accounting for the clustering of patients. I then examined                mymodel<-gee(Category ~ Predictor 1 + Predictor 2.. , family=binomial(link="logit"),        corstr = "AR-M",        data=mydata, id=Person)      in order to account for the time ordering of the observations on the individual people. This didn't change much. Then I tried to model them using the following set of MCMCPack commands:                mymodel<-MCMCglmm(category~  Predictor1 + Predictor2..,        data=mydata, family=binomial(link="logit"))      An examination of the output was thrilling, showing statistical significance for many predictors. I hailed myself as a newly converted bayesian, until I realised I hadn't accounted for repeated measures within patients. I understand that I have to account for that. I understand that this may mean fitting a hyperprior for each individual - is that right? What form will this take in BUGS? Here's a basic log reg model: (kudos to Kruschke, J., Indiana)                   model {       for( i in 1 : nData ) {         y[i] ~ dbern( mu[i] )         mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))       }       b0 ~ dnorm( 0 , 1.0E-12 )        for ( j in 1 : nPredictors ) {         b[j] ~ dnorm( 0 , 1.0E-12 )       }     }      However, no hyperprior here for the individual. Here's my best attempt so far at a within-individual design, accounting for repeated measures within people: Here's Jackman's model for JAGS               1 model{     2 ## loop over data for likelihood     3 for(i in 1:n){     4  y[i] ~ dbern( mu[i] )         mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))     6 }     7 sigma ˜ dunif(0,20) ## prior on standard deviation     8 tau <- pow(sigma,-2) ## convert to precision     9     10 ## hierarchical model for each state’s intercept & slope     11 for(p in 1:50){     12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal     13 }     14     15 ## means, hyper-parameters     16 for(q in 1:2){     17 mu[q] ˜ dnorm(0,.0016)      } Here's my bastard-child model for BUGS               1 model{     2 ## loop over data for likelihood     3 for(i in 1:n){     4 mu.y[i] <- alpha + beta[s[i],1] + beta[s[i],2]*(j[i]-jbar)     5 demVote[i] ˜ dnorm(mu.y[i],tau)     6 }     7 sigma ˜ dunif(0,20) ## prior on standard deviation     8 tau <- pow(sigma,-2) ## convert to precision     9     10 ## hierarchical model for each state’s intercept & slope     11 for(p in 1:120){     12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal     13 }     14     15 ## means, hyper-parameters     16 for(q in 1:2){     17 mu[q] ˜ dnorm(0,.0016)       }      Can somebody let me know if I'm heading in the right direction. My understanding of this is growing, but slowly. Please be gentle. I'm a medic, not a statistic! I have used R quite a bit, but I'm new to BUGS, and new to Bayes. Thanks, R


title: How to detect structural change in a timeseries

text: Is there a specific method to detect change points(structural breaks) in a timeseries?(stocks prices) Thanks


title: Effective sample size for test of independence

text: I am interested in applying a test of independence such as chi-square or log- likelihood in the case where the observations are not independent (observations are sequential in time and there is a temporal correlation). I have read about computing the 'effective sample size' for use in variance estimates of dependent data, and intuitively it seems like the same idea should apply to these tests (since I have fewer indepedent samples than actual samples). I have struggled to derive anything concrete though. Does anyone have ideas about this or pointers to existing related work?


title: Generate random strings based on regular expressions in R

text: I'm trying to generate random textual data based on regular expressions. I'd like to be able to do this in R, as I know that R does have regex capabilities. Any leads? This question has come up before in forums (StackOverflow Post 1, StackOverflow Post 2, etc.), but they always mention solutions based on other programming platforms (perl, .NET, ...), not R.


title: How to calculate confidence interval when only a part of the samples are valid?

text: I will simplify our problem in this way. Say there are 100,000 cases in total to examine. Due to the time limitation, we randomly selected 2,000 of them. Then we found 1,000 of them are invalid, so we have only 1,000 valid cases left. Finally we categorize these 1,000 valid cases into 2 categories, A and B; they have 300 and 700 cases respectively. We want to calculate the confidence interval, in order to check the statistical significance of the results. In other words, if the result shows that there are 30% of the cases in category A, how trustworthy this percentage is when talking about the whole population. We used this website, http://www.surveysystem.com/sscalc.htm, to calculate the confidence interval, so the percentages will be like $30\pm 3.7\%$ and $70\pm 3.7\%$. So there are two ways of deciding the population and the sample size. (1) Population is 100,000 (real population) and sample size is 2,000; (2) Population is 50,000 (estimated) and sample size is 1,000. I think we should use option (2), because we actually found _roughly_ only 50% of the original cases are valid cases. But this ratio is actually _estimated_ by the 2,000 cases we sampled. How does this estimation affect the confidence of the result? Would somebody recommend other ways to check the statistical significance of the result in our case? Thanks for your help!


title: Detecting patterns of cheating on a multi-question exam

text: QUESTION: I have binary data on exam questions (correct/incorrect). Some individuals might have had prior access to a subset of questions and their correct answers. I don’t know who, how many, or which. If there were no cheating, suppose I would model the probability of a correct response for item $i$ as $logit((p_i = 1 | z)) = \beta_i + z$, where $\beta_i$ represents question difficulty and $z$ is the individual’s latent ability. This is a very simple item response model that can be estimated with functions like ltm’s rasch() in R. In addition to the estimates $\hat{z}_j$ (where $j$ indexes individuals) of the latent variable, I have access to separate estimates $\hat{q}_j$ of the same latent variable which were derived from another dataset in which cheating was not possible. The goal is to identify individuals who likely cheated and the items they cheated on. What are some approaches you might take? In addition to the raw data, $\hat{\beta}_i$, $\hat{z}_j$, and $\hat{q}_j$ are all available, although the first two will have some bias due to cheating. Ideally, the solution would come in the form of probabilistic clustering/classification, although this is not necessary. Practical ideas are highly welcomed as are formal approaches. So far, I have compared the correlation of question scores for pairs of individuals with higher vs. lower $\hat{q}_j -\hat{z}_j $ scores (where $\hat{q}_j - \hat{z}_j $ is a rough index of the probability that they cheated). For example, I sorted individuals by $\hat{q}_j - \hat{z}_j $ and then plotted the correlation of successive pairs of individuals’ question scores. I also tried plotting the mean correlation of scores for individuals whose $\hat{q}_j - \hat{z}_j $ values were greater than the $n^{th}$ quantile of $\hat{q}_j - \hat{z}_j $, as a function of $n$. No obvious patterns for either approach. * * * UPDATE: I ended up combining ideas from @SheldonCooper and the helpful Freakonomics paper that @whuber pointed me toward. _Other ideas/comments/criticisms welcome._ Let $X_{ij}$ be person $j$’s binary score on question $i$. Estimate the item response model $$logit(Pr(X_{ij} = 1 | z_j) = \beta_i + z_j,$$ where $\beta_i$ is the item’s easiness parameter and $z_j$ is a latent ability variable. (A more complicated model can be substituted; I’m using a 2PL in my application). As I mentioned in my original post, I have estimates $\hat{q_j } $ of the ability variable from a separate dataset $\\{y_{ij}\\}$ (different items, same persons) on which cheating was not possible. Specifically, $\hat{q_j} $ are empirical Bayes estimates from the same item response model as above. The probability of the observed score $x_{ij}$, conditional on item easiness and person ability, can be written $$p_{ij} = Pr(X_{ij} = x_{ij} | \hat{\beta_i }, \hat{q_j }) = P_{ij}(\hat{\beta_i }, \hat{q_j })^{x_{ij}} (1 - P_{ij}(\hat{\beta_i }, \hat{q_j }))^{1-x_{ij}},$$ where $P_{ij}(\hat{\beta_i }, \hat{q_j }) = ilogit(\hat{\beta_i} + \hat{q_j})$ is the predicted probability of a correct response, and $ilogit$ is the inverse logit. Then, conditional on item and person characteristics, the joint probability that person $j$ has the observations $x_j$ is $$p_j = \prod_i p_{ij},$$ and similarly, the joint probability that item $i$ has the observations $x_i$ is $$p_i = \prod_j p_{ij}.$$ Persons with the lowest $p_j$ values are those whose observed scores are conditionally least likely -- they are possibly cheaters. Items with the lowest $p_j$ values are those which are conditionally least likely -- they are the possible leaked/shared items. This approach relies on the assumptions that the models are correct and that person $j$’s scores are uncorrelated conditional on person and item characteristics. A violation of the second assumption isn’t problematic though, as long as the degree of correlation does not vary across persons, and the model for $p_{ij}$ could easily be improved (e.g., by adding additional person or item characteristics). An additional step I tried is to take r% of the least likely persons (i.e. persons with the lowest r% of sorted p_j values), compute the mean distance between their observed scores x_j (which should be correlated for persons with low r, who are possible cheaters), and plot it for r = 0.001, 0.002, ..., 1.000. The mean distance increases for r = 0.001 to r = 0.025, reaches a maximum, and then declines slowly to a minimum at r = 1. Not exactly what I was hoping for.


title: Different number of samples (observations) per class (one vs. all classification)

text: I have 20 observations labeled for each class. The number of classes is 5, so I have a total of 100 observations. I want to classify one class vs. other classes (one vs. all). For this, I first labeled '1' for one of five classes (= 20 observations), and '2' for the other classes (= 80 observations). After that, I spiltted the whole data into training data set and test data set for each class ('1' or '2'). Since I used 10-fold cross-validation method, 20 observations for class '1' were separated into 18 and 2 for training data set and test data set, respectively (62 training data and 18 test data for class '2'). I trained a LDA classifier for '1' or '2' labeled observations using training data set, and estimated classification accuracy using test data set. My question is that isn't there any bias problem when using different number of observations to train a classifier. In my case, the difference is quite big (20 vs 80). I am worrying about this issue. If so, please tell me a proper solution (e.g., a specific classifier) and reference paper to overcome this problem.


title: How to calculate the number of trials for a 'meaningful' result

text: Apologies for newbie question but have failed to understand what I've read about in this area so pointers to the correct methods would focus my attention. I have a piece of software which when tested repeatedly in the same way fails occasionally - say once in 20 runs. If I make a change to fix the problem, how many times should I test it to be confident (say 95%) that the changed software is actually better than the original. More generally: if a test fails a out of b times. After the **fix** , how many times (c) must it run without any failures to be confident to d%.