title: Generating random variables satisfying constraints

text: I need to generate a list of random variables $\bf{x}$ subject to constraints that can be expressed in the form $\bf{E}x=b$ where $\bf{E}$ is an $m \times n $ matrix if $\bf{x}$ has $n$ entries. In all the cases I'm dealing with, $n >> m$, for example $n$ will be around 14,000 and $m$ will be 50. I'm not sure what method I will use for random sampling, either normal or uniform, it is not clear which is best for the problem I'm trying to solve, but I need all the variables to be sampled from distributions with the same mean and range/variance. What I've been doing to solve this is reducing $\bf{E}$ to row-echelon form, setting all the variables corresponding to columns to the right of the last pivot to random values, and then solving the remaining square matrix equality. There is a problem however, to solve the square matrix equality, I subtract the already set values from the right hand side. Unfortunately, the variances add as well, so my last 50 values tend to vary hugely, which is unfortunately unacceptable in this problem. Is there a better way to do this? I cannot think of a way to fix the current method I am using. I use R.


title: What is the difference between a stationary test and a unit root test?

text: What is the difference between the Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test and the augmented Dickey-Fuller (ADF) test? Are they testing the same thing? Or do we need to use them in different situations?


title: Understanding scipy Kolmogorov-Smirnov test

text: I'm trying to understand the Kolmogorov-Smirnov test using a very simple example. I generate a set of random, uniform values between 0 and 1.0. I then test that these values are from a uniform distribution by using the scipy kstest function. I'm expecing a very small D value and a pvalue close to 1.0, but instead I get wildly varying pvalues every time I run the code. What am I missing?               import numpy as np     import scipy     a = np.random.uniform(size=4999)     print(scipy.stats.kstest(a, 'uniform'))      Here are the outputs of a few consecutive runs:               (0.0075523161200627964, 0.93798952050647577)     (0.013787195268362473, 0.29799260741344774)     (0.014359046616557847, 0.25402403230845855)     (0.012521820948675988, 0.41329007558099806)     (0.011159003477582918, 0.56216895575676396)


title: Understanding why p value is uniformly distributed

text: I know that the mathematically, p values should be uniformly distributed under the null hypothesis. However, what is this null hypothesis (in this context)? Can someone give me an example please ?


title: Estimability and less than full rank model matrices

text: So I have some data...               > head(df)               y  a   b   c  d   e  f     g     1  9.393915 No  No  No No Yes No WAW9H     2  9.936906 No  No  No No Yes No WAW9H     3  7.573025 No  No Yes No  No No WAW9H     4  9.123976 No Yes Yes No Yes No WAW9H     5 13.569253 No Yes  No No Yes No WAW9H     6 10.468927 No  No  No No Yes No WAW9H      and I fit this model... > mod <\- lm(y ~ a + b + c + d + e + f + g + a:g + b:g + c:g + d:g + e:g + > f:g, data = df) but some of the coefficients are NA:               > summary(mod)          Call:     lm(formula = y ~ a + b + c + d + e + f + g + a:g + b:g + c:g +      d:g + e:g + f:g, data = df)          Residuals:      Min       1Q   Median       3Q      Max      -16.2662  -1.8112   0.0096   1.8455  28.3563           Coefficients: (55 not defined because of singularities)                   Estimate Std. Error t value Pr(>|t|)         (Intercept)   9.855044   0.128512  76.686  < 2e-16 ***     aYes          0.152247   0.341803   0.445 0.656016         bYes          1.536731   0.156493   9.820  < 2e-16 ***     cYes          0.956878   0.130900   7.310 2.71e-13 ***     .     .     .     aYes:gBJK2U         NA         NA      NA       NA       I suppose it has to do with there not being enough data to estimate these coefficients (check out BJK2U)               > table(data.frame(df$a, df$g))     df.a  AS2CD BJK2U . . .       No   2096     2 . . .           Yes    84     0 . . .          How do I fit a reduced model without these coefficients? How do I remove some of the levels of the factor g? The reason I want to do this is so that I don't get warnings when I call the predict() function on this model. I don't want to accidentally try to predict a non-estimable thing. Apologies for the lack of reproducibility in this example.               > levels(df$g)      [1] "AS2CD" "BJK2U" "CAT2N" "CER"   "CG3N"  "CH2XF" "CMU5P" "CV7D"  "DDS2S" "DGP3A" "DRB5P" "DRG3E" "DTL4N" "EAR9W" "ES6T"  "GEM3M" "GRU6N" "JHM7Q" "JJA9H"     [20] "JPR5K" "JSH5N" "JSR7N" "MAK4Z" "MKW2Y" "MS7HA" "NPC8A" "PAR8B" "PRK5R" "RCR3E" "RDP5P" "REO4X" "RM5S"  "RR5WF" "RS5FG" "SAO4R" "SAS3K" "SDL4V" "SK7BQ"     [39] "TAS8N" "WAW9H" "WB4Z"


title: Statistics question on interpreting my results for sampling distribution

text: The weights of pennies minted after 1982 are approximately normally distributed with mean 2.46 grams and standard deviation 0.02 grams. What is the probability that in a simple random sample of 10 pennies minted after 1982, we obtain a sample mean of at least 2.465 grams? I found out that the answer is 0.2148, but I don't know how to interpret my result. I know that if I take 100 of these pennies about 21 of them will result in a sample mean of at least 2.465, but when I say this it makes little sense. Am I saying that I am averaging 21 of these pennies to be at least 2.465 grams? Are they averaged randomly? What's a better interpretation? The probability that a randomly selected penny will have a sample mean of 2.465 grams is 21.48%? When I say this, I still can't make sense of the wording.


title: Predicting trend component using theta method

text: I have a gold price series data. I want to make a forecast of gold price series. First, the series is decomposed with STL, and then each component is predicted using GRNN and theta method. The theta method will be used to forecast the trend component. If I want to use thetaf() function, can I apply the function directly to the trend component to get the forecast of the trend component? Thank you.


title: Sufficiency of order statistics

text: I am told the following proof is incorrect, but I cannot understand why. Consider $X_{(1)}, \ldots, X_{(n)}$ are the order statistics of a random sample of size $n$. I want to show that the order statistics are sufficient. So I wrote down: $$P(X_1, \ldots, X_n|X_{(1)}, \ldots, X_{(n)})= \tfrac{1}{n!}$$ as given the vector of order statistics, there are $n!$ possibilities for the sample $X_1, \ldots, X_n$. As we are in an i.i.d. case, then each vector is equiprobable, and so the equality follows. I am told this is not true specially in the case of discrete random variables. I don't see how it is wrong though. Any explanation would be great.


title: Is there an efficient way to compute the ordered probabilities of a weighted sample of events?

text: For this example, I have 4 different events **a ... d** , each with a different probability of occurrence. My sample space is the following: { **a a b b b b c c c d d d** } p( **a** ) = 2/12 p( **b** ) = 4/12 p( **c** ) = 3/12 p( **d** ) = 3/12 Samples are drawn one at a time without replacement. As soon as a particular event is drawn, all other identical outcomes are removed from the sample space. For example, if I draw event **b** , then the other 3 **b** 's are removed from the sample space leaving me with { **a a c c c d d d** }. Is there an efficient way to calculate the probability of any particular event occurring in any particular position? In this example, I have 4 distinct events that can be drawn and 4 possible positions for each event. I'd like to be able to compute the probability of drawing event **a** in the 1st, 2nd, 3rd, and 4th positions, event **b** in the 1st, 2nd, 3rd, and 4th positions, and so forth. Since the likelihood of a particular event occurring in particular position depends on the events selected before it, I am unsure how to come up with a formula that generalizes to larger datasets. I can compute this fairly straightforwardly with basic conditional probability in this example, but I'm wondering if there is a more efficient formulation out there when the numbers of possible events increase beyond a reasonable standard for manual calculation. Any help with this would be greatly appreciated. Thanks!


title: Missing values in GLMM

text: I am creating a GLMM based on a experiment where each subject has 2 repeats. In some instances though there is only data for one of a given subjects repeats for most there is data for both. Can I still justify having subject as a random effect? Thanks, Jonathan


title: What's a stationary VAR?

text: * What is a stationary VAR (vector autoregression)?   * Can a VAR with non-stationary variables be stationary?   * How do you test whether a VAR is stationary or non-stationary? (Example in `R` language if possible/applicable).


title: How can I get annual rates of change for combined trend estimates?

text: I would like to combine trend indices (gained with different methods referring to the same subject, assuming they do not differ significantly) of two different time series and to derive the combined annual rate of change and standard error. To combine annual indices ($Y$) and standard errors ($S$) in the overlapping period I use the following formulas: $$ \begin{aligned} Y^*&=(Y_1/S_1^2+Y_2/S_2^2)\cdot S_1^2\cdot S_2^2/(S_1^2+S_2^2) \\\ S^*&=\sqrt{S_1^2+S_2^2/S_1^2\cdot S_2^2} \end{aligned} $$ Here is an example dataset with trend estimates gained from generalized estimation equations (I do not have the raw data from which the estimates are calculated):               year      Y1        S1        Y2        S2     1990      91.3      12.0      NA        NA      1991      81.8      8.7       NA        NA      1992      65.9      6.9       NA        NA      1993      80.8      8.0       NA        NA      1994      81.7      8.1       NA        NA      1995      90.0      8.5       NA        NA      1996      113.7     10.0      NA        NA      1997      117.5     10.2      NA        NA      1998      139.9     11.6      NA        NA      1999      120.0     9.9       NA        NA      2000      130.1     10.6      NA        NA      2001      115.1     9.4       NA        NA      2002      105.2     9.2       NA        NA      2003      83.9      7.5       NA        NA      2004      105.8     8.9       NA        NA      2005      111.0     9.5       85.7      7.2         2006      100.0     0         100.0     0           2007      105.3     9.4       104.0     6.1         2008      93.5      8.6       102.3     5.5         2009      117.0     10.8      111.5     6.2         2010      107.0     10.4      108.5     6.1          Now, I would like to get an estimate of the combined annual rate of change and its standard error for the time period of 1998 to 2008. Does anyone have any idea to share?


title: multivariate dirichlet for multiple imputation

text: I dealing with 3 covariates {x1, x2, x3} all three are discrete and contain missing data.                                  x1  | x2| x3 |                        --------------                        3   | 2 |  . |                        .   | 1 |  3 |                        1   | . |  2 |                        3   | 2 |  . |                        2   | . |  1 |                        .   | 3 |  1 |                        --------------      I am exploring the option of using multiple imputation. I plan on generating observations using joint multivariate dirichlet distribution. Am i heading in the right direction ? Please advise. Below is my whole dataset ![enter image description here](http://i.stack.imgur.com/NE2Z3.png) / ********************** pi***************************/               y   x1  x2  x3  Freq    pi     0   1   1   1   15  0.02078     1   1   1   1   7   0.0097     0   2   1   1   13  0.01801     1   2   1   1   9   0.01247     0   3   1   1   10  0.01385     1   3   1   1   13  0.01801     0   1   2   1   16  0.02216     1   1   2   1   11  0.01524     0   2   2   1   14  0.01939     1   2   2   1   18  0.02493     0   3   2   1   6   0.00831     1   3   2   1   26  0.03601     0   1   3   1   13  0.01801     1   1   3   1   16  0.02216     0   2   3   1   18  0.02493     1   2   3   1   16  0.02216     0   3   3   1   17  0.02355     1   3   3   1   15  0.02078     0   1   1   2   12  0.01662     1   1   1   2   18  0.02493     0   2   1   2   8   0.01108     1   2   1   2   17  0.02355     0   3   1   2   16  0.02216     1   3   1   2   11  0.01524     0   1   2   2   17  0.02355     1   1   2   2   10  0.01385     0   2   2   2   8   0.01108     1   2   2   2   14  0.01939     0   3   2   2   17  0.02355     1   3   2   2   11  0.01524     0   1   3   2   11  0.01524     1   1   3   2   13  0.01801     0   2   3   2   16  0.02216     1   2   3   2   23  0.03186     0   3   3   2   11  0.01524     1   3   3   2   21  0.02909     0   1   1   3   15  0.02078     1   1   1   3   8   0.01108     0   2   1   3   12  0.01662     1   2   1   3   11  0.01524     0   3   1   3   10  0.01385     1   3   1   3   11  0.01524     0   1   2   3   8   0.01108     1   1   2   3   13  0.01801     0   2   2   3   23  0.03186     1   2   2   3   11  0.01524     0   3   2   3   7   0.0097     1   3   2   3   14  0.01939     0   1   3   3   8   0.01108     1   1   3   3   5   0.00693     0   2   3   3   8   0.01108     1   2   3   3   16  0.02216     0   3   3   3   16  0.02216     1   3   3   3   19  0.02632


title: How should I do a grid search for the gamma in SVM?

text: I know that you can do a grid search for the C parameter (-c) in libSVM by going through value that go from 10^-5 to 10^5. How should I go about finding the optimal epsilon parameter (-p) ? Is performing a grid search from the default p = 0.1, in different powers of 10 the best way to go? Is it critically kernel dependent? (I think so, since for a linear kernel an epsilon value might give a different result than an RBF). To give some context, my dataset is image feature related in the domain of computer vision, so feature spaces oscillate between 512 to 4096 in dimension per exemplar.


title: What does the standard error of my IV estimate tell me?

text: I have computed the IV estimate and the standard error of the IV estimate. Why do I care about its standard error?


title: What if results from a Breusch-Pagan test for heteroscedasticity contradicts those of a White's test?

text: Testing for heteroscedasticity I get these results: Breusch-Pagan / Cook-Weisberg test for heteroskedasticity Ho: Constant variance Variables: fitted values of log_expdu                    chi2(1)      =    21.41          Prob > chi2  =   0.0000      White's test for Ho: homoskedasticity against Ha: unrestricted heteroskedasticity                    chi2(47)     =     48.91          Prob > chi2  =    0.3964      The Breusch-Pagan test tells me that there is heteroskedasticity, while the White's one the opposite. Which result should I use?


title: Compute expected standard deviation from proportion data

text: This example is from the first chapter of Gelman & Hill 2007: > > The following are the proportions of girl births in Vienna for each month > in 1908 and 1909 (out of an average of 3900 births per month): >> >> .4777 .4875 .4859 .4754 .4874 .4864 .4813 .4787 .4895 .4797 .4876 .4859 .4857 .4907 .5010 .4903 .4860 .4911 .4871 .4725 .4822 .4870 .4823 .4973 >> >> The data are in the folder girls. von Mises (1957) used these proportions to claim that the sex ratios were less variable than would be expected by chance. >> >> Compute the standard deviation of these proportions and compare to the standard deviation that would be expected if the sexes of babies were inde- pendently decided with a constant probability over the 24-month period. * * * The standard deviation is: 0.064 I'm wondering how to calculate the expected standard deviation? Should I use a binomial model? but where do the probabilities come from then? EDIT: The next step would be to assess whether the observed (0.0064) and chance (0.0080) deviations significantly differ. The exercise gives a hint: "under the randomness model, the actual variance should have a distribution with expected value equal to the theoretical variance, and proportional to a chi^2 with 23 degrees of freedom" I'm wondering how I could assess whether the two deviations are significantly different. The only solution I can come up with is (observed-expected^2)/expected = (0.0064-0.0080)^2/0.0080 = 0.00032 and then compare this to a chi^2 distribution with 23 degrees of freedom. Is this the correct way to do it?


title: Distribution of starting points for blocks within DNA

text: Given a DNA sequence of $L$ I am looking for a distribution to model the starting positions of $k$ blocks with length $w$. I am starting to investigate in a discrete form of dirichlet distribution but I was curious if any of you ever faced the same problem and how you solved it. So, any ideas?


title: How to compute simulated differences of means repetitively?

text: I would like to run a 1000 simulations of $X-Y$ where $X$ is the average of 100 random observations with mean of 69.5 and SD 2.9 and $Y$ is the average of 100 random observations with mean of 63.9 and SD 2.7. I did the following:               x <- rnorm(100, 69.5, 2.9)     y <- rnorm(100, 63.9, 2.7)      And now I would like to somehow run a 1000 simulations of $X-Y$.


title: performing parameter updates on Gamma distribution

text: I have the following form for a joint distribution $$ P(w, \lambda, \phi \vert y) = P(\phi) \times P(w \vert \lambda) \times P(\lambda) \times \prod_{i=1}^{N}P(y_i \vert w_i, \phi, \lambda) $$ I have some approximation for $w$ and what I am trying to do is update the parameters for the distributions for $\lambda$ and $\phi$. To take advantage of the conjugate prior properties, I have the following setup. $\lambda$ is modelled as a Gamma distribution where $a$ and $b$ are the scale and rate parameters of the distribution respectively. $$ P(\lambda; a, b) = \frac{b^a \lambda^{a-1} \exp{(-b\lambda)}}{\Gamma(a)} $$ The log-posterior on $\lambda$ can be then written as: $$ \log P(\lambda; a, b) = (a-1)\log(\lambda) -b\lambda + const(\lambda) $$ $\phi$ which is the noise _precision_ is also modelled as a Gamma distribution and the log posterior can be written similarly as: $$ \log P(\phi; s, t) = (s-1)\log(\phi) -t\phi + const(\phi) $$ The log transformed prior on $w$ (which depends on $\lambda$ is: $$ \log P(w \vert \lambda) = \frac{1}{2}(\log(\lambda \Lambda) - w^T(\lambda \Lambda)w) + const (\lambda, w) $$ This is just a zero mean MVN with the _precision_ being defined as $\lambda \Lambda$. Finally, the log-likelihood term can be written as: $$ \log P(y \vert x, w, \phi) = \frac{N}{2}\log(\phi) -\frac{\phi}{2} \sum_{i=1}^{N}(y_i - t(x_i, w_i)) + const(w, \phi) $$ x and y are observed and t is a non-linear transformation parameterised by $w$, which is the parameter of interest. So, what I wanted to do was for a given $w$ obtain the posterior marginal distribution for $\lambda$. Now my, log posterior is given as (ignoring constant terms for the moment): $$ \mathcal{L} = \log P(\lambda; a, b) + \log P(\phi; s, t) + \log P(w \vert \lambda) + \log P(y \vert x, w, \phi) $$ So, I thought that we can condition on $w$ treating it as a constant and integrate over $\phi$ and we should get the parameters for $\lambda$ as the posterior should also be Gamma distributed due to conjugate properties. So, the posterior density on $\lambda$ is: $$ P(\lambda) = \mathcal{L}\:d\phi $$ Now, doing the integrating I get something complicated as: $$ \phi((a-1)\log\lambda - b\lambda) + (s-1)\phi(log \phi - 1) -\frac{t\phi^2}{2} + \phi (\frac{1}{2}(\log(\lambda \Lambda) - w^T(\lambda \Lambda)w)) + \frac{N}{2}\phi(\log\phi -1) -\frac{\phi^2}{2}(y-t(x, w))^T(y-t(x, w)) $$ So basically it becomes something complicated and it is not clear to me at all how the posterior Gamma distribution can be inferred from it which makes me think I have made a mistake with my update strategy. Can someone comment whether updating like this should be feasible? i thought it should be as this is whole point of using conjugate priors but I am at a loss as to what I have done wrong. I am inclined to think that the problem lies in the step $$ P(\lambda) = \mathcal{L}\:d\phi $$ but I am not sure what is the main issue.


title: Calculating significance and uplift on Revenue A/B Tests

text: I'm currently running an AB test on a bunch of websites for a change we've made on all the websites and am measuring revenue/site/group so my test groups look like so: Website a | Test group a | unique visitors | revenue Website a | Test group b | unique visitors | revenue Website b | Test group a | unique visitors | revenue Website b | Test group b | unique visitors | revenue e.t.c. Initially I want to show there's a difference between the two test groups overall and so have used a null hypothesis of: "The test group makes no difference to the revenue/cookie across all sites". I can then happily sum all the unique visitors (across sites) and revenue for each test group. However, I don't really know what to do from here. I know that conversions can be treated binomially but revenue can't be. We also can't treat average order value as an exponential or normal distribution (ideas I looked at). Initially, I thought about using a Mann-Whitney U test but most things I've read seems to suggest these work best for < 20 entries in your array. When I create a vector containing the amount every customer spent (a very sparsely populated vector - most people don't buy anything), each test group has a different length vector of size ~100,000. My next attempt involved bootstrapping (I think) - I basically created a normal distribution for each test group by sampling from the aforementioned sparsely populated revenue vector. Creating around 10,000 vectors like this and running a normality test on the sum of each of the vectors gives me a normal distribution for each test group. From this I believe I can perform a t-test on the two normally distributed vectors to find the significance of the change made and hence support/disprove the null hypothesis? My question is, are either of the above attempts the right thing to do? I can't imagine that this is that uncommon a thing to want to do that there isn't a generally accepted way of doing this. Additionally, I'd also like to be able to report something along the lines of: Test group x is y% better than test group z with 5th and 95th percentile on that percentage uplift of q and j. What would be the best way to calculate these bits of information? Thanks, Matt


title: Fitting a log (or generalized?) linear model

text: I would like to fit a model of the form $z = k x^\alpha y^\beta$ to some data I have (it's a spatial gravity model). Now I know you could take logs of both sides and fit a linear regression $\log z = \log k + \alpha \log x + \beta \log y$ But there are two problems:   1. Does this violate assumptions of OLS regression if $\log x$ and $\log y$ don't have a normal distribution?   2. Some of my $z$s are zero - how to deal with that? **EDIT** In response to whuber's question. This is a model of how many people from each of a number of districts travel to an event. $x$ is district population and $y$ is (mean) distance (from each district) - thus I expect $\alpha \approx 1$ and $\beta < 0$. $z$ is sampled from a survey of under 1% of over 100k visitors. I am interested in obtaining estimates for where the remaining 99% of visitors came from, for which I do not expect $z=0$ for any district, but in the limited sample data $z$ is sometimes zero because the survey people didn't happen to meet anybody from that district. If anyone has comments about how the above modifies expectations about the error distribution on $z$ and hence the best way to fit the model I would be interested to hear it. As I was expecting $a \approx 1$, I tried nonlinear least squares with $a$ fixed at 1, but that led to extremely bad variance on $k$ ($\sigma_k > k$). I have since fitted GLM with a poisson link function which gives sensible $k$ and $\beta$ but $\alpha \approx 0.85$ - not great, but it's only a rough model so I can live with it.


title: Consistency of OLS in presence of deterministic trend

text: For consistency of OLS estimator for linear model $$ y_i = \beta^T x_i + \epsilon_i, \; i = 1,\cdots, n, $$ the model assumptions are usually (the ones I am familiar with)   1. The sequence of random vectors $\\{ (x_i, \epsilon_i) \\}$ is (jointly) strictly stationary and ergodic.   2. $\mathbb{E}[x_i \epsilon_i] = 0$ for all $i$. LLN then gives consistency. Now if one of the regressors, say $(x_i)_1$, is $i$---deterministic linear time trend if $i = t$ and we have a time series, is consistency still possible? I suppose one trick is to run the regression on the first differences $ y_i - y_{i-1} = \beta^T (x_i - x_{i-1}) + (\epsilon_i - \epsilon_{i-1}), \; i = 1,\cdots, n, $ and assume the additional orthogonality condition $\mathbb{E}[\epsilon_i x_{i \pm 1}] = 0$. Are there standard/more sophisticated ways to deal with this? What about more general time trends---cyclic, exponential, etc?


title: Mann-Whitney U test and K-S test with unequal sample sizes

text: I want to compare two distributions, to see if they are significantly different. They represent task time completions (so they range from 1 to around 1000 seconds) in two different months. They are not normally distributed. I want to see if their central tendencies are significantly different (at a first glance the mode, mean and median between the two months seem very close, just 3-4 seconds difference), but also to see if their shapes are similar (again, at a first glance, they look similar). I am currently carrying this analysis with SPSS 20. I have the Mann-Whitney test for testing central tendencies and the Kolmogorov-Smirnov test for the shape of the distribution, (although I have read that the K-S test is an overall comparison test for the distributions). Also, in the first month I have 300,000 observations and in the second month 122,000 observations. So, a lot of data ... but disproportionate. Is this an impediment to running these tests, the fact that the sample sizes are not equal? I ran both Mann-Whitney and K-S and they both seem to reject the null. How much should I trust the results given my sample sizes? Do you suggest any alternative tests? Thanks


title: How to identify variable (from many variables) which is able to discriminate between groups?

text: I currently have a data frame with 98 observations and 107 variables. All of the variables are numeric, but one variable is binary (yes or no). My goal is to determine which correlation and/or variable give the greatest segregation between the yes and no samples. I have been using the pairs () function to do this, but I can only do a few variables at a time. Is there a way to determine which correlation gives the greatest discernment between yes and no? To Clarify - My table is 98 observations and 107 variables, but doing a correlation matrix with the pairs function is not able to fit all of the variables. I have used this function:               pairs(x[70:80], ch=21, bg=c("red","green")[unclass(x$outcome)])


title: What regression formula would best fit this curve?

text: I am using Matlab to try and find a good fit for this curve: ![cumulative distribution](http://i.stack.imgur.com/pF3lz.png) None of the built-in formulas seem to work well. Any suggestions?


title: Detecting 'causality' in Likert-time series data

text: ### _[Note] I've decided to re-write my question for the sake of brevity. The original question can be found below._ Suppose a number of individuals fill in a questionnaire at a multiple number of time points. In other words, for each person $i \in \\{1, \ldots, m\\}$ and for each time point $t \in \\{1,\ldots,T\\}$, we have responses $\mathbf{y}_i^{(t)} = (y_{i1}^{(t)},\ldots,y_{ip}^{(t)})^{T}$ where $p$ is the number of items on the questionnaire. I am interested in finding out which item responses cause/predict other item responses, preferably using network models. I need a model that can do this under the following conditions:   * all responses $y_{ik}$ are in some finite range (they're Likert responses)   * the data may not be stationary   * I would like to account for feedback loops. That is, at time 1 a response to item A might influence an response to item B at time 2, which might influence item A at time 3.   * An estimation of direction is required, i.e. I'm not only interested in correlation, but also direction. * * * ### Original formulation of my question: I'm trying to find a suitable model for a problem I've been asked to investigate. Suppose a number of individuals fill in a questionnaire at a multiple number of time points. In other words, for each person $i \in \\{1, \ldots, m\\}$ and for each time point $t \in \\{1,\ldots,T\\}$, we have responses $\mathbf{y}_i^{(t)} = (y_{i1}^{(t)},\ldots,y_{ip}^{(t)})^{T}$ where $p$ is the number of items on the questionnaire. I am interested in providing these persons with feedback regarding their responses, e.g. the course of outcomes to specific items, and relations between item-responses. In particular, I'm interested in estimating and visualizing a network model which is fitted to this data. The hope is that this network model can indicate which item-responses are predictive for other item-responses. Granger-causality came to mind but it does not seem the best method here. For one, the data is categorical (Likert-scale) and not necessarily stationary. Also, normality usually cannot be assumed. I am not interested in predicting anything, the most important feature is that one get's a sense of "what item-responses `predict' other item-responses?" What would be a proper model for this setting? Thanks in advance.


title: Interpretation of power and detectable difference

text: I know power is the likelihood of correctly rejecting the null hypothesis, but I just want to check that my interpretation of it with the detectable difference is correct. Consider a 2 sample test with the null hypothesis that there is no difference between group means. My understanding is that, by putting a minimum detectable difference of say 5%, with power 80%, into the sample size calculation, then in conjunction with the computed sample size this means that: > if the true difference is 5% then my test will have a likelihood of 80% of > rejecting the null hypothesis. My test will also have a slightly higher > likelihood of rejecting the null hypothesis if the true difference is 6%, > etc. Is this correct?


title: Summation of a product

text: I need to calculate the following expression: $$\sum_{k=1}^N a_k b_k$$ ${a_k}$ and $b_k$ are real positive numbers. N and k are integers. I know the average values of $a_k$ , defined as $\overline {a} = {\sum_{k=1}^N a_k \over N } $ and $b_k$ , defined as $\overline {b} = {\sum_{k=1}^N b_k \over N } $. I also know the standard deviation of $a_k$ and $b_k$, $\sigma_a$ and $\sigma_b$. If only knowing these quantities, I have to make some approximation, I would like to know how much error I am producing with that approximation. $N$ is relatively big. Any help is appreciated.


title: Reference Request: Regularization

text: Lately, a wealth of regularizers have come into being. The area of model selection has generated a great amount of interest. Many times, we would like to control the complexity of the model, but do so also by promoting some kind of desired structure (sparsity, group sparsity, low rank, ...) I was wondering if there are any good literature survey of regularization. Any relevant references will be appreciated. It is hard to keep up with all the existing regularizers (sparsity-inducing regularizers, mixed norm regularizers, matrix norms regularizers, atomic norm regularizers, ...).


title: Assessing role of a count variable in regression... Do you need a zero?

text: Is it acceptable to run a logistic regression on a yes/no DV and include a predictor variable that is a count of the number of times something happened previously, but none of the cases has a zero count? It seems to me you would be testing if more than 1 event matters, but not whether the overall number matters compared to no events. Thanks. To give some context, it's a study to see how interaction with government officials affects future use of government services. If the participants were limited to those who had at least one interaction, could you then use a regression to identify the potential effect of each additional interaction? It would seem to be a better design and more robust results if you had some "zero" participants too. Thoughts? Thanks for the thoughtful responses. It's not my study, I'm assessing someone else's. It just struck me as odd to talk about the effect of something without including the comparison of people who had no experience with X. It may not be the number of times something happened that creates the effect, but rather that it happened at all. you may detect an intensification of the effect with more of X, but without zero, do you know if it would hold if there was no X. And, yes, there are many other issues, as it's apparent the number of previous X is definitely not random. Thanks again


title: Is there a better test?

text: I have data of the following sort: There are numerous individuals. This total set of individuals is classified according to three treatments.The first two treatments each divide the total set of individuals into two groups. When both are applied, the total set is divided into four subgroups. Each individual is measured before and after the third treatment. This last treatment is actually to be a dependent or outcome variable. The change, if any, that occurs before and after is of interest. If it helps, the first treatment is Sex (male and female levels), and the second treatment is genotype with two different levels (two kinds of otherwise the same gene). The third treatment is a surgery. I am interested in the difference in the result before compared after the surgery. And I am interested in whether Sex and Genotype are significant factors. What I have done: I have done so far two things. First, I took the difference in result before and after surgery, and used that as a dependent variable in an ANOVA using Sex and Genotype and Sex*Genotype as factors. Second, I have included the before and after results in the model. The ANOVA factors were Sex, Genotype, and procedure (two levels: before result and after result), plus interactions. Here, each individual is included in the analysis twice: one time for before surgery result, and another time for after surgery result. Is there a better test?


title: Latent Dirichlet Allocation - understanding the posterior

text: I have a problem understanding the posterior for computing LDA, stated in page 7 of Blei (2007). From my point of view, it's not exactly consistent with Bayes' theorem, as described here. Could anyone give me a simple explanation of how this formula was derived? I don't really understand how the $p(\beta, \theta, z, w)$ corresponds to $P(B|A)P(A)$ in Bayes theorem. I will be extremely grateful for any help.


title: Please help I am confused to decide which statistical test?

text: I have a data set, have two interval variables and lots of categorical (ordinal and nominal) sociodemographic variables.   Lets say one of may interval variable is $x$ and the other is $y$, I want to test effect of $x$ and a demographic variable such as gender or educational level on $y$. In other words, how can I compare or find a relationship of $x$ scores with different educational levels and $y$ scores (dependent variable)? I want to say that there is a difference between $x$ scores from high school and $y$ scores. Or I want to say that there is a relationship between $x$ scores from high school and $y$ scores. Or I want to say that x scores from university educational level is related to y scores.


title: Multiple time series methods for trend identification, forecasting etc

text: I have several time series consisting of aggregated macro-economic indices and I am trying to choose one or several appropriate techniques in order to answer a number of questions. First, I need to identify one or several seasonal trends in the indices. Second, I need to be able to identify anomalities or trend breaks in one or several of the indices. Lastly, I need to be able to forecast the indices or some aspects of them, or at least find out whether the indices are forecastable or not. I know there are several linear, parametrical methods that can used in order to answer these questions. Numerous non-linear and non-parametric methods are of course also available. However, I am quite novice when it comes to time series analysis and non-parametric methods specifically. Therefore, I would appreciate if you could give me tips on which method, or mixture of methods, to use.


title: Partition of sums of squares (ANOVA)

text: Can anyone explain the theory (or the formula) about computing Sum Sq (bold highligh below) related to regression items? The Wikipedia link gives an introduction on how to calculate the total, model, and regression sum of squares. Is it similar to the Sum Sq computation? Is the regression sum of squares equal to (0.000437+ 0.002545+ 0.060984+ 0.062330+ 0.060480)?               TraingData <- data.frame(x1 = c(3.532,2.868,2.868,3.532,2.868,2.536,3.864),                              x2 = c(1.992,1.992,1.328,1.328,1.328,1.66,1.66),                              y  = c(9.040330254,8.900894412,8.701929163,9.057944749,                                     8.701929163,8.74317832,9.10859913)                              )     lm.sol <- lm(y~1+x1+x2+I(x1^2)+I(x2^2)+I(x1*x2), data=TraingData)     anova(lm.sol)          Analysis of Variance Table          Response: y                 Df   **Sum Sq**     Mean       Sq F    value Pr(>F)     x1          1   0.000437  0.000437    0.1055    0.8001     x2          1   0.002545  0.002545    0.6141    0.5768     I(x1^2)     1   0.060984  0.060984   14.7162    0.1623     I(x2^2)     1   0.062330  0.062330   15.0409    0.1607     I(x1 * x2)  1   0.060480  0.060480   14.5945    0.1630     Residuals   1   0.004144  0.004144


title: What type of multivariate linear regression is this?

text: I'm trying to reproduce a result from a book (see bottom) and it doesn't work. I would like to do some further readings about this method but he doesn't specifically give the method other than a formula. I've already tried identifying the method using Wikipedia with no success. This is the linear model: $Z = cX + dY$ He provides an equation for $c$ : $$ c = \frac{{\rm corr}(X,Z) - {\rm corr}(Y,Z){\rm corr}(X,Y)}{1-{\rm corr}(X,Y)^2} $$ The weight $d$ is calculated equivalently. He then writes that $c$ and $d$ can be used to calculate the squared error. Using trial and error, I figured out that the correlation coefficient ${\rm corr}()$ is very likely Spearman's $\rho$ (At least that's the method he used so far to calculate correlation coefficients.). Additionally, he mentions that the means of $X$ and $Y$ are assumed to be vanishing. I'm relatively new to linear regression, so at first I thought it's least squares but this equation doesn't look like it to me. Does anyone recognize this method and can give me a name, so I can read more about it? This all comes from a popular science book on football/soccer statistics. The book is in German. The formula can be found in Appendix A7.3 on pp. 297. The particular example, I'm trying to reproduce can be found on pp. 140.


title: Formula to calculate variance for the sampling distribution of sample statistics (i.e. sample mean, sample median, and sample proportion)

text: What would the formula to calculate variance for the sampling distribution of sample statistics (i.e. sample mean, sample median, and sample proportion) be? At first I thought it would just be the var() command in R with denominator 1/(n-1) but someone pointed out the denominator could be 1/n instead, could someone help explain this one to me? So my main question is what the denominator would be, is it 1/n or 1/(n-1)? Not in particular the entire variance formula.


title: k-Nearest-Neighbor for Factorial Based Data

text: This week I read about using the kNN algorithm to predict the outcome of a continuous variable based on the input of one or more continuous variables. The example was predicting the price of a wine bottle based on rating and age. Is it possible for this to be used to predict a continuous outcome based on factorial based data? For example, predicting the cost of a car based on manufacturer, model, warranty, etc? Or perhaps predicting the attendance of a basketball game based on opponent, day of week, time of day, etc? In similar situations in the past I have converted each factor to a number, but I'm curious to know if there's a way to skip that step.


title: Advice on resampling scheme for a small sample/ costly computation situation?

text: I am testing combinations of preprocessing activities on a small data set (n=48, p=30). The script generates 3200 different versions of the original data and measures how they perform in a classification task against 5 classifiers. Classifier parameters are tuned with 10 times bootstrap. I assume a good sampling scheme would be to repeat 10-fold cross-validation 10 times, but this is computationally too much, at least for my resources:   * 3200 times    * 10 * 10 (cross-validation)   * for each round tune and train the five models 10 * 5  The idea at this moment is to   1. For each data version    2. Tune the models' parameters with 10 times bootstrap sample (2/3 tuning set - 1/3 tuning test set)   3. Train the models with new bootstrap sample given the parameters   4. Test the learned models by taking 3/3 bootstrap repeated 100 times This is a rather complicated sampling scheme and mixes training and test sets. What would be a good yet computationally feasible sampling scheme in this situation?


title: Delta function in monte carlo sampling

text: I am confused by the dirac delta function in the formulation of monte carlo sampling. http://www.cs.ubc.ca/~arnaud/doucet_johansen_tutorialPF.pdf, for instance, defines in section 3.1 page 8 the marginal as $\pi(x_k)=\frac{1}{N}\sum_{i=1}^N\delta_{X^i_k}(x_k)$. What does it mean here? That the density of $x_k$ is the average of the densities of the samples equal to $x_k$? Thanks in advance.


title: Weighted sample and logistic regression

text: I am doing logistic regression in which I have a weighted sample data. Should I should have the weight option on or off when doing logistic regression?


title: Long-run behavior in coin tossing experiment

text: I have read the following information on this site and I am not sure if it is actually true (or if I misinterpret it): > Events that are random are not perfectly predictable, but they have long- > term regularities that we can describe and quantify using probability. In > contrast, haphazard events do not necessarily have long-term regularities. > Moreover, we limit the term probability to events for which we can specify > all possible outcomes. How a fair coin lands when it is tossed vigorously is > a canonical example of a random event. One cannot predict perfectly whether > the coin will land heads or tails; however, in repeated tosses, the fraction > of times the coin lands heads will tend to settle down to a limit of 50%. > The outcome of an individual toss is not perfectly predictable, but the > long-run average behavior is predictable. Thus it is reasonable to consider > the outcome of tossing a fair coin to be random. Let's assume 100,000 coin tossing trials in total and a perfectly fair coin. If we observed that in the first 1000 trials, Heads occurred 600 times and Tails 400 times, then what the above paragraph tells us is (the way I understand this): we could predict that at some point, in the next 99,000 trials (and maybe for a relatively short number of trials), Tails will start coming more often than the Heads so that the final result (of 100,000 tossings) will be approximately 50,000 Heads and 50,000 Tails. Is that interpretation correct or is such predictability rule valid? If yes, is there a way to show this by Bayesian methods? (This is not a student exercise. I am a postdoc in Bioinformatics and, probably, a bit confused.)


title: precision vs covariance

text: I have seen that when modelling the likelihood term with a multivariate Gaussian tend to parameterise the Gaussian with the inverse covariance matrix (precision) rather than the covariance matrix, claiming it might be easier for certain calculations typically done in Bayesian inferencing. Does anyone know what they are alluding to?


title: Unsupervised, supervised and semi-supervised learning

text: In the context of machine learning, what is the difference between unsupervised learning, supervised learning and semi-supervised learning? And what are some of the main algorithmic approaches to look at?


title: Best statistic to compare frequencies in two datasets

text: I am looking for a statistical test to compare frequencies in a particular dataset. In a long string, a particular sequence of characters occurs at a frequency of x per 1000 characters. In some portions of the string, the frequency is higher at y per 1000 characters. Can we tell if y and x are significantly different or not?


title: Evaluating a regression model

text: For classification problems I have been using Neural Networks and measuring Type I and II error using the confusion matrix and its measures as per this resource which is pretty straight forward. When faced with an estimation problem, how would one assess the model performance? Assuming that there are no classes and the output is interpreted in real form. Beyond averaging distance metrics, which does not lend much insight.


title: R package for smooth transition regression models

text: Is there a R package that I can use to specify Smooth Transition Models. I'm looking specifically for something that allows me to specify a TAR model for a given time series. In 2008, a package, RSTAR, was presented in the R User Conference by Mehmet Balcilar. However, the packages doesn't seems to be published and I was wondering if there is any other alternatives in R.


title: How to deal with 0 values on the diagonals for a chi-squared test

text: I am performing a chi-squared test in which differently oriented lines cross- cut each other. Because you cannot have two similarly-oriented lines cross- cutting each other, the diagonals on both my observed and expected tables is 0. How can I then perform a Chi-squared test as to not get a Div/0 problem? Thanks


title: How to estimate the absolute expected difference?

text: Suppose we have two random variables $X$ and $Y$ with unknown distributions. I am looking for an unbiased estimator for the absolute expected difference: $$ | E \\{ X - Y \\} | . $$ For instance, suppose we have unbiased independent observations $x_1, \ldots, x_n$ and $y_1, \ldots, y_m$ ($n,m \geq 1$), how can we use these data to construct an unbiased estimator for the above quantity? Preferably, I would like to find an unbiased estimator for the general case, including perhaps only mild assumptions such as the existence of the first (few) moment(s). However, I am happy with any progress including:   * solutions for specific distributions with unknown parameters (e.g., under the assumption that $X$ and $Y$ are both distributed according to a normal distribution, but the means and variances are unknown),   * efficient biased estimators, for instance that minimize the mean squared error (again, suggestions for the general case and for specific distributions are welcome), and   * any thoughts on good estimators for special cases. Thanks for any input!


title: Finding the covariance matrix to find the best linear predictor (AR(1) model)?

text: I need to find the covariance matrix of two given estimates of an AR(1) model $$X_t = \phi X_{t-1} + Z_t$$ to find its best linear predictor of $X_2$, given $X_1$ and $X_3$. Let W = ($X_1, X_3$)' How could I calculate the covariance matrix of W? Cov(X1, X1), Cov(X1, X3), etc. Thank you


title: What is the relationship between sample size and the influence of prior on posterior?

text: If we have a small sample size, will the prior distribution influence the posterior distribution a lot?


title: Does a within-subjects factor with unequally represented levels mess up a repeated-measures ANOVA?

text: Is there any downside to using a repeated-measure ANOVA when you have a within-subjects 2-level factor where 1 level is represented by more trials than the other. I know in ANOVA we are comparing means, but does the the tighter standard error on the one level throw off the ANOVA? For example, if my IV is stimulus congruency and I present congruent trials 80 times while only presenting incongruent trials 20 times, can I just run the the ANOVA including all of the trials or should I run the analysis by selecting the same amount of trials for each factor level (i.e. select only 20 of the congruent trials to pair with the 20 incongruent trials)? If this is OK for a strictly within-subjects design, does it cause problems if a balanced or unbalanced between-subjects variable is added to the design?


title: Inference in linear model with conditional heteroskedasticity

text: Suppose I observe independent variable vectors $\vec{x}$ and $\vec{z}$ and dependent variable $y$. I would like to fit a model of the form: $$y = \vec{x}^{\top}\vec{\beta_1} + \sigma g\left(\vec{z}^{\top} \vec{\beta_2}\right) \epsilon,$$ where $g$ is some positive-valued twice- differentiable function, $\sigma$ is an unknown scaling parameter, and $\epsilon$ is a zero-mean, unit-variance Gaussian random variable (assumed to be independent of $\vec{x}$ and $\vec{z}$). This is essentially the setup of Koenker's test of heteroskedasticity (at least as far as I understand it). I have $n$ of observations of $\vec{x}, \vec{z}$ and $y$, and I would like to estimate $\vec{\beta_1}$ and $\vec{\beta_2}$. I have a few problems, though:   1. I am not sure how to pose the estimation problem as something like least squares (I assume there is a well-known trick). My first guess would be something like $$min_{\vec{\beta_1}, \vec{\beta_2}} \left(\sum_{i=1}^n \frac{\left(y_i - \vec{x_i}^{\top}\vec{\beta_1}\right)^2}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)\left(\sum_{i=1}^n \frac{1}{g\left(\vec{z_i}^{\top}\vec{\beta_2}\right)^2}\right)^{-1},$$ but I am not sure how to solve that numerically (perhaps an iterative quasi-Newton method might do).   2. Assuming I can pose the problem in a sane way and find some estimates $\hat{\beta}_1, \hat{\beta}_2$, I would like to know the distribution of the estimates so that _e.g._ I can perform hypothesis tests. I would be fine with testing the two coefficient vectors separately, but would prefer some way to test, _e.g._ $H_0: \vec{w_1}^{\top} \vec{\beta_1} + \vec{w_2}^{\top} \vec{\beta_2} \le c$ for given $\vec{w_1}, \vec{w_2}, c$.


title: What is a statistical test of whether a variable is uniformly distributed with kinks?

text: I am trying to come up with a way to statistically test this question, and could use some input, as it's far afield from what I normally do. I have a hierarchical network, like this one: ![enter image description here](http://i.stack.imgur.com/ryKh1.png) And a centrality measurement of each node. What I want to do is test if the measurement for each "layer" (i.e. nodes 1 step from the diamonds, or 2 steps from the diamonds, etc) are uniformly distributed. The problem is that the layers aren't actually independent - the layers closer to the diamond will inherit some of the distortions from the layers below them, so what I actually want to test is that the distribution for layer X is _no more distorted_ than the layer before it. Any ideas how to go about that?


title: Increasing the number of bins in a SAS 9.2 contour plot

text: Borrowing from this post on SAS and R: http://sas- and-r.blogspot.com/2011/07/example-92-transparency-and-bivariate.html I've got a large collection of data that's best approached (or at least well approached) using some sort of scatter plot with smoothing or binning. R's smoothScatter function works marvelously, but it has a large number of bins (128 is the default I believe). The problem is the contour plot in SAS 9.2. There appears to only be 10 levels possible in the SAS contour plot. The data itself has _tons_ of values around zero, and then a spread of more positive values that are all getting lumped together in one category, and thus hidden in the plot. Does anyone know how to increase the number of contours SAS will plot, or is it fixed at 10? For example, the plot in R showing where I want to be roughly, and the closest I have gotten in SAS: ![Smoothed R Plot](http://i.stack.imgur.com/iLeYR.png) ![SAS](http://i.stack.imgur.com/ymLLY.png)


title: Statistically follow up interactions?

text: There is quite some discussion about whether one should statistically follow up interactions. In some fields it is unusual not to do so. I have a three-way interaction with one continuous variable, one dichotomous and one with three categories and asked for advice on this forum. I received usefull advice from Glen de Palma who said I should plot a 2 x 3 interaction once at plus one standard deviation of the continuous variable, and once at minus one standard deviation.I did this. In the paper I give the F-value, df and p-value and than describe without further testing what can be seen in the two plots. Does anybody have a paper to which I can refer that supports this approach? I am curious to hear others opinion on this topic as well. Best wishes, Mariska


title: Can I do a t test if I have little to no variance in one group?

text: I have 4 groups that I am comparing to a criterion. In one of my groups, all the participants answered the same on every item, i.e. there is no variance. How do I deal with that in my ANOVA? Also, what do I make of that in a t test I am running comparing it to a criterion since I will get no error term? If I do include one participant that I am not sure I am including in my student, the variance is not totally uniform with 1 different observation out of 37 but when I run it, it is not significant because the variance is too small. I understand that there is nothing I can do computation-wise. I am asking how one handles it conceptually.


title: Mathematician wants the equivalent knowledge to a quality stats degree

text: I know people love to close duplicates so I am **not** asking for a reference to _start_ learning statistics (as here). I have a doctorate in mathematics but never learned statistics. What is the shortest route to the equivalent knowledge to a top notch BS statistics degree and how do I measure when I have achieved that. If a list of books would suffice (assuming I do the exercises lets say), that's terrific. Yes, I expect working out problems to be an implicit part of learning it but I want to fast track as much as realistically possible. I am not looking for an insanely rigorous treatment unless that is part of what statistical majors generally learn.


title: Question about independence assumption for ANOVA, t-test, and non-parametric tests

text: I'm a novice in statistics and I have some confusion about the assumption of independence for statistical tests.   1. I searched the Internet and some information says that for the t-test, the observations in the two groups should be independent (that is, measurements in sample 1 and measurements in sample 2 should be different). Some other information says that all observations (even in the same group) should be independent. Which one is correct?    2. Is the independence assumption for ANOVA and the independence assumption for the t-test the same?    3. Do non-parametric tests, such as the Wilcoxon signed rank test, also need to satisfy the independence assumption?


title: Bias of more than one endogenous variables

text: I may have a model with omitted variables that are correlated with my predictor variables. If I have in my model, let's say, two endogeneous variables X1, X2, but I am interested in obtaining an unbiased coefficient of just one variable X1, do I have to find instrument for both variables X1, X2? Or is it sufficient to instrument just for X1 and leave X2 alone? X2 is just a minor control variable that I don't have an instrument for. However, my main interest is the coefficient on X1, for which I have a good instrument for. Is this okay to proceed? Many thanks!


title: Multistage sampling for Population Median

text: Has there been any research done on how to estimate the population median using multistage sampling? Simple random sampling is not possible in my case I would really like to use something like multistage sampling.


title: Interaction between dummies and a continuous variable

text: Lately I have had some trouble interpreting the interaction term I use in a regression. I have two dummy variables and a continuous variable. For the explanation, BLACK takes the value 1 if the respondent is black and 0 otherwise. MAR takes the value 1 if married and 0 otherwise. I also include a continuous variable for the education of the respondent EDUC. I create the variable MARBLACK which takes the value one for married blacks. When building my model, I estimate: $$ EARNINGS= \alpha + \beta_1 BLACK + \beta_2 MAR + \beta_3 MARBLACK + \beta_4 BLACK*EDUC + \beta_5 MAR*EDUC + \beta_6 MARBLACK*EDUC + \beta_7 EDUC $$ Here is my vision of it: $\alpha + \beta_1 BLACK + \beta_2 MAR + \beta_3 MARBLACK + \beta_4 BLACK*EDUC + \beta_5 MAR*EDUC + \beta_6 MARBLACK*EDUC + \beta_7 EDUC $ gives the predicted value for married blacks $\alpha + \beta_1 BLACK + \beta_4 BLACK*EDUC + \beta_7 EDUC $ gives the predicted value for non-married blacks $\alpha + \beta_2 MAR + \beta_5 MAR*EDUC + \beta_7 EDUC $ gives the predicted value for married non-blacks $\alpha + \beta_7 EDUC $ gives the predicted value for non-married non-blacks (my referent class). So, if I want to calculate the _premium_ of being married and black, I would keep $\beta_1 BLACK + \beta_2 MAR + \beta_3 MARBLACK + \beta_4 BLACK*EDUC + \beta_5 MAR*EDUC + \beta_6 MARBLACK*EDUC $ and not simply $ \beta_3 MARBLACK + \beta_6 MARBLACK*EDUC $ My intuition is that being black and married, three effects accrue to the respondent: the fact of being black, the fact of being married, and the fact of belonging to both categories. However, discussing with another researcher, she claims that it is only the MARBLACK coefficients that I should consider. She explains it taking conditional expectations: Let $E[EARNINGS|BLACK=1,MAR=1]= \alpha + \beta_1 BLACK + \beta_2 MAR + \beta_3 MARBLACK + \beta_4 BLACK*EDUC + \beta_5 MAR*EDUC + \beta_6 MARBLACK*EDUC + \beta_7 EDUC $ and $E[EARNINGS|BLACK=0,MAR=1]= \alpha + \beta_2 MAR + \beta_5 MAR*EDUC + \beta_6 + \beta_7 EDUC $ She then compute the effect of the change, i.e. $E[EARNINGS|BLACK=1,MAR=1]-E[EARNINGS|BLACK=0,MAR=1]= \beta_1 BLACK + \beta_3 MARBLACK + \beta_4 BLACK*EDUC + \beta_6 MARBLACK*EDUC $. She calls it EARNINGSbis. She then goes on and do the following: $E[EARNINGSbis|MAR=1]=\beta_1 BLACK + \beta_3 MARBLACK + \beta_4 BLACK*EDUC + \beta_6 MARBLACK*EDUC$ and $E[EARNINGSbis|MAR=0]=\beta_1 BLACK + \beta_4 BLACK*EDUC$ She again computes the effect of the change and she is left with $\beta_3 MARBLACK + \beta_6 MARBLACK*EDUC$ Could any one of you please clarify what my method and her method mean ? And especially where they deviate in terms of interpretation ? Thank you very much. Kind regards, Olivier.


title: How to draw random samples from a non-parametric estimated distribution?

text: I have a sample of 100 points which are continuous and one-dimensional. I estimated its non-parametric density using kernel methods. How can I draw random samples from this estimated distribution?


title: Visualize process data

text: I have a dataset of order processing with 8 million rows with following columns:   * `HistoryId` \- Identity column of the records   * `ItemId` \- Id of the item tracked   * `Previous status` \- status of Item before the history record created   * `New status` \- status of Item when history record created   * `Time till status change` \- The time difference between the previous record of the item till this history record created. I want to visualized the these data with a graph / chart that displays the flow similar to this:               Status 1 --> Status 2 --> Status 3 --> Status 4 --> Status 5                                            Status 5                               Status 4 --> Status 5     Status 2 --> Status 3 --> Status 5 --> Status 2 --> Status 4 --> Status 5                                                         Status 5                           --> Status 4 --> Status 5                  Status 4 --> Status 5      In example above - assume that all of the items have:   1. First status is either `Status 1` or `Status 2`.    2. All item start with `Status 1` have 2nd status is `Status 2`, while for item start with `Status 2` partially `Status 3`, and others are `Status 4`   3. Moving on similar logic of (2.) with break down by each flow. I want to explore the data and visualize all possible flows of the data. Then summarize them by counting the number of items go through each `Status` of each flow. Currently, I am doing them manually. Is it possible to automate this and visualize it on graph? Sample Data:               HistoryID     ItemId     Previous_status     New_Status        1             1          NA                  status_1          2             2          NA                  status_2        3             1          status_1            status_2        4             1          status_3            status_4       5             2          status_1            status_3        6             1          status_4            status_5         7             2          status_3            status_5          And here is a sample output even this output is not completely what I wanted - Explanation:   * `X` \- is the index of status that an Item have in it life cycle.   * `Y` \- is the status name   * `Size` \- the number of item go through status Y at index X.   * Skip the color as it was exclude from the data in example.   * As you can see most of item have the status Y as the beginning status (the fist left column).   * Then moving to the right they are breaking down to other status (sometime they may come back the beginning status but in later index)   * The limitation of this is not show the flow by detail even I can see the flow by the changing in size of the point overtime.   * What I really want is something like a decision tree where you can see how item flow through the status, and each flow is separate from each others. ![Sample out put but not exactly](http://i.stack.imgur.com/811af.png) * * * I managed to used the riverplot package. However I have some prolems as image below. Anyone know how to display the lable at the side instead of right on the spots? as I have a about 30 status and it is very confusing when having them display like the image below. For your reference, here is the link to the package riverplot tutorial ![enter image description here](http://i.stack.imgur.com/uksfg.png)


How to calculate in-control ARL for control charts? For example, to design an EWMA Control Chart, how we could conclude that in-control ARL is approximately 500 for Lambda=0.1 and L=2.7? Thank you!


title: Drawing multiple barplots on a graph in R

text: I would like to plot four barplots on a single graph in R. I have used the following code. Here, how can keep a legend on top of the graph, specifically the legend should be between 2 and 3 barplots. I also tried with `par(mar=c(4.1,4.1,8.1,4.1)` but there is no success. Moreover, I also tried to run `legend()` after the second barplot, but there is no use. The legend is for all the four barplots. Please help me in this.                   par(mfrow=c(1,4))         barplot(t(A), beside=T, ylim=c(-100,100),..)         barplot(t(B), beside=T, ylim=c(-100,100),..)         barplot(t(C), beside=T, ylim=c(-100,100),..)         barplot(t(D), beside=T, ylim=c(-100,100),..)         legend(...)


title: How can this inequality be otherwise expressed?

text: Consider an ordered pair of random variables $(x_i,y_i)$, and the following relationship: $$\frac{\left(\sum_i^n x_i\right)^2}{\sum_i^n x_i y_i}\leq \sum_i^n\frac{x_i}{y_i}$$ Is there some useful way of restating the relationship in terms of statistical functions? Can the relationship between the two variables be summarized in some intuitive, informal way? The terms on LHS appear in deriving the least squares fit equation, but I don't recognize the RHS term.


title: How to analyze and interpret A/B test results via Bootstrap method?

text: First of all, apologies if this is a stupid question, I'm fairly new to statistics, so still finding my way a bit in terms of the fundamentals! We've run a split test of a new product feature and want to measure if the uplift on revenue is significant. Our observations are definitely not normally distributed (most of our users don't spend, and within those that do, it is heavily skewed towards lots of small spenders and a few very big spenders), so we've decided on using bootstrapping to compare the means, to get round the issue of the data not being normally distributed. So my results show that we do have an uplift of around 8% vs. control. I now want to calculate how confident I can be in this uplift. Is it as simple as measuring the proportion of the probability density function below zero, for the PDF that is test group PDF minus control PDF? (e.g., that portion reflects the % chance that my 2 PDFs are not different?) Any help would be much appreciated.


title: Is bootstrapping the right method for extreme distributions?

text: Wondering if anyone has an opinion on whether bootstrapping the difference in means is the right method given I have a situation with extreme data points. I've decided to use this as I don't think a t test is appropriate I have about 30k observations per group (3 groups) My situation is about spend, and I have extreme outliers: the outliers aren't quite like an "income" distribution. That is, most users (95%+) will spend zero, a subset of users will spend 5 - 10 dollars. some will spend about 20 or 50 dollars and then a select few will spend 500+, with a couple of users spending 5000 or 10000+ I am trying to test which group brought in the most revenue per user. Can anyone offer any advice on which statistical test is best suited?


title: Longitudinal predictive models

text: I have a predictive model to construct and am looking for some ideas on the approach. I have a large training dataset of customers' balances (think a savings account) at time t=1 to 36 representing monthly averages over three years. I will want to construct a model to predict the value of a customer's account at t=4..36 using only information available at the end of t=3. _Hence, this is a forecasting type problem, but not in the time series sense where you observe a sequence of data significantly longer than the time horizon you are trying to predict._ I am thinking of approaching this like a discrete hazard model, where a row in the dataset is created for every customer/month combination (so here each customer has 33 rows in the data and a variable called t takes on the values 4 to 36). Repeated for every month is the "initial state" variables for the customer known at the end of month 3. This is the setup of the so-called person period data. I would then learn a regression model on this data using the initial state data and 't'. There would be no accounting for the fact that the customers are repeated in the data. This works for a discrete time hazard model, but my question is if this is OK for a linear regression or a machine learning algorithm (e.g. neural net)? Is there a better way? **ADD** : Specifically I am wondering if such a model will fail if it does not explicitly account for the repeated measures and the correlation between rows (same customers at various values of 't'). When dealing with a predictive model, how could this be accounted for?


title: P value for bootstrapping two means in laymen's terms

text: I've bootstrapped two means - each sample had 10,000 observations so I resampled each mean with replacement of 10,000 numbers from each sample. I did this over 50,000 iterations I've plotted the distributions of the bootstrapped means and it looks like they are different. However I'm lost in the next step. Is someone able to tell me in laymen's terms how I can calculate a P value and how to know if one bootstrapped distribution is significantly different from the other? I am using R to do my bootstrapping.


title: Pooling regression results in SPSS

text: I have to solve the following issue:   1. I run my linear regression model many times (let's say 1000 times) with two variables: y - continuous dependent variable, x - continuous independent variable (mean of several consequent measurements).   2. The independent variable in each model was randomly drawn using its mean and standard deviation   3. I have the regression coefficient and standard error for this independent variable in each of the models. Somehow I have to combine these results into one regression result. As far as I know the regression coefficients of 1000 models can be just averaged. However, this is not really clear to me how can I estimate the total variance of 1000 models.


title: How to find Variance from rpart object in R

text: I am fairly new to R and am creating an rpart object using the method "anova" (Since I am dealing with continuous columns)               EstimationModel <- rpart(OutpoutColumn ~ ., data=DATAFRAME,method="anova")      Now I want to find out the variance of each node. The rpart object gives me the mean but how do I get the variance?


title: Normality test of latent residuals (Heckman) Probit?

text: I am running a Heckman Probit. Both steps assume normal errors of the latent variable, correct? My data is from a survey and a lot of variables are strongly skewed, so I am worried whether this assumption holds. Is there any test I can use? The sample is 10 000, does it help? Also, can I make the errors closer to normal by transforming the data in some way? Say, taking logs? Thanks!


title: power/sample size for a complex set of equations

text: I run an experiment in which several teams perform several tasks each (not all teams perform all tasks, but in general, I can make sure that enough teams perform any given task, and that each team performs enough tasks). I have 30 different tasks, but for the sake of the question let's say its t tasks, g teams and that each team performs k tasks. The outcome depends on some team level known variables (controls), but mainly on 9 variables that I need to estimate. in each task the weight of the different variable is different (but known). This can be summarized in the following equation: $T_i = controls_{team(i)} + \sum_{j=1}^9 w_j^{task(i)}*v_j + \epsilon_i$ my question would be how can I compute power in this setting. specifically I want to be able to calculate the sample size (e.g. the number of teams and/or the number of repetition per team and/or the number of different tasks) in order to estimate the $v_j$ variables reliably. I know that if I had only one task this would look like a regular linear regression, but as the weights are fixed, it would be non-identified model. with 30 tasks it is identified model (at least in theory), but I can't figure out how to calculate power. It seems related to SEM in that each tasks yields a different equation based on its known weights. it also seems relevant to panel data analysis because of the team*task structure. I'm not expert in either, so I may be wrong. Any help/insights/references are welcome.


title: Hypothesis testing - Wilcoxon test, bootstrapping, or something else?

text: A colleague has developed a treatment for to "prevent falls" in cognitively impaired, psychiatric patients. Since this would be very useful treatment in this population, we especially do not want to make a Type II error (i.e., fail to reject the null, when we should reject it). Since the data is not normally distributed, another colleague evaluated the full data using (appropriately, I believe) the Wilcoxon test, and did not find significance. There may be valid methodological reasons for this, which I may follow-up on later with another question. I was concerned about committing a Type II error, and obtained some preliminary data, which I have below. These data reflect pre/post scores (# of "falls") for the same patients (no control group), so should be considered "paired" and not independent:               pre <- c(9,8,37,12,8,3,4,4,3,5,4,8,4,8,9,11,2,4,0,0,5,12,10,2,8,3,0,22,1,0,0,5,0,3,1,5)          post <- c(10,8,6,4,5,2,4,4,2,2,1,7,2,1,3,9,2,2,0,0,6,16,4,3,4,7,0,10,3,0,0,4,0,1,1,5)      When I ran a bootstrapping procedure on this preliminary data (adapted from Crawley, The R Book, p. 385)               preBoot <- numeric(10000)     for (i in 1:10000) {preBoot[i] <- mean(sample(pre, replace=T)) }     quantile(preBoot, c(0.025, 0.975))      and compared the post mean to the bootstrap estimate of the sampling distribution of the mean, I found that the treatment did have a significant beneficial effect. To evaluate significance, I simply took the quantiles for the sample estimate at 0.025 and 0.975; is this correct or am I confusing what I would do with a normal distribution with the distribution of the sample estimates of the mean? Also, using wilcox.test in R on the preliminary data (i.e.,               wilcox.test(pre, post, paired=T, exact=F)      shows this to be significant. I would like to know, before I go further, did I use the bootstrapping procedure correctly, and is this a legitimate test for these type of data? Are there other tests we should consider, and what would be the best way to report this? I am especially interested in methods that would allow us to obtain confidence intervals. Also, I see in this previous question Wilcoxon one-tailed test the response indicated that "keep in mind that it's generally not advisable to use one- tailed tests," but if I'm interested specifically in fewer falls after the treatment intervention, wouldn't a one-tailed test be appropriate? **Additional Info Update:** I just found a wonderful review on count data and analysis by Neal Alexander, _Review: analysis of parasite and other skewed counts_ , accessed via PubMed http://www.ncbi.nlm.nih.gov/pubmed/22943299 which discusses the issues I've been facing in a very accessible manner (and it is free online). Other reading this question may also find this quite helpful. I'm still digesting this info. This probably belongs in a new question, but in essence, I believe in my field (clinical psychology) the standard way of looking at this data would be with a Wilcoxon test, with perhaps a square root transformation and a t-test running in second place. Most people currently don't use R, and so don't seem to be aware of or use bootstrapping, which I actually believe would be better than the above two methods. If anybody has further info/or information to the contrary, I would appreciate it).


title: Repeated measures vs. cluster in logistic regression

text: I have data from an experiment in which participants played a repeated game 50 times. That is, I have 50 data points per participant. I want to run a logistic regression and I have of course to take into account that I have data points that are correlated, but I don't know **if I should use repeated measures or if I should work with clustered sample**. I have found some very clear articles on why clusters are better than other techniques (Williams, 2000; Petersen, 2009), but I have found nothing on repeated measures or on the comparison between the two. I would appreciate it if you could help me on this one! Thanks!


title: Best method to estimate the mean of a normal distribution?

text: Let $X = ( x_1, ..., x_n ) $ be $n$ samples from a normal distribution with unknown mean. What is the best estimator for this mean? I can think of at least 2 unbiased estimators:   * The empirical mean $ \hat \mu_1 = \frac{\sum_i x_i }{n} $    * A bayesian approach $ \hat \mu_2 = E[ P( \theta | X ) ] $, where $ P( \theta | X ) $ is the posterior distribution over the mean. The prior distribution over the mean is $ \cal N ( m, s^2) $ for some meta parameters $(m, s^2)$. Are there other biased or unbiased estimators for this problem? Which one is the best?


title: Building prediction model with estimated predictor variables

text: I'm planning to use logistic regression with multiple (~5) predictor variables to predict whether something happens or not. I have two types of predictor variables: known (measurable) variables and unknown (not measurable) variables. From the measurement perspective only the former can be measured and the later have to be estimated. To rate how good my estimations are, I created a special setting (with the help of another perspective) so that for every estimated value I have also the "real" value. If I had only measured variables, I would divide my data into a training set and a test set in order to build a regression model and to obtain its performance characteristics. However, now I have 2 training sets and 2 test sets:   a) training set with real values   b) training set with estimated values   c) test set with real values   d) test set with estimated values Since my special setting is not suitable for daily practice, the regression model I build will most likely operate on estimated data only. Which combination of a) or b) and c) or d) should I use and why?


title: Taylor expansion to contain sample mean, sample variance, sample skewness, and sample kurtosis

text: I have the following expression: $$\frac{1}{p} \ln\left(1+\frac{p^1}{1!n}\sum_{i=1}^n x_i + \frac{p^2}{2!n} \sum_{i=1}^n x_i^2 + \frac{p^3}{3!n} \sum_{i=1}^n x_i^3 + \frac{p^4}{4!n} \sum_{i=1}^n x_i^4 + \cdots \right)$$ Now let $$Y = \frac{p^1}{1!n}\sum_{i=1}^n x_i + \frac{p^2}{2!n} \sum_{i=1}^n x_i^2 + \frac{p^3}{3!n} \sum_{i=1}^n x_i^3 + \frac{p^4}{4!n} \sum_{i=1}^n x_i^4 + \cdots$$ Then we have $$\frac{1}{p}\ln\left(1+Y\right)$$ Using the Taylor series expansion on log, we have $$\frac{1}{p}\ln(1+Y) = \frac{1}{p}\sum_{n=1}^{\infty} (-1)^{n+1} \frac{Y^n}{n} = \frac{1}{p}\left[Y - \frac{Y^2}{2} + \frac{Y^3}{3} - \frac{Y^4}{4} + \frac{Y^5}{5} - \cdots\right] $$ My question is, how can I expand the above expression so that my final approximation contains the following variables: Sample mean: $\displaystyle \overline{x} = \frac{1}{n} \sum_{i=1}^n x_i$ Sample variance: $\displaystyle s^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i - \overline{x})^2 $ Sample skewness: $\displaystyle g_1 = \frac{\sum_{i=1}^n (x_i - \overline{x})^3}{(n-1)s^3}$ Sample kurtosis: $\displaystyle g_2 = \frac{\sum_{i=1}^n (x_i - \overline{x})^4}{(n-1)s^4}$ * * * To illustrate exactly what I mean, I know how to obtain the approximation so that it includes the sample mean and sample variance. Start with the expression as derived above: $$\frac{1}{p}\left[Y - \frac{Y^2}{2} + \frac{Y^3}{3} - \frac{Y^4}{4} + \frac{Y^5}{5} - \cdots\right]\ \ \ \cdots \ \ (1)$$ Ignore the terms from $\frac{Y^3}{3}$ onwards and substitute in the original expression for $Y$: $$\frac{1}{p}\left[ \left(\frac{p}{n}\sum_{i=1}^n x_i + \frac{p^2}{2!n} \sum_{i=1}^n x_i^2 + \cdots \right) - \frac{1}{2}\left(\frac{p}{n}\sum_{i=1}^n x_i + \frac{p^2}{2!n} \sum_{i=1}^n x_i^2 + \cdots\right)^2 + \cdots \right]$$ $$=\frac{1}{p} \left[\frac{p}{n} \sum_{i=1}^n x_i + \frac{p^2}{2n} \sum_{i=1}^n x_i^2 - \frac{p^2}{2n^2}\left(\sum_{i=1}^n x_i \right)^2 + \cdots\right] $$ $$\approx \frac{1}{n} \sum_{i=1}^n x_i + \frac{p}{2n} \sum_{i=1}^n x_i^2 - \frac{p}{2} \left( \frac{1}{n} \sum_{i=1}^n x_i\right)^2 $$ $$ = \overline{x} + \frac{p}{2} \left[\frac{1}{n}\sum_{i=1}^n x_i^2 - \left(\frac{1}{n}\sum_{i=1}^n x_i \right)^2\right]$$ $$= \overline{x} + \frac{p}{2} \left[\frac{n-1}{n} s^2 \right] \ \ \ \cdots \ \ (2)$$ * * * As can be seen above, the final approximation contains the sample mean and sample variance. However, I am not sure how exactly I can manipulate the expansion to contain sample skewness and sample kurtosis. Any help would be appreciated.


title: How do you calculate the likelihood ratio test statistic in this problem?

text: I am trying to solve this problem but I am not sure how to proceed to get the formula for the test statistic and critical values. Suppose $X_1,X_2,\ldots X_n$ are i.i.d. observations from a multivariate normal distribution $N(\mu,\Sigma)$ where $\Sigma$ is known. Further assume that **R** is a given matrix and **r** a given vector. Use the likelihood ratio procedure to produce a test statistic for $$H_0\colon R\mu = r \\\ H_1\colon R\mu \neq r$$ Give explicit formula for the test statistic and the critical values. Can someone give me a hint how to proceed?


title: Evaluation of groups of parameters across multiple data sets in logistic regression

text: I have a logistic regression problem where I want to evaluate the aggregate contribution of groups of parameters (independent variables) across different sets of observations. For instance, let's say I have 10 total parameters, 5 of which belong to Group1, for a given data set. I can easily compare the full model to the model without Group1, using a pseudo R^2, likelihood ratio test, AIC, or similar, to describe the contribution of Group1 to DataSet1. What I need to do is aggregate these comparisons across different data sets, where the solutions are not necessarily related (ie the 'true' coefficients are not the same per parameter across data sets) - I am interested in the amount of information conveyed by the group of parameters in aggregate, not the coefficients themselves. So, I cannot train a single model using all the observations. One example solution I came up with is to add the model deviances (deviance = -2*loglikelihood) after solving each model independently and calculate the likelihood ratio/AIC/pseudo R^2 etc as if it were one large model. One complication of this, however, is that the number of parameters in the group or full model might change in each data set. For instance, in DataSet1, Group1 may have 5 parameters (or 6 parameters but no observations for parameter 6 - the params are usually categorical - etc), but in DataSet2, Group1 may have 3 or 7 or 150 parameters, etc. So, to calculate AIC or p-value for the likelihood ratio test, I have calculated the effective #params as the weighted average #params based on #observations in each data set. Alternatively, I could fix the #params as the max from each data set, but I fear this would artificially squash my results. I'd love some feedback on this approach or pointers to alternatives.


title: Is it necessary to detrend and decycle time-series data when using machine learning methods?

text: For example: I want to forecast future values of a time-series based on previous values of multiple time-series' using a ANN and/or SVM. Inputs will be lagged values from each time series, and the outputs will be one-step-ahead forecasts (forecasts with further horizons will be done by "rolling" the predictions forward using previous predictions). Now, shouldn't SVMs and ANNs be able to learn trends and cycles? Would they not be able to learn things like, "with all else being equal the output of this series should be 2x the previous output?" Or, if I provide a categorical variable for month, "since it's January, divide the prediction I would've made by 2?" Would attempting to decycle and detrend the data result in imposing more bias than necessary?


title: How to determine a user's favorite content producer from individual ratings?

text: Consider the following scenario: Alice subscribes to a video rental service that allows her to watch movies. Every time A watches a movie, she rates it either thumbs up (1) or thumbs down (0), and she then chooses the next movie she wants to watch. Every movie belongs to exactly one director, and a director can have directed many movies. The question is, what is the best way to determine who A's "favorite" director is? My initial thought was do something like:   * For each director of at least one movie that A watched, calculate the lower bound of some binomial confidence interval (e.g. Wilson score interval) as A's "favorability" score for that director However, this binomial approach seems flawed because it ignores a seemingly crucial piece of information: Alice has an entire universe of movies to choose from, and if she consistently chooses to watch movies from a certain director, then doesn't that tell us something about her preference for that director, even if she then rates that director's movies below her average? I feel like there must be some element of "voting with one's feet" that is ignored if we only consider ratings on the movies that were watched. What is the best way to combine both the selection of movie/director with the ratings on individual movies to determine who is A's favorite director? It seems like A's preference for director D has to be a function of A's ratings on D's movies that she watched, and also the % of all of D's movies that A chose to watch. UPDATE: I should make clear, the problem I'm dealing with isn't quite as simple as the thumbs up/thumbs down case, it's really more like "A watches a movie and then checks a box if she liked it." So each viewing does result in a 0 or a 1, but the absence of checking a box isn't quite the same thing as a "thumbs down" because the viewer might only feel compelled to check the "approve" box if she _really_ likes something. All the more reason that the choice of what to watch has to factor into preferences


title: Confidence interval: "$v_2$ is at most $x$ larger than $v_1$, with a confidence of $95\%$"

text: I have two unknown values $v_1$ and $v_2$ and two sets of measurements $S_1$ and $S_2$ of these values. The measurements can be assumed to have a normal distribution. I would like to say something like _"$v_2$ is at most $x$ greater than $v_1$ with a confidence of 95%"_. So, how small may $x$ be, in order to be 95% sure that that $v_2$ is greater than $v_1 + x$?


title: Mixture distributions moments if one distribution has undefined/infinite moments

text: Consider probability density functions $f_{1}\left(x\right)$ and $f_{2}\left(x\right)$ and the mixture distribution $$f_{3}\left(x\right)\equiv pf_{1}\left(x\right)+\left(1-p\right)f_{2}\left(x\right)$$ Assume that the moments associated with $f_{1}\left(x\right)$ are all clearly defined, but the moments of $f_{2}\left(x\right)$ may not necessarily be (e.g. a Cauchy or Levy distribution). To calculate the expected value one would find $$\int_{-\infty}^{\infty}xf_{3}\left(x\right)dx=p\int_{-\infty}^{\infty}xf_{1}\left(x\right)dx+\left(1-p\right)\int_{-\infty}^{\infty}xf_{2}\left(x\right)dx$$ Presumably if $\int_{-\infty}^{\infty}xf_{2}\left(x\right)dx$ is infinite or undefined, then the expected value only exists when $p=1$. Is this logic correct, even if $p$ is close to but less than $1$? Does it extend to higher moments as well?


title: Calculate probability (area) under the overlapping area of two normal distributions

text: I have two normal distributions defined by their averages and standard deviations. Sample 1: Mean=5.28; SD=0.91 Sample 2: Mean=8.45; SD=1.36 You can see how they look like in the next image: ![enter image description here](http://i.stack.imgur.com/ZfxV9.png) How can I get the probability to obtain an individual from the overlapping area (green)? Is the probability the same as the area?


title: Using MICE in R: is it possible to impute only sub-sections of the data?

text: When using the `mice` library in `R` to impute data I encounter the following problem. I have a data matrix with missing information on `y1` and `y2` and predictor variable `x`. However, instead of imputing all missing data `NA` on `y1` and `y2`, I would like to impute only the subset of data identified by `s1==1` and `s2==1`. Due to the data structure I cannot just drop cases with `s1==0` or `s2==0`. The data look similar to this:               n=20     y1=c(rnorm(n))     y2=c(rnorm(n))     x=c(y1+y2+rnorm(1))     y1[c(1:10,16:20)]<-NA     y2[11:20]<-NA     s1=c(rep(0,10),rep(NA,5),rep(1,5))     s2=c(rep(NA,10),rep(1,5),rep(0,5))     (data<-data.frame(y1,y2,x,s1,s2))      The question is now to have `mice` impute the missing data on `y1` and `y2` for `s1==1` and `s2==1`, but not for `s1==0` and `s2==0`. It is obvious that selecting on `data$s1==1` will drop all observations on `y2`. The reason I want to do this, is that in a larger data file the number of necessary imputations leads to covergence problems in `mice` with `meth='polyreg'`. By setting the `s`indicators I am trying to reduce the complexity to the parts of data that need imputation for relevant inference.


title: How to separate groups using PCA?

text: I have two groups: `A`) controls: $25\times2000$, `B`) patients: $12\times2000$. Group `A` has one PC explaining a lot of the variance. (Suggesting a homogenous population?) Group `B` cannot be explained by the first 2 or 3 (each of the PCs seem to be dominated by single patients). If I use the group `A` PC1 as an example of "normal control", is there a way to show how "far" or "close" each of the patients are from this? I.e. How can I establish that some patients are very similar to the controls and some aren't? * * * I have 25 controls but only 12 patients. Each of them have a spectrogram (time x freq x channels) that is collapsed down to 1 vector.I have tried PCA on just the controls and get a good 1st PC (it's not dominated by any one subject) and looks like the individual data. With patients alone, I do not get this. If I do the PCA on both patients and normals [controls] together, the 1st PC still looks like the normals but it doesn't separate the 2 groups. So **it was suggested to me to take the PCA of just normals and then project the patients into this space but how to do that is unclear.** Just want to add - I don't think there are two populations; i.e., controls and patients are not very different (at least on the measure that I am using here). But, some patients are "farther apart" than others from the controls; i.e., some look like normals but others are very far from them.


title: Validation - correctly compare and validated imputation models

text: I've seen a lot of interesting questions here about multiple imputation and also great answers that helped me a lot to impute my data. I've used Predictive Mean Matching, EMB and I would like to use Random Forest but I'm not being able to run the model (so sad). I also have results from the complete dataset. Overall, the dataset is composed by 8364 observations and only 2860 are complete observations without any missing data. The imputed models are somewhat similar results, the imputed models have narrow confidence intervals and the model obtained by PMM through Mice does not select one variable due to the use of predictors _(Explanation: This function creates a predictor matrix using the variable selection procedure described in Van Buuren et al.~(1999, p.~687--688). The function is designed to aid in setting up a good imputation model for data with many variables.)_ I know I can validate the model with the function validate from the package rms (I plan to do that but I need to understand it better). But who can I compare imputed models and draw some conclusions. Is the output of validate() enough? For instance, if I wanted to repeat another analysis with the same data, could I chose the "best" method or it's better to repeat the analyses for all methods? I understand this question can be quite basic but I'm having some issues figuring out the best way to proceed. Thank you very much!


title: Failed convergence in neuralnet() function

text: I'm having a problem with the neuralnet() function due to lack of convergence. The example code below illustrates the problem. I get the same kind of result using my own real data as well. nnet() is able to find an answer with the same data (see the nnet example in the link). The defaults in the documentation appear reasonable so I left them alone. Does this seem like it's 'just the data' or could the method be tuned by setting parameters? It's not obvious which ones. The example is from http://beckmw.wordpress.com/2013/11/14/visualizing-neural- networks-in-r-update/               > library(clusterGeneration)     > seed.val<-2     > set.seed(seed.val)     > num.vars<-8     > num.obs<-1000          > #input variables     > cov.mat<-genPositiveDefMat(num.vars,covMethod=c("unifcorrmat"))$Sigma     > rand.vars<-mvrnorm(num.obs,rep(0,num.vars),Sigma=cov.mat)          > #output variables     > parms<-runif(num.vars,-10,10)     > y1<-rand.vars %*% matrix(parms) + rnorm(num.obs,sd=20)     > parms2<-runif(num.vars,-10,10)     > y2<-rand.vars %*% matrix(parms2) + rnorm(num.obs,sd=20)          > #final datasets     > rand.vars<-data.frame(rand.vars)     > resp<-data.frame(y1,y2)     > names(resp)<-c('Y1','Y2')     > dat.in<-data.frame(resp,rand.vars)          > #neuralnet function from neuralnet package, notice use of only one response     > library(neuralnet)     > form.in<-as.formula('Y1~X1+X2+X3+X4+X5+X6+X7+X8')     > set.seed(seed.val)     > mod2<-neuralnet(form.in,data=dat.in,hidden=10)      > Warning message: algorithm did not converge in 1 of 1 repetition(s) within > the stepmax stepmax is at the default of 1e5


title: When is Likelihood Function Positive Semidefinite

text: This may be a very misinformed question, but I cant figure out why its not true. Here goes: According to Wikipedia and this post, the hessian of a likelihood function equals the information matrix, or the covariance matrix of the score functions, ie: $$I(\theta)_{i,j} = \mathrm{E}_\theta \left[ \left(\partial_i \log f_{X\mid\Theta}(X\mid\theta)\right) \left(\partial_j \log f_{X\mid\Theta}(X\mid\theta)\right)\right] \, , = -E\left[\frac{\partial^{2} \log(f(X|\theta))}{\partial \theta_{i} \partial \theta_{j}}\bigg|\theta\right]$$ If this is true, wouldnt these conclusions hold true:   1. The Hessian of the likelihood functions is always positive semidefinite (PSD)   2. The likelihood function is thus always convex (since the 2nd derivative is PSD)   3. The likelihood function will have no local minima, only global minima!!! These results seem too good to be true, but I cant seem to understand why they are false. Thanks!


title: Confidence interval for cross-validated classification accuracy

text: I'm working on a classification problem that computes a similarity metric between two input x-ray images. If the images are of the same person (label of 'right'), a higher metric will be calculated; input images of two different people (label of 'wrong') will result in a lower metric. I used a stratified 10-fold cross-validation to calculate the misclassification probability. My current sample size is around 40 right matches and 80 wrong matches, where each datapoint is the calculated metric. I'm getting a misclassification probability of 0.00, but I need some sort of confidence interval / error analysis on this. I was looking into using a binomial proportion confidence interval (where I'd use the results of the cross-validation as a correct labeling or incorrect labeling for my number of successes). However, one of the assumptions behind the binomial analysis is the same probability of success for each trial, and I'm not sure if the method behind the classification of 'right' or 'wrong' in the cross-validation can be considered to have the same probability of success. The only other analysis I can think of is to repeat the cross-validation X times and calculate the mean/standard deviation of the classification error, but I'm not sure if this is even appropriate as I'd be reusing the data from my relatively small sample size several times. Any thoughts? I'm using MATLAB for all my analysis, and I do have the Statistics toolbox. Would appreciate any and all assistance!