title: Computing conditional expectation of ordered normal random variables

text: There are $m$ normally distributed, independent random variables $N_1, \ldots, N_m$ with distinct means $\mu_1, \ldots \mu_m$ and standard deviations $\sigma_1, \ldots, \sigma_m$. Then, we observe a permutation of the numbers $\pi = \{1, \ldots, m\}$. How can we efficiently compute the conditional expectation of the random variables in same ordering as this permutation? _Added bonus: how can we compute the conditional variance?_ An example:   1. we have four independent random variables $N_1, N_2, N_3, N_4$, all with different means and variances.    2. We are given the permutation $\pi = (3, 1, 2, 4)$.   3. What's $\mathbf{E}((N_1, N_2, N_3, N_4) \mid N_3 > N_1 > N_2 > N_4)$? One way to do this is by Gibbs sampling: we can sample from the conditional distribution $(N_1, N_2, N_3, N_4) \mid N_3 > N_1 > N_2 > N_4$ by doing the following:   1. Initialize some arbitrary values $x_i^0$ such that $x_3^0 > x_1^0 > x_2^0 > x_4^0$   2. For $t = 1 \ldots M$, choose a random index $i$.   3. $x_{\pi_i}^t \sim \text{TruncatedNormal}(\mu_{\pi_i}, \sigma_{\pi_i}, x_{\pi_{i-1}}^{t-1}, x_{\pi_{i+1}}^{t-1})$, with lower/upper bounding values $-\infty$ or $\infty$ if $i=1$ or $i=m$, respectively, and $x^t = x^{t-1}$ otherwise   4. Compute the mean and variance of the $x^t$ after some break-in period. The above is basically an MCMC method that walks over the conditional distribution to obtain an average computationally. The question is, is there a more efficient, faster, or better way to compute this expectation? _This question is related tothis one on math.stackexchange about computing the probability of a particular ordering, which it turns out can be done using the multivariate normal CDF (for which there exist numerical methods.)_


title: Can Mean be smaller than SD if the scores have negative numbers

text: Many of the articles discuss about the relationship between mean and standard deviation of numbers above zero stating that mean should be greater than standard deviation But if the numbers or the data has negative integers with scores ranging from -10 to +10, is it necessary that mean should be greater than standard deviation.


title: Is the logit function always the best for regression modeling of binary data?

text: I've been thinking about this problem. The usual logistic function for modeling binary data is: $$ \log\left(\frac{p}{1-p}\right)=\beta_0+\beta_1X_1+\beta_2X_2+\ldots $$ However is the logit function, which is an S-shaped curve, always the best for modeling the data? Maybe you have reason to believe your data does not follow the normal S-shaped curve but a different type of curve with domain $(0,1)$. Is there any research into this? Maybe you can model it as a probit function or something similar, but what if it is something else entirely? Could this lead to better estimation of the effects? Just a thought I had, and I wonder if there is any research into this.


title: Incorporating validation data into training set

text: Suppose that I divide my data for modelling purposes into training, testing, and validation which will then be deployed for an application (as in the response to this question). Why not incorporate the data from the validation phase to train my model again for application (or testing)? I have been told this is to avoid overfitting which makes sense, but I seem to have already incorporated data from the validation phase implicitly by my choice of the model. Is the bias that is introduced not worth the benefit of the additional data? Since I've selected my model based on the results of validation phase and I won't adjust the model parameters further, what harm is there in giving the model more data? Apologies for the vague and confused nature of the question but it has bothered me for a while.


title: Forecasting beyond one season using Holt-Winters' exponential smoothing

text: I am using the Holt-Winters' exponential smoothing technique to forecast expenditure data 2 years into the furture. The monthly data has an increasing trend and annual seasonality. I'm using MS Excel with the Solver add-in to calculate the optimal values of $\alpha$, $\beta$ and $\gamma$ to give the smallest MSE for the forecasts. The optimal values found for $\alpha$ and $\beta$ lie in (0,1) and $\gamma$ is found to be 1. I am able to calculate the forecasts for the next year (season) because the seasonals from the previous year exist. However, the forecasts for the second year are calculated to be zero because the seasonals do not exist (because _m_ is greater than 12). I have discovered that if $\gamma$ is zero, then the seasonals will be periodic, and so could be replicated after the last observed values. Is this the best way to forecast beyond one season after the last observed values? Any advice would be appreciated. Example data is below. Forecasts are needed for every month up to December 2011. I cannot see how this is possible unless $\gamma$ is zero. Numbers of Tourists   Period Month No. Tourists (Yt)   1 Jan-99 500   2 Feb-99 543   3 Mar-99 899   4 Apr-99 835   5 May-99 900   6 Jun-99 881   7 Jul-99 1154   8 Aug-99 1586   9 Sep-99 743   10 Oct-99 1104   11 Nov-99 799   12 Dec-99 560   13 Jan-00 514   14 Feb-00 665   15 Mar-00 949   16 Apr-00 975   17 May-00 924   18 Jun-00 724   19 Jul-00 1155   20 Aug-00 1541   21 Sep-00 746   22 Oct-00 944   23 Nov-00 786   24 Dec-00 652   25 Jan-01 479.4   26 Feb-01 644.4   27 Mar-01 815.8   28 Apr-01 1035.4   29 May-01 1000.9   30 Jun-01 793.8   31 Jul-01 1347.3   32 Aug-01 1378   33 Sep-01 798.1   34 Oct-01 1070.5   35 Nov-01 625.3   36 Dec-01 654   37 Jan-02 477.5   38 Feb-02 656.2   39 Mar-02 888.7   40 Apr-02 926.6   41 May-02 1000.1   42 Jun-02 1030.8   43 Jul-02 1123   44 Aug-02 1473.5   45 Sep-02 717.8   46 Oct-02 974.7   47 Nov-02 761.2   48 Dec-02 641.5   49 Jan-03 501.6   50 Feb-03 588.3   51 Mar-03 917.6   52 Apr-03 990   53 May-03 1051   54 Jun-03 764.4   55 Jul-03 1014.2   56 Aug-03 1313.6   57 Sep-03 736.3   58 Oct-03 1042.9   59 Nov-03 685.9   60 Dec-03 621.5   61 Jan-04 492.8   62 Feb-04 722   63 Mar-04 869.9   64 Apr-04 927.9   65 May-04 1028.1   66 Jun-04 883   67 Jul-04 1097.4   68 Aug-04 1398.9   69 Sep-04 834.4   70 Oct-04 1072.3   71 Nov-04 801.9   72 Dec-04 711.2   73 Jan-05 616.1   74 Feb-05 774   75 Mar-05 1088.5   76 Apr-05 956.2   77 May-05 1175.6   78 Jun-05 949.5   79 Jul-05 1120.8   80 Aug-05 1426.2   81 Sep-05 841.5   82 Oct-05 996.6   83 Nov-05 908   84 Dec-05 696.7   85 Jan-06 606.4   86 Feb-06 771.6   87 Mar-06 967.1   88 Apr-06 1235   89 May-06 1216.1   90 Jun-06 945.1   91 Jul-06 1194.4   92 Aug-06 1433.4   93 Sep-06 830.6   94 Oct-06 984.7   95 Nov-06 880.2   96 Dec-06 668.3   97 Jan-07 644.9   98 Feb-07 808   99 Mar-07 998.2   100 Apr-07 1283.9   101 May-07 1080.9   102 Jun-07 989.9   103 Jul-07 1167   104 Aug-07 1568.9   105 Sep-07 951.7   106 Oct-07 1121.4   107 Nov-07 859   108 Dec-07 660.9   109 Jan-08 647.9   110 Feb-08 911.1   111 Mar-08 1201.2   112 Apr-08 1258.1   113 May-08 1177.8   114 Jun-08 1067.6   115 Jul-08 1349.4   116 Aug-08 1702.1   117 Sep-08 982.8   118 Oct-08 1116.5   119 Nov-08 904.7   120 Dec-08 655.9   121 Jan-09 733.75   122 Feb-09 852.67   123 Mar-09 1049.88   124 Apr-09 1377.11   125 May-09 1344.05   126 Jun-09 1030.95   127 Jul-09 1242.56   128 Aug-09 1542.24   129 Sep-09 1016.42   130 Oct-09 2301.41   131 Nov-09 1138.9   132 Dec-09 1032.87


title: Creating an index that measures how accurate economsits are at predicting the US economy

text: I am interested in building an index that tracks how accurate economists are at predicting several US economic statistics:   * US Jobless claims number: comes out weekly, ranges from 100k to 400k, usually in eight year cycles   * U Michigan Consumer Confidence: comes out bi-weekly, ranges from 50 to 115.   * Case-Shiller Housing index: comes out once per month, ranges from 60 to 200 Each of these indexes have an economists estimate number and then the actual number and I am specifically interested in the accuracy of the estimates versus the actual. there are usually 30-40 economists that provide an estimate. My questions are as follows:   1. Should I use the average of the economists estimates in my index or should I use the median estimate? Usually the economists "cluster" their guesses near last months numbers while there is always the odd-ball economist who will throw out an extreme number.   2. I want to equally weight the accuracy of the estimates on each of the numbers, however the data comes out in different frequencies. If I build a spreadsheet with weekly dates, should I "fill in" the last data observation in the blank cells in my spreadsheet. For instance should I collect four different weekly US Jobless Claims numbers while repeating the same number four times for Case-Shiller (since I have no data for three weeks)?   3. When I build my index, should I calculate the % deviation from estimate vs actual and then get the z-score of that column or should I calculate the difference between estimate and actual and then get the z-score for the entire column of that difference? In the situation of Case-Shiller the housing index grew from ~60 to 183 so a 3-point difference in estimate was much bigger back then than it was when the index was at its peak. I am curious to see how you would set this up. Thanks!


title: Sample size and power calculation for randomized cluster study

text: Stuck at a pretty weird research design and would like to get some input on how to go about calculating power. A researcher would like to randomize four groups into four different interventions (one of them being the proposed new standard, the other three are currently used practices), followed by comparing the magnitude of improvement in each intervention. The chart below shows the scheme. The four groups are randomized as a whole cluster (aka group and treatment do not cross). However, the communities within each group are not randomly chosen, but rather enlisted based on logistical reason. All eligible individuals within each community will be enrolled, and I do have the predicted number of participants in the community level. We have also come up with expected differences between the four interventions. ![enter image description here](http://i.stack.imgur.com/xse0R.png) I am unsure on how to approach its power calculation. Particularly, does this design qualify as a multi-level model? What would be an ideal approach to calculate the power? Any comments/referral to other materials are welcome. Many thanks.


title: Inference and noninference parts of statistical theory

text: I was wondering from the view of dividing the topics of statistical theory into inference part and non-inference part, what inference topics and non- inference topics statistical theory is covering? By inference, I mean the task in logic to reach some conclusion from some premises. Probabilistic inference is a way of logical inference using probability. Can we say statistical theory is just providing probabilistic ways to accomplish logic inference? Or statistical theory and logic inference are overlapping, but neither falls into the other? Relating to my previous post regarding decision theory and statistical theory, how are logical inference, decision theory and statistical theory are related together and differ from each other? Thanks and regards!


title: You observe k heads out of n tosses. Is the coin fair?

text: I was asked this question with $(n, k) = (400, 220)$ in an interview. Is there a "correct" answer? Assume the tosses are i.i.d. and the probability of heads is $p=0.5$. The distribution of the number of heads in 400 tosses should then be close to Normal(200, 10^2), so that 220 heads is 2 standard deviations away from the mean. The probability of observing such an outcome (i.e. more 2 SDs away from the mean in either direction) is slightly less than 5%. The interviewer told me, essentially, "if I observe something >= 2 SDs from the mean, I conclude that something else is going on. I would bet against the coin being fair." That's reasonable -- after all, that's what most hypothesis tests do. But is that the end of the story? For the interviewer that seemed to be the "correct" answer. What I'm asking here is whether some nuance is justified. I couldn't help but point out that deciding that the coin is not fair is a bizarre conclusion in this coin-tossing context. Am I right to say that? I'll try and explain below. First of all, I -- and I would assume most people as well -- have a strong prior about coins: they're very likely to be fair. Of course that depends on what we mean by fair -- one possibility would be to define "fair" as "having a probability of heads 'close' to 0.5, say between 0.49 and 0.51." (You could also define 'fair' as meaning that the probability of heads is exactly 0.50, in which case having a perfectly fair coin now seems rather *un*likely.) Your prior might depend not only on your general beliefs about coins but also on the context. If you pulled the coin out of your own pocket, you might be virtually certain that it's fair; if your magician friend pulled it out of his, your prior might put more weight on double-headed coins. In any case, it's easy to come up with reasonable priors that (i) put a large probability on the coin being fair and (ii) lead your posterior to be quite similar, even after observing 220 heads. You'd then conclude that the coin was very likely to be fair, despite observing an outcome 2 SDs from the mean. In fact, you could also construct examples where observing 220 heads in 400 tosses makes your posterior put _more_ weight on the coin being fair, for example if all unfair coins have a probability of heads in $\{0, 1\}$. Can anyone shed some light on this for me? Edit: After writing this question I remembered that I had heard about this general situation before -- isn't it Lindley's "paradox"? Another edit: Whuber put a very interesting link in the comments: You Can Load a Die, But You Can't Bias a Coin. From page 3: > It does not make sense to say that the coin has a probability p of heads, > because it can be completely determined by the manner in which it is tossed— > unless it is tossed high in the air with a rapid spin and caught in the air > with no bouncing, in which case p = 1/2. Pretty cool! This ties into my question in an interesting way: suppose we know that the coin is being "tossed high in the air with a rapid spin and caught in the air with no bouncing." Then we definitely shouldn't reject the hypothesis that the coin is fair (where "fair" now means "having p=1/2 when tossed in the way described above"), because we effectively have a prior that puts all probability on the coin being fair. Maybe that justifies to some degree why I am uncomfortable rejecting the null after 220 heads are observed.


title: How to prove the number of distinct distributions in a group of distributions?

text: Let's say we had 5 distributions: A,B,C,D,E. An ANOVA test would tell us whether or not all of the means are equal, and thus a low p-value would mean at least one of the means is unlikely to be equal to the mean of the whole population of distributions. Is this a correct interpretation? If so, I'm looking for something similar to an iterated ANOVA test, which tests all combinations of possible distributions to determine the likeliest number of distributions. So, if all 5 of the above had distinct means, it would return 5. If all 5 did not have distinct means and would not pass ANOVA, it would return 1. If all were distinct except D and E, then it would return 4. Does such a statistical function exist with an implementation in R?


title: Recommendations for mathematical multivariate statistics with exercises

text: I need a graduate level, mathematically rigorous textbook on multivariate analysis and inference to enhance myself. I've been reading the Elements of Statistical Learning and doing the problems within, but I need a book with other focuses. Topics like famous distributions (Wishart, Wilks lambda, etc.), hypothesis testing, some theory on estimation (point, interval), and other more modern materials are welcome. I've checked this question, but the OP there was looking for something useful for his psychological analysis, which is not exactly what I want here. Currently I have the An Introduction to Multivariate Statistical Analysis, so I would also like to hear the comments about this book and its exercises. Thank you.


title: Test for equal variability in mixed model setting

text: I have a setting where I normally would model the variability in measurements by a linear mixed model which would look in R as follows.               require(lme4)     lmer(measure ~ equipment + (1|operator) + (equipment|batchid), data=mydataset)       So basically   * a fixed effect of equipment   * random intercept per operator   * random intercept per batch, differing for each type of equipment Now in this model, I would like to define a test, which checks whether the variability explained by the equipment (fixed effect + random effect) changes by equipment level. There are 2 equipment levels and a unit based on which data is collected is called a batch. Where can I find the specification of such a test?


title: ZScore threshold and low values time-series

text: Example of z-score computation: 1 - E.g.               Time-series: [0, 0, 0, 0, 1]     Current: 1     Mean: 0.2     Std: 0.44721     Z = (1 - 0.2) / 0.44721 ~= 1.7888      2 - E.g.               Time-series: [100, 100, 100, 100, 1000]     Current: 1000     Mean: 280     Std: 402.49224     Z = (1000 - 280) / 402.49224 ~= 1.7888      How can we represent that in case of example 2 - the variation is much more significant than the example 1 because of a higher mean? Is the only solution is to have z-score thresholds based on the mean value?


title: Calculating the variance of the histogram of a grayscale image

text: I am doing image processing and I want to calculate the variance of a histogram of pixel intensities. The first method I have tried: The images store the pixels values using double precision numbers, however to make a histogram, they need to be scaled so that the can be grouped in to bins. For example into bins of 0 to 255. And I get a histogram looking like this ![hist](http://i.stack.imgur.com/tArfe.png) I can then calculate variance in the usual way > $\mu = \frac{1}{n} \sum_x xf(x)$ > > $ \sigma^2 = \frac{1}{n} \sum_x {(x- \mu)^2 f(x)}$ But I have suspicions that this is not very usefull for comparing variances between different images. For example, when I get an image with a wider spread of pixel intensities, the histogram is scaled over a wider range. This means that pixels that would otherwise be grouped into a different bin, get placed in the same bin. As you can see here ![squashed](http://i.stack.imgur.com/wfXnb.png) It seems to defeat the purpose of a histogram to place most of the pixles into one or two bins. Is this a job for kernel density estimation? I need some way of comparing apples to apples.


title: Suitable analysis for my data...?

text: everyone, I'm new at the statistical analysis so I need some help from more experienced people. I'm working with biological data from gene expression analysis. I have 5 treatments (1 parental toxin; 3 compounds resulting from the metabolism of the parental toxin -let's calll them A, B and C), and a control, each one with 3 time points (6, 12 and 18h). I was trying to relate them somehow, although they are independent experiments, and summarize like this "the response of the toxin A after 12h is 23% similar to the parental toxin after 6h". So here's a list of what I've tried so far: 1. Pearson correlation/Cosine similarity. I created a matrix comparing each treatment. However, it didn't add much, because I can't turn the coefficient into a numerical valueor percentage. 2. Linear regression. Since Pearson does not discriminate between the axis, I used linear regression shifting the position of the treatment that I wanted to be the "dependent variable" (in this case, x=A,B or C toxins and y=parental toxin). I was planning to use the slope (a), as measure. If the angle a was between 0 and 22.5º degrees, I would have a response of y = 0.5 times toxin A; if the angle was 45 º, I would have y = toxin A; and if the angle >45º, I would have y = 2 times toxin. But I'm sure that's not a good idea. My supervisor insists on that, but I don't know if that's viable. What do I do with the euqation intercepts? 3. Multiple linear regression. I used the treatments that showed the highest correlation coefficients to build a multiple linear regression equation, and so I had "parental toxin 6h = 0.754 toxin A 12h + 0.246 toxin B 12h + 0.847 toxin C 18h". That's almost what I want, but, again, I feel that this is not the best way to analysze my data because it doesn't fulfill the necessary assumptions to use linear regression. 4. Binary analysis, percent agreement, logistic regression, matrix determinants, etc. So, if anyone have any suggestion, publication suggestion, anything that could shed some light to my (so far) insoluble problem, I would appreciate very much. Thanks in advance. ADDITIONAL INFO: My data is composed of 2280 genes, which have different expression values among treatments. In some treatments they are significant (expression ratio >1.5 or <-1.5 times and p-value cutoff 0.05), in others, don't. Thus, I have a lot of positive and negative values, ranging from -8.0 to +8.0, but not linearly.


title: C-Test for comparing poisson means in Scipy

text: I'm relatively new to stats, and I have a large dataset of strings (n=114541) for which I'm interested in a binary event: whether X occurs or not (specifically, whether a specific language construct is present in a given string). X occurs rarely (443/n). X is poisson-distributed. The strings in the dataset are arranged in a time series, and I'm trying to determine whether there's a significant change in the occurrence of X before and after a given moment in time. Because X follows a poisson distribution, I can easily calculate lambda before and after that moment in time (let's say that `lambda1 = 127/35465`, and `lambda2 = 316/79076`). `lambda1/lambda2` (or `lambda1-lambda2`) gives me an indication of whether the probability of X changed. What I need to know is if that change is significant (i.e., in the example above, if the difference between 0.00399 before the switch point and 0.0035 after is significant). I've read that the C-test can test for significance between two lambda values in a poisson distribution, but I'm unsure if that's the right test. And if it is: is there a way of running it in Scipy (I've looked through the documentation, and couldn't find any mention of it)?


title: How to plot a fan (Polar) Dendrogram in R?

text: I'm referring to something like this: ![alt text](http://i.stack.imgur.com/bQ8b1.png) suggested dataset for showing a solutions:               data(mtcars)     plot(hclust(dist(mtcars)))


title: Why representation of AR process comes up in estimation

text: Let ${X_t}$, $t=...-2,-1,0,1,2...$ be a stochastic process that satisfies: $X_t=\rho X_{t-1}+\varepsilon_t$ with $|\rho|<1$ and $\varepsilon_t$ is a white noise. In that case, we also know that there is a $\beta$ such that $|\beta|>1$ and another white noise $\eta_t$ such that ${X_t}$ also satisfies: $X_t=\beta X_{t-1}+\eta_t$ Suppose I have some observations of $X_t$. If I run OLS by data generated by those equivalent process, which one should I expect? Why?


title: How do I estimate a time series regression using GMM in the way proposed by Acosta-Ormaechea and Morozumi (2013)?

text: In their paper Acosta-Ormaechea and Morozumi (2013) propose a use of GMM for estimating a regression in which they try to find the impact of reallocating public expenditure from some unproductive to those productive, and I am trying to write a monograph, so I would like to know how can I use the same method of estimation, but in my particular case for time series. N.B. Sorry, my English is bad, but I hope you understood my point. I hope to find a solution as soon as possible. And if you mind leave your email so you can guide me on estimation. Alexandre Baptista - Mozambique alexvbaptista@gmail.com    **References**   Acosta-Ormaechea, S., and Morozumi, A. (2013). Can a Government Enhance Long- Run Growth by Changing the Composition of Public Expenditure? IMF Working Paper. International Monentary Fund.


title: How to deal with features only available for a few instances?

text: I am working with a regression problem. I have some features which are only available for a few instances. But with those few instances based on those features we can build a model which gives reasonable results. At the same time I have other features which are available for almost all instances. How should I combine those two type of features?


title: Explaining the difference between baseline data and subsequent experiment data if 1.00 is baseline

text: I have what I guess is a simple question really and it's probably as much about the English language as about stats. I should explain that this isn't my strong point at all but I have to write up results on some experiments I've carried out and I don't know how best to explain the way I'm comparing the results. I have a table comparing data from a number of different experiments and my baseline data. The trouble I'm having is explaining what I mean by that. If I say that my baseline data (we'll say the number is 800) represents 1.00 and my first experiment data is 600 that is .75. If my third experiment is 1200 it would be 1.50 and so on. How would I explain what is happening here? Or is there even a name for what I'm doing here? If I at least knew what to call this comparison method I could look it up Thanks for your help and sorry if its irrelevant. Any help would be really greatly appreciated Edit: Clarification The baseline result is the average result from running an experiment with no parameters 10 times as a control. Each subsequent experiment result is also an average of 10 experiments. Subsequent experiments introduce different parameters. For example Baseline = no removal Experiment 1 - remove parameter A   Expt. 2 - remove parameter B   Expt. 3 - remove parameter A+B


title: Is ARMA model possible for series with non significant ACF/PACF?

text: I am playing around with SP500 data from MASS package and from observation of the ACF and PACF there seem to be no significant autocorrelation. Now I want to model the volatility of the returns but I also want a mean equation. I'm not sure how I will get my mean equation as the ACF and PACF suggest that I can't use MA or AR terms in my mean equation. Also applying the auto.arima function from the forecast package yield the following (3,0,5), I am confused with such results as the ACF/PACF do not even show significance. I am using the box jenkins approach in my project and it wouldn't make sense to say that there are 3 AR and 5 MA terms when the ACF and PACF don't support what I say. I need to be able to explain why I am choosing this rank. Thanks


title: Maximum entropy priors for hypotheses

text: Suppose I have a random variable X which can take integer values 0 through 99. It can be, for example, numbers written on balls in a large urn. I have three hypotheses about the distribution of these numbers on the balls: $H_0$ says each number is just as likely as every other number; $H_1$ says that there's a 90% chance that the number on any given ball is between 50 and 59 including; and $H_2$ says that there's a 90% chance that the number on any given ball is 55. I want to know, based on the Principle of Maximum Entropy, what the priors for each of these hypotheses should be, and how to calculate that. It feels intuitive to me that $P(H_1):P(H_2)$ should equal $10:1$ or something like that but I don't know how to show that.


title: Differences between planning, prediction and forecasting

text: What are the differences between planning, prediction and forecasting in relation to data mining and business intelligence? It would be great if you could provide some concrete examples of these three concepts. I don't have enough experience to understand it.


title: Bootstrapping of RNA-Seq data: normal distribution?

text: I have an RNA_Seq dataset (FPKM values) obtained from the ENCODE Caltech datasets for K562 cells; I also have 6 different lists of genes and I want to show that the expression values of my panel of genes is not random. Therefore, I decided to use R and perform bootstrapping. Bootstrapping, as far as I know, requires normally distributed data. My RNA- Seq data are not normally distributed and they are strongly skewed on the left. I used BoxCox transformation:               box.cox.powers(RNA_Seq[,2])      to normalize them. I obtained a negative value (-0.03701075 to be precise) and raised my FPKM values to that. I then used:               scale()      to obtain the z-scores and see if the z-scores of my genes are bigger/smaller than +/-1.96. that's because I am interested to see if my panel of genes is more highly expressed -or it has a lower expression- than what you would expect by chance. However, the values with the **highest** FPKM became those with the **lowest** z-score and vice versa since I elevated to a negative number. I understand why this happens mathematically; however, I don't really know how to handle the bootstrap now because the data show the exact opposite of what I would reasonably expect! So, how would you guys handle this situation to test my genes of interest?   I looked up online but I could not find anything specific for RNA-Seq, or I didn't know how to apply it to RNA-Seq data. I asked this question on Stack Overflow but it was marked as off-topic and I should have posted it here. So here we go :)


title: General theorems for consistency and asymptotic normality of maximum likelihood

text: I'm interested in a good reference for results concerning asymptotic properties of maximum likelihood estimators. Consider a model $\{f_n(\cdot \mid \theta): \theta \in \Theta, n \in \mathbb N\}$ where $f_n(\mathbf x \mid \theta)$ is an $n$-dimensional density, and $\hat \theta_n$ is the MLE based on a sample $X_1, \ldots, X_n$ from $f_n(\cdot \mid \theta_0)$ where $\theta_0$ is the "true" value of $\theta$. There are two irregularities I'm interested in.   1. The data $X_1, \ldots, X_n$ are not iid and, as a result, the Fisher information about $\theta$ accrues at a rate slower than $n$.    2. $\Theta$ is a bounded set, and with positive probability $\hat \theta_n$ lies on the boundary. The boundary corresponds to a "simpler" model, and so there is particular interest in whether or not $\theta_0$ lies on the boundary.  My particular questions are   1. Letting $J_n(\theta)$ denote the observed Fisher information corresponding to $\theta$, and suppose $\theta_0$ lies in the interior of $\Theta$. Under what conditions is $$\left[J_n(\hat \theta_n)\right]^{1/2}(\hat \theta_n - \theta_0)$$ asymptotically normal as $n \to \infty$? In particular, are the regularity conditions similar to the usual ones, with the relevant modification being $J_n(\hat \theta_n) \to \infty$ in some sense?   2. Suppose instead that $\theta_0$ is on the boundary, and again recall that $\hat \theta_n = \theta_0$ happens with positive probability - for concreteness, in a mixed effects model $Y_{ij} = \mu + \beta_i + \epsilon_{ij}$ we can have $\hat \sigma_{\beta}^2 = 0$. Under what conditions does $\hat \theta_n \to \theta_0$ (almost surely or in probability) and under what conditions will $\hat \theta_n = \theta_0$ eventually (this probably fails for the mixed effects model, but corresponds to "oracle" properties for the LASSO and related estimators, so perhaps it is too much to ask for general results)? Again, just a pointer to a text with results at this level of generality would be greatly appreciated.


title: How to construct a 95% Confidence Ellipse in R

text: See the cross-post at http://stackoverflow.com/questions/16659749/getting-the- ellipse-function-to-match-the-dataellipse-function-in-r. * * * So I am sampling from a multivariate normal distribution in R, and am trying to figure out how to calculate its 95% confidence ellipse using the ellipse() function in the package car. Here is some code I am running:               mu = c(0,0)     sigma = matrix(c(20,0,0,45),nrow=2)          z = rmvnorm(10000,mu,sqrt(sigma))          par(mfrow=c(1,2))     plot(z)     ellipse(mu,sqrt(sigma*qchisq(.05,2)),radius=1)     dataEllipse(z,levels=.95)      So basically I want the ellipse command to replicate the dataEllipse command. If anyone has any suggestions that would be greatly appreciated!


title: Independent, Dependent Variables, Nonlinear

text: So, I am not a statistician, but a burgeoning mathematician. Here is something that is perplexing me in statistics, however. My apologies if the question is silly, on one hand, or just not solvable, on the other. Also, any online resources would be helpful. I am looking to learn: Let's say we have independent variables $ x_1, x_2, ..., x_n $, the dependent variable $y$ and the function $ y = f(x_1, x_2, ..., x_n) $ Say, further, that the variables $x_1, ..., x_n $ are a function of time-- $x_j = x_j(t) \forall j$ -- meaning $y = f(x_1(t),...,x_n(t))$ Given $m$ time samples of all variables, if we believe the relationship is linear, we may form the $ m \times n $ matrix $A$ whose $jth$ column is given by the values of $x_j(t)$ at the $m$ time-values. Using ordinary least- squares, we may find an $\bf{z} \in \mathbb{R}^n$ that minimizes $$ \| \bf{y} - A \bf{z} \|_2 $$ where $\bf{y}$ is the vector of time-samples of $y$. However, what if we don't believe the relationship betwen $y$ and the $x_j$s are linear (or only linear with some), or what if a linear fit is terrible? Given the data and no other information, how can one test for such a thing? What are alternatives to ordinary least squares when we don't have information about the underlying pdfs or variances? Are there functions built in R or scipy(python) for such a question? Thanks.


title: Why do I get same results for OLS and GLS in R?

text: When I run this code:               require(nlme)          a <- matrix(c(1,3,5,7,4,5,6,4,7,8,9))          b <- matrix(c(3,5,6,2,4,6,7,8,7,8,9))          res <- lm(a ~ b)          print(summary(res))          res_gls <- gls(a ~ b)          print(summary(res_gls))      I get the same coefficients and the same statistical significance on the coefficients:               Loading required package: nlme          Call:     lm(formula = a ~ b)          Residuals:         Min      1Q  Median      3Q     Max      -2.7361 -1.1348 -0.2955  1.2463  3.8234           Coefficients:                 Estimate Std. Error t value Pr(>|t|)       (Intercept)   2.0576     1.8732   1.098   0.3005       b             0.5595     0.2986   1.874   0.0937 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 2.088 on 9 degrees of freedom     Multiple R-squared: 0.2807, Adjusted R-squared: 0.2007      F-statistic: 3.512 on 1 and 9 DF,  p-value: 0.09371           Generalized least squares fit by REML       Model: a ~ b        Data: NULL            AIC      BIC    logLik       51.0801 51.67177 -22.54005          Coefficients:                     Value Std.Error  t-value p-value     (Intercept) 2.0576208 1.8731573 1.098477  0.3005     b           0.5594796 0.2985566 1.873948  0.0937           Correlation:        (Intr)     b -0.942          Standardized residuals:            Min         Q1        Med         Q3        Max      -1.3104006 -0.5434780 -0.1415446  0.5968911  1.8311781           Residual standard error: 2.087956      Degrees of freedom: 11 total; 9 residual      Why is this happening? In what cases do OLS estimates are the same as GLS estimates?


title: Overfitting when using corrected AIC for model selection

text: I am using the corrected AIC to select the lag order in a simple AR(p) model. I chose the the AICc since my sample is fairly small (n=135). The AICc minimal model is the AR(15). To me it seems like an overfit to include 15 lags in such a small sample. Do you agree? And does anyone know a rule of thumb for max lag order relative to sample size? From the `ar`function in the `stats` package in R: > `order.max`   >  Maximum order (or order) of model to fit. Defaults to the smaller of N-1 > and 10*log10(N) where N is the number of observations except for > method="mle" where it is the minimum of this quantity and 12. To me these default make little sense. Where do they come from?


title: Classification trees to pick predictive variables

text: Running a logistic regression we get p-values for all the input variables which helps us choose significant inputs. Similarly can we use the classification trees to pick variables that are split, and use those variables in the model? I think the fact that splitting the dataset on a variable leads to lower classification error should be a good indicator of the predictive power of the variable, is that true?


title: Positive log likelihood values and penalty of more complex models when ranking models using AIC

text: I am ranking candidate models using Akaike information criterion (AIC). All my models have positive -2*LL (log likelihood) values which as far as I understand is expected under certain circumstances and not much of a problem... ...unless - and this is just me wondering - I want to rank models using AIC=-2*LL+(2k) with k being the number of parameters in the model. With large positive -2LL there appears little penalty for more complex models, which suddenly rank surprisingly high. Can I still use AIC for model comparison? Thanks heaps!!


title: How to fit Bradley–Terry–Luce model in R, without complicated formula?

text: The Bradley–Terry–Luce(BTL) model states that $p_{ji} = logit^{-1}(\delta_j - \delta_i)$, where $p_{ij}$ is the probability that object $j$ is judged to be "better", heavier, etc, than object $i$, and $\delta_i$, and $\delta_j$ are parameters. This seems to be a candidate for the glm function, with family = binomial. However, the formula would be something like "Success ~ S1 + S2 + S3 + S4 +...", where Sn is a dummy variable, that is 1 if object n is the first object in the comparison, -1 if it is the second, and 0 otherwise. Then the coefficient of Sn would be the corresponding $delta_n$. This would be fairly easy to manage with only a few objects, but could lead to a very long formula, and the need to create a dummy variable for each object. I just wonder if there is a simpler method. Suppose that the name or number of the two objects being compared are variables (factors?) Object1, and Object2, and Success is 1 if object 1 is judged better, and 0 if object 2 is.


title: auto.arima warns NaNs produced on std error

text: My data is a time series of employed population, L, and the time span, year.               n.auto=auto.arima(log(L),xreg=year)     summary(n.auto)     Series: log(L)      ARIMA(2,0,2) with non-zero mean           Coefficients:              ar1      ar2      ma1     ma2  intercept    year           1.9122  -0.9567  -0.3082  0.0254    -3.5904  0.0074     s.e.     NaN      NaN      NaN     NaN     1.6058  0.0008          sigma^2 estimated as 1.503e-06:  log likelihood=107.55     AIC=-201.1   AICc=-192.49   BIC=-193.79          In-sample error measures:                ME          RMSE           MAE           MPE          MAPE      -7.285102e-06  1.225907e-03  9.234378e-04 -6.836173e-05  8.277295e-03               MASE       1.142899e-01      Warning message:     In sqrt(diag(x$var.coef)) : NaNs produced      why does this happen? Why would auto.arima selects the best model with std error of these ar* ma* coefficients Not a Number? Is this selected model valid after all? My goal is to estimate the parameter n in the model L=L_0*exp(n*year). Any suggestion of a better approach? TIA. data:               L <- structure(c(64749, 65491, 66152, 66808, 67455, 68065, 68950,      69820, 70637, 71394, 72085, 72797, 73280, 73736, 74264, 74647,      74978, 75321, 75564, 75828, 76105), .Tsp = c(1990, 2010, 1), class = "ts")     year <- structure(1990:2010, .Tsp = c(1990, 2010, 1), class = "ts")     L     Time Series:     Start = 1990      End = 2010      Frequency = 1       [1] 64749 65491 66152 66808 67455 68065 68950 69820 70637 71394 72085 72797     [13] 73280 73736 74264 74647 74978 75321 75564 75828 76105


title: Random generation of scores similar to those of a classification model

text: Hello fellow number crunchers I want to generate n random scores (together with a class label) as if they had been produced by a binary classification model. In detail, the following properties are required:   * every score is between 0 and 1   * every score is associated with a binary label with values "0" or "1" (latter is positive class)   * the overall precision of the scores should be e.g. 0.1 (<- parameter of the generator)   * the ratio of scores with label "1" should be higher than overall precision in the top-section and lower in the bottom section (<- the "model-quality" should also be a parameter of the generator)   * the scores should be in such a way, that a resulting roc curve is smooth (and not e.g. that a bunch of scores with label "1" are at the top and the rest of the scores with label "1" is at the bottom of the list). Does anyone have an idea how to approach this? Maybe via generation of a roc- curve and then generating the points from that cure? Thanks in advance!


title: How do I study the correlation between two categorical variables?

text: Say I have 1000 unique websites in my sample. Each website has two categorical variables assigned to it:   * PolicyDomain (e.g. 0 = "health", 1 = "defence", etc... up to 42).   * Community (e.g. "Community1", "Community2" ... up to "Community 6"). How do I find out whether there is a correlation between PolicyDomain and Community in my sample? That is, is there a relationship between a website's PolicyDomain versus its Community? Are some PolicyDomains more likely to be in certain Communities?


title: Classification problem using imbalanced dataset

text: I am working on a pattern identification/classification problem on an imbalanced dataset, with target to non target proportion in population approx as 1%:99%. There are around 0.5 million records in my dataset. I am restricted to use SAS E-Miner for this analysis. Currently I am using the following approach:   1. Give appropriate decision weights (profit matrix)    2. Undersample on the majority class or "good" records.   3. Running a decision tree on the sample.  My questions are:   1. How can I undo the effect of undersampling of majority class?   2. Does giving appropriate decision weights actually helps in removing the bias introduced because of undersampling? Or are these really two independent things.   3. Even if these decision weights are applied, how do we determine the optimal decision threshold for basing our decisions. I have tried a boosting algorithm (without adjusting prior probabilities and without using decision weights) but the number of rules/patterns that get thrown up are around 20+ which seems like a mild concern to me. Would appreciate any inputs from CV community folks.


title: Kaplan Meier survival estimate at t=infinite

text: Since the KM estimator does not allow for any other state other than the event of interest is it safe to assume that $S(t=\infty) = 0$ ? Or can it be that $S(t=\infty) > 0$


title: How to measure agreement between a set of raters?

text: I am running an experiment where I have a finite set of raters and a finite set of items, and raters have to provide their subjective judgment about each item. The goal is to measure the importance of those items for the raters. For every item, each rater uses a Likert-like scale (1= Unimportant, 2= Of Little Importance, 3=Moderately Important, 4=Important, 5=Very Important) Knowing that the judgments are subjective I want to measure how raters agree in their ratings, and eventually observe new patterns in theirs judgments. The question is: Which statistical method/tool is more appropriate for such an analysis?


title: Multidimensional quantiles

text: I have 1000 observations with 2 continuous variables : > Observation ID | X | Y


title: Evaluating the quality of a classifier which provides probability per class

text: I am looking for a metric to evaluate the quality of a classifier which answers the question for a class of an object with a probability per known class. Say I use the classifier to classify an object of class A and the classifier classifies the object with 75% probability as A and 25% probability as C. One option is, just to pick the probabilities, in this case the 75%. Are there any other metrics? The Classifier uses a database of classes with different objects but not all objects of one class share the same properties found in an query object, so not all objects of a class may be relevant.


title: Can the Mantel test be extended to asymmetric matrices?

text: The Mantel test is usually applied to symmetric distance/difference matrices. As far as I understand, an assumption of the test is that the measure used to define differences must be at least a semi-metric (meet the standard requirements of a metric but not the triangle inequality). Can the assumption of symmetry be relaxed (giving a pre-metric)? Is it possible to apply the permutation test in this case, using the full matrix?


title: Fit a distribution to a combinatorial problem

text: In my previous question titled "Conditional combinations of balls in bowls", is there a distribution to fit $k$ when $d \gg M$? I mean, when $d$ is so large, what is the distribution of the total number of balls? $M$ can be a number between $4$ to $64$ and $d$ (number of bowls) is very large, about $10000$. Each bowl can have $0,1,2,\ldots,M$ balls in it with equal probability. What is the distribution of total number of balls $k$? $\Pr(k=k_0)=$?


title: Survey Weighted Random Effects Logit Model in R

text: I am trying to predict a binary outcome with a model that includes a random effect using survey data. I've included a description of the sampling design below, so feel free to comment on my survey weighting approach. My primary question is how to include a random effect in the survey weighted model. Here is the code up to this point:               # Libraries     library(survey)     # Make dataframe object where d is the working dataframe     dfobj <- svydesign((id = ~cluster+household, strata = ~interaction(region, urban)  weights = ~chweight, strata = ~strata, data = d)          # Run a logit model     formula1 <- stunting ~ modern_toilet + diarrhoea + fever + insurance + sex + age + region_code     model1 <- svyglm(formula=formula1,design=dfobj,family = quasibinomial)      I would like the random effect to be on the region. Thanks, Sampling Description: The MICS 2006 used a two-stage stratified sample design. At the first stage of sampling, 300 census enumeration areas (124 urban and 176 rural EAs) were selected. These are a subsample of the 660 EAs (281 urban and 379 rural) selected for the GLSS 5. The clusters in each region were selected using systematic sampling with probability proportional to their size.


title: Kolmogorov-Smirnov test using half the sample to estimate mean and variance and testing against the other half moments for normality

text: using a monte carlo simulation Im just trying to test the null hypothesis for normality in a sample that splits N values into 2, using the first N/2 values to estimate the sample mean and variance, and the second N/2 values to apply the test based on these moment estimates (the Normal distribution tested against will use the sample mean and variance estimated from the first N/2 values). I am using a 95 % confidence interval and when I do 10000 iterations and test all their p values I am getting an acceptance rate of approx 85% and was wondering could some one explain the reason why this is not 95~% and what the reason for the rejection rate is?


title: High order Markov Chain transition probability

text: I'm modelling a time series using high order markov chain and estimate the transition probability of the equivalent first order markov chain. The problem is that some of transition in the equivalent first order markov chain from state i to state j doesn't occur in the training dataset at all. How to handle this problem and give a estimate on this transition probability? Thanks a lot!


title: About linear combinations of independent Brownian motions

text: Let $B(t)$ and $W(t)$ be two independent Brownian motions. Show that $X(t)=\frac{B(t)+W(t)}{\sqrt2}$ is also a Brownian motion. Find the correlation between $B(t)$ and $X(t)$.


title: population statistics and deviation from mean

text: Suppose $X \sim N (50, 20.5)$ represents a population. If $X_1, X_2, …, X_n$ is a random sample from this population, find $n$ such that $\mathbb{P}(X_1 + X_2 + … + X_n > 2000) = 0.95 $


title: Papers on resource prediction over time series with asymmetric information?

text: I'm looking for resources to help me solve an issue of resource prediction. We can make a number of observations of the resource over time, but the way in which the resource changes is affected by both us and other players, furthermore the amount of resource as it appears to us can be false in some circumstances. I'm looking for ways of being able to sort the signal from the noise and be able to build some kind of understanding in terms of how the resource changes over time in terms of a model. Imagine, but on a much larger scale.   1. Resource = 10   2. Resource = 9   3. Resource = 9   4. Resource = 8   5. Resource = 8 (We purchase 2 lots)   6. Resource = 6 (We know that it has definitely decreased by 2)   7. Resource = 5   8. Resource = 5   9. Resource = 6 (Okay, there must have been at least one more than is being reported) 10.Resource = 3 11.Resource = 0 (There may be more in reality, but we still could not purchase regardless) 12.Resource = 0 13.Resource = 1 (Okay, there were definitely 2 more than being reported) 14.Resource = 0


title: HIstogram of oriented gradients (HOG) features descriptor theoretical problems

text: I'm going to implement HOG as my features descriptor.   But there are some things that make me confused: For example:   If we have an image with size of 10 x 20   1. If we want to compute the HOG of that image (not overlay), and we determine that the   cellsize would be 1x1 and the blocksize will be 2x2. So I will have a block with size 4 (2x2= 4 cells x 1 pixels = 4). And then I will have 50 blocks. Is that right?   2. And the features that I get will be computed for each cell / block? If it's for each cell, so I will have features of 4x50=200? The next questions will be based on this link                1. "The same function can also be used to generate a pictorial rendition of the features,          although this unavoidably destroys some of the information contained in the feature          itself. To this end, use the render command:"              imhog = vl_hog('render', hog, 'verbose') ;         clf ; imagesc(imhog) ; colormap gray ;              Why they are using render? Also what does "verbose" mean?                2.  "In this case the feature has 31 dimensions. HOG exists in many variants. VLFeat          supports two: the UoCTTI variant (used by default) and the original Dalal-Triggs          variant (with 2×2 square HOG blocks for normalization). The main difference is that          the UoCTTI variant computes both directed and undirected gradients as well as a four          dimensional texture-energy feature, but projects the result down to 31 dimensions.          Dalal-Triggs works instead with undirected gradients only and does not do any          compression, for a total of 36 dimensions."            What do they mean by normalization, directed and undirected gradients as well as a four        dimensional texture-energy feature? And how do they perform compression?          3. Is 'Flipping HOG from left to right" that they mean is overlay the blocks?          4. "Another useful option is BilinearOrientations switching on the bilinear orientation            assignment of the gradient (this is not used in certain implementation like            UoCTTI)."             What does BilinearOrientations mean?


title: Do I need a planned contrast or post-hoc analysis when doing LMM in R?

text: I am doing a LMM in R and would like to know if I need to do a planned contrast or a post-hoc analysis. From my understanding, the LMM in R already provideds me a planned contrasts and if I have many types in each variable, you can just change the baselines of model to see all angles of it.


title: Modelling points in basketball match as a Markov chain with Poisson regression

text: I am modelling a simple model of how many points are scored in total in basketball matches. I have a dataset of 500 matches. Pre game I know which team are favorite to win, for all the matches in my dataset. For every minute of the matches, I observe the difference in score. The difference in score is a factor: X for draw, H for home team leading and A for away team leading. The favorite team variable is also a factor, either Home or Away. I actually have specified more factors in my model for different types of home and away lead, but for simplicity this is how I present it. So I want to use Poisson regression to model how many points are scored, with difference in score and which team are favorite as explanatory variables: $\log(\mu)=\log(t_x)+\beta_0 +\beta_1\text{Favorite} +\beta_2\text{Diff_score} +\text{interaction_favorite/diff_score}$, or,               basketball<-glm(points~favorite*diff_score, offset=log(t_x) family="poisson", data=data),      where $t_x$ is the offset. So Y is scaled with a time vector $t_x$ for minute $x$, such that we are looking at points scored after minute $x$ until end of match. For example, if $t_4$=1.1, and points scored after minute 4 is 120, then the observation for minute 4 would be ($120/1.1,X_4$), where $X_4$ is the state of the match (explanatory variables) at minute 4. More specifically, I treat the basketball matches as Markov chains. For every minute, of every basketball match, I look at the state in which the game is in (favorite, score difference) and look at the outcome at the end of the game. So, for every game, I have 48 observations (48 minutes in a basketball match) and 48 outcomes (points) on which I perform the poisson regression. For example, if I look at just match number 1: $Y/t_1$ is the number of points at the end of the match, scaled with some $t_1$, and $X_1$ is the state at minute 1. I have this set of observations (for just match number 1): ($Y/t_1$,$X_1$),($Y/t_2$,$X_2$),...,($Y/t_{48}$,$X_{48}$). I do this for all the matches, and perform poisson regression on 48*500 datapoints. Obviously the states of the matches are dependent on each other, but I include this with the interaction term in the glm function in R. And I view the outcome (points at the end) as only dependent on the state which you are currently in. But I am not sure if I am "allowed" to do this. ![Here are some plots of my results](http://i.stack.imgur.com/E2gu6.png) The middle and right plot seems to look "OK", but the left plot worries me... Thank you


title: Estimating distribution given top k order statistics and unknown n

text: This problem occurred to me a couple of days ago, in the context of a game with a leaderboard. I wondered, given only the leaderboard, could I estimate parameters for the distribution of scores? Estimation of order statistics is a reasonably well understood problem, but as far as I can see all require knowledge of n, the total number of games giving rise to the leaderboard. If we know properties of the order statistics (for example, if we know that they tend to a uniform distribution) we might estimate the parameters of the order statistics distribution and work back from there. However, this feels rather unsatisfying.


title: Subset models in auto.arima function in forecast package

text: I wanted to ask whether it was possible to use the auto.arima function to identify subset ARIMA models rather than those of pure lags? I have identified a model in Stata in subset lags that performs well and wanted to cross check this with the `auto.arima()` function but I can't seem to figure out if subset lags are supported.


title: How will you deal with "don't know" and "missing data" in survey data?

text: As title, I am thinking of merging both into "missing data", which is to name it as NA in R. Since I don't see it will make much sense (or even any sense), to separate the "don't know" row out and to compare the information with other rows. Is it OK for me to do so? Thanks!


title: Justification of aggregation in multi-level model - ICC in SPSS

text: I am running a multi-level model on data which is set out with one row per individual survey return, each with Individual ID and Organisational ID as columns. There are variables which can be aggregated in to group, but I first have to justify this aggregation, and believe I have to use ICC(1), ICC(2) and Rwg. Although I can find a lot explaining what these are, I cannot seem to figure out how to do this with my output in SPSS. I at first believed it came from the output of a one-way ANOVA, but have since been told it should be two-way for ICC(2) as this is a two-level multi-level model. Can anyone help me on this one? I can't see how I can do a two-way ANOVA as I surely only need to use the Organisational ID and aggregating variable in question and two-way requires three variables? Any help would be very gratefully received. Thanks!


title: Checking whether accuracy improvement is significant

text: Suppose I have an algorithm that classifies things into two categories. I can measure the accuracy of the algorithm on say 1000 test things -- suppose 80% of the things are classified correctly. Lets suppose I modify the algorithm somehow so that 81% of things are classified correctly. Can statistics tell me anything about whether my improvement to the algorithm is statistically significant? Is the concept of statistical significance relevent in this situation? Please point me in the direction of some resources that might be relevant. Many Thanks.


title: Probability of Unique Minimum (Discrete)

text: This is a discrete problem concerning integers. If there are $n$ independent random variables $X_1,...,X_n$ that each take on a value from $\{1,...,x\}$ uniformly at random ($x$ distinct values), what is the probability that the minimum (call it $l$) is _unique_? The following are approaches I've taken: **Lower Bound:**   Treat the problem as a birthday paradox problem. This is a lower bound because the birthday paradox tells us the probability that no _pair_ share the same birthday. We only care if anyone shares a birthday with the _earliest_ birthday (whatever the earliest might be)... **Upper Bound:**   Sort the values in ascending order and select the first entry ($l$). Note: there may be ties. All $n-1$ remaining values shouldn't match $l$ if it's a unique lowest. This is trivial if we ignore the fact that selecting the lowest value $l$ gave us information about the remaining values (namely, they are "squished" into the range $\{l,...,x\}$) $$ p=n(\frac{x-1}{x})^{(n-1)} $$ ...But we can't actually ignore that fact so this is an upper bound. **Exact?**   The upper bound calculation nudged me into the direction of the exact solution but I can't simplify it past: $$ \frac{n}{x}\sum_{j=1}^{x}(\frac{x-j}{x-j+1})^{(n-1)} $$ Too be clear, here is what each variable means:   x = number of discrete values that each RV might take on   j = the assumed value of the minimum for that iteration   n = the number of variables I got to this form by starting with the law of total probability for a unique minimum given that minimum is $j$. The probability that the minimum is $j$ is a constant with respect to $j$ and can be pulled out. The probability that there is a unique minimum given that the minimum is $j$ is what's left inside the summation. I very well may have made a mistake somewhere. Example:   x = 10. Values drawn from $\{1,...,10\}$   n = 4. Values: $[3, 5, 3, 10]$ (no distinct minimum)


title: Bootstrapping a sample with unequal selection probabilities

text: I want to "blow up" a sample, taken with replacement, for which I know the overall sampling probability $\pi_i$ for each item $i$. Is it valid to use bootstrapping and apply inverse probability weighting during the selection (as in the Horvitz-Thompson estimator: weight each item with $1/\pi_i$), or are there any pitfalls? A quick search on Google suggests that this hasn't been fully investigated yet, and the `boot` package in R allows weights but does not comment on where they are supposed to come from. The purpose of "blowing up" is, among others, to be able to resample _with uniform probability_ from the blown-up population.


title: Analysing over 300,000 rows in Excel to make pretty graphs

text: I'm doing a research module for my Computer Science degree, and for my topic, I have collected over 500,000 tweets using the Twitter Streaming API, using a ruby script to store them in a Mongo database (BSON/JSON). I started recording the tweets on Tuesday 7th Feb, and stopped the following Tuesday, so there is a week's worth of tweets. Here is what the spreadsheet looks like. I have successfully exported around 300,000 tweets to an excel spreadsheet (I can hear groans already). I would like to make some time series charts, for example volume of tweets over time and eventually include followers_count as a weighting. But I'm unsure as to how I would calculate this. I think I need to make the created_at column more meaningful to excel but converting it to a date/time it can understand. I've also had a go with Rapid miner and managed to import a spreadsheet and convert the `created_at` field into something the program can understand, but I didn't really have any idea what I was doing after that! I'd really appreciate some hints as I'm a bit stuck right now.


title: When parsing text for n-grams - should punctuation be included?

text: I want to start working on data-mining by parsing text. It seems like the best place to start is by processing n-grams out of text to try sentiment analysis.               Muffins are fine, I wouldn't say I like them though.      However, I'm curious to know if I should include punctuation or not. (I plan on starting with 3-grams and working up since I'm not sure 2-grams include enough information for accurate results.)               Muffins are fine | are fine [,] | I wouldn't say | ....      Since a "," was found, start over at next word after the ",". Instead of including the punctuation like normal.               Muffins are fine | are fine , | fine , I | , I wouldn't | ...      Can anyone tell me if this is a bad idea?


title: How can I simulate new observations to a fitted model object while specifying the covariates in R?

text: I have a model that is fitted with glm and I want to be able to calculate the response of this model when I specify the independent variables, while simulating the error terms. For example, if my formula is `y~x1*x2`, then this corresponds to: $y_{i} = \mu + \beta_{1} x_{1i} + \beta_{2} x_{2i} + \beta_{3} x_{1i} x_{2i} + \epsilon_{i}$ where $\epsilon_{i}$ is an error term. I want to be able to simulate $y^{\star}_{i} = \hat{\mu} + \hat{\beta}_{1} x^{\star}_{1i} + \hat{\beta}_{2} x^{\star}_{2i}+ \hat{\beta}_{3} x^{\star}_{1i}x^{\star}_{2i} + \epsilon_{i}$, with $\hat{\mu}$, $\hat{\beta}_{1}$, $\hat{\beta}_{2}$, and $\hat{\beta}_{3}$ being the values estimated by glm and $x^{\star}_{1 i}$ and $x^{\star}_{2 i}$ being the specified covariates / independent variables. This is similar to what predict() does, however, it sets $\epsilon_{i}=0$. It looks like what I want to do is similar to the functionality provided by simulate() and predict(). Although it doesn't look like simulate() allows you to specify the covariates. What is the best way to do this?


title: Variable importance for SVM regression and averaged neural networks

text: Is there any way to calculate variable importance in R for SVM regression and averaged neural networks? I've been using **caret** package, that has **varImp** function in it               > m <- best.tune(svm, train.x = descr[rownames(tr[[i]]),2:ncol(descr)],      train.y = tr[[i]][,1],  data = df, cost = 2^(seq(0,10,5))),      tunecontrol = tune.control(sampling = "cross"))     > varImp(m)      Error in UseMethod("varImp") :       no applicable method for 'varImp' applied to an object of class "svm"      According to the developer, this approach wasn't realized for SVM method However, **rminer** package suggests such function as _Importance_. Though, it throws an error:               VariableImportance = Importance(svmFit, data=descr[rownames(tr[[i]]), 2:ncol(descr)],                                         method="1D-SA")     Error in Importance(svmFit, data = descr[rownames(tr[[i]]), 2:ncol(descr)],  :      duplicate 'switch' defaults: 'lm == func...' and 'NULL'


title: Fitting a non-linear random effects model with binomial error structure

text: I've got a non-linear model that I've been applying to some data of repeated binary outcomes. I have data for multiple years, and I'd like to add random effects (by year) for two of my parameters. Looking at the `lme4` package in R, it seems that it supports non-Gaussian error structure and nonlinear models, but not simultaneously. Is there a way around this, or another package I could be using? I hear ADMB has a steep learning curve, so I'd prefer not to tackle that. Without random effects, my model is quite simple--only 4 parameters and I've been able to fit it using a formula call to the `mle2` function from the `bbmle` package.


title: Binomial trials required to achieve a given probability of at least one success

text: If p is the probability of success of a binomial trial i would like to calculate the number of trials n required that if performed would give a probability x of at least one success. Is there a way to obtain this n in R? I have a large vector of probabilities and would like to extract this n for a given x. Thank you


title: Cross validation to test the performance of two spatial interpolation methods

text: I have 14 weather stations' temperature data for the period between 2010–2013. I need to evaluate the performance of two spatial interpolation methods. I suggest to select 10 days from this period of temperature data and estimate temperature data for these days to each station and then calculate the RMSE of these days in both methods. Does this mean I perform a cross validation test using RMSE?


title: How to extract the importance of predictors in Adaptive boosting?

text: I have trained a boosting classifier using `ada` package in `R`. Now I want to see the importance of my predictors in constructing the model. I have used `varplot` and got the following plot of variable importance ![enter image description here](http://i.stack.imgur.com/bgPya.png) I have also used `listAdaVarsUsed` to extract the number of times variables have been used in the model building process:                              tx         total.profit                  p1x                   102                   65                   37                     T            GR.change            DY.change                    29                   24                   18          loyal.member            FV.change            MT.change                    14                   13                   10             BH.change            DL.change Socio.Economic.level                     6                    4                    2       > The question is, which of these two measures can show me the importance of > my predictors? As you see for instance `FV.change` is considered to be an > important variable according to the plot but it has been used only 13 times. > How this is possible? Your help would be highly appreciated.


title: What is autocorrelation function?

text: Can somebody explain autocorrelation function in a time series data? Applying acf to the data, what would be the application?


title: Comparing nested binary logistic regression models when $n$ is large

text: To better ask my question, I have provided some of the outputs from both a 16 variable model (`fit`) and a 17 variable model (`fit2`) below (all predictor variables in these models are continuous, where the only difference between these models is that `fit` does not contain variable 17 (var17)):               fit                    Model Likelihood     Discrimination    Rank Discrim.                                  Ratio Test            Indexes          Indexes             Obs        102849    LR chi2   13602.84    R2       0.173    C       0.703           0          69833    d.f.            17    g        1.150    Dxy     0.407           1          33016    Pr(> chi2) <0.0001    gr       3.160    gamma   0.416          max |deriv| 3e-05                          gp       0.180    tau-a   0.177                                                     Brier    0.190                      fit2                 Model Likelihood       Discrimination    Rank Discrim.                                  Ratio Test            Indexes          Indexes             Obs        102849    LR chi2   13639.70    R2       0.174    C       0.703           0          69833    d.f.            18    g        1.154    Dxy     0.407           1          33016    Pr(> chi2) <0.0001    gr       3.170    gamma   0.412          max |deriv| 3e-05                          gp       0.180    tau-a   0.177                                                     Brier    0.190                I used Frank Harrell's `rms` package to build these `lrm` models. As you can see, these models do not appear to vary much, if at all, across _Discrimination Indexes_ and _Rank Discrim. Indexes_ ; however, using `lrtest(fit,fit2)`, I was provided with the following results:                L.R. Chisq         d.f.            P      3.685374e+01     1.000000e+00    1.273315e-09       As such, we would reject the null hypothesis of this likelihood ratio test; however, I would assume this is likely due to the large sample size ( _n_ = 102849) as these models appear to perform in a similar fashion. Furthermore, I am interested in finding a better way of formally comparing nested binary logistic regression models when _n_ is large. I greatly appreciate any feedback, R scripts, or documentation that can steer me in the right direction in terms of comparing these types of nested models! Thanks!


title: High concordance in cox PH model even though PH assumption is violated

text: I am using Cox proportional hazards models for survival analysis. The specific reason I am interested in them is because they give a nice way to quantify effect size between groups via the hazard ratio, assuming the PH assumption is not violated. I use R's `survival` package for the modeling. The models I am making have no interaction terms between predictors. I am currently in the following situation for some models:   * the model is a good fit, with concordance over 80% (obtained via `summary` of the `coxph` fit).    * the PH assumption is violated for one or more predictors (e.g. very small $p$-values for some predictors in `cox.zph`). This leads me to a couple of questions:   1. What is the main consequence hereof?    2. Is there still some value in the hazard ratios for predictors for which the PH assumption is invalidated?   3. How can the fit be so good despite a key assumption being violated?


title: Discrete unemployment duration with additional retrospective information

text: I am trying to run a discrete duration (logistic) regression on unemployment. However, I face the following challenge: The survey collects information on whether the person is employed or not for ten consecutive months. In addition, at the first interview those who are unemployed are also asked about their past unemployment duration. As such, after restricting my data set to those who have experienced unemployment, I know the labor force status of individuals for up to ten months (as some of them would drop from the sample due to attrition) plus for the ones that were unemployed at the first interview their past unemployment duration. I also have information about the education, marital status, as well as other key covariates. The challenge is to model this duration. I thought of running a discrete duration analysis, where among the covariates I include dummies for the ten consecutive months plus some explanatory variables, such as age, occupation, etc. (some of which are time-varying). For the ones that became unemployed during the interview period, I know when they started their unemployment spell and when they exited it or were censored, so my model should do well with them. The issue is with the ones that were already unemployed during the first wave, as they have already been unemployed when they entered the interview period. As I have information about their unemployment duration, I could potentially recalculate their duration from the time they became unemployed. The issue is that I do not observe some of the key covariates for this retrospective period and especially for the time-varying covariates this may be a problem. I am inclined to go with the initial model (where I include dummies for the 10 periods I observe) and then somehow include the information about the past employment duration for those who were already unemployed during the first wave of the interview in the model, but I am struggling to find a promising way to do that. I thought to include another regressor which indicates the initial unemployment duration status, but I am not convinced how this would work. I highly appreciate any suggestions.


title: Log-linear or poisson model with R

text: I have a data.frame (`myData`) with 6 variables which are:               Fish.species | anemone.species | substrate | Depth | Turbidity | anemone.size      Some of these variables are categorical (`Fish.species`, `Anemone.species`, `substrate`) and some are ordinal and quite a lot of different values (`Depth`, `Turbidity`, `anemone.size`). I'd like to know:   1. Which (if any) of the 5 variables has an effect on `Fish.species`   2. Which (if any) of the 5 variables has an effect on `Anemone.species` I'll focus only on the first point as I think the second point is essentially the same analysis reversing simply the two variables. To answer this question, I should use regression for count data. Somehting like a poisson or a log-linear model. Is it correct? What about:               myTable = xtabs(~Fish.species + anemone.species + substrate + Depth + Turbidity + anemone.size, data=myData)      Note: There are some 0 counts that might be an issue I guess! There are even more zeros than non-zero when looking at               table(Fish.species, Anemone.species, myData)      Should I cut my continuous variables like this?               Depth = cut(Depth, breaks=c(0,3,6,9,12,15,18,21,24))      Note: There would still have many zero counts!               loglm(~Fish.species + anemone.species + substrate + Depth + Turbidity + anemone.size, data=myTable)      I get only one p.value and I don't really get what it means and how I can know what variables has an effect on `Fish.species`. Can you help me with that? what about:               glm(Fish.species ~ anemone.species + substrate + Depth + Turbidity + anemone.size, data=otherData, family=poisson)      where `Fish.species` is not a categorical variable but a count as it was in the `loglm` function Would it make more sense?


title: Correct model to use in travel speed analysis

text: I am attempting to analyse downstream travel speed of fish within two different rivers, however, I am getting lost on which model I should use. Please see working below data:               data4<-structure(list(code = structure(1:32, .Label = c("10888", "10889",      "10890", "10891", "10892", "10893", "10894", "10896", "10897",      "10898", "10899", "10900", "10901", "10902", "10903", "10904",      "10905", "10906", "10907", "10908", "10909", "10910", "10914",      "10916", "10917", "10919", "10920", "10922", "10923", "10924",      "10925", "10927"), class = "factor"), speed = c(0.0296315046039244,      0.0366986630049636, 0.0294297725505692, 0.048316183511095, 0.0294275666501456,      0.199924957584131, 0.0798850288176711, 0.0445886457047146, 0.0285993712316451,      0.0715158276875623, 1.53078364415863, 1.32386730007184, 0.0879404548120395,      0.0884436446391665, 0.142680147253483, 0.256122292011316, 0.029228002169911,      0.0535706245570686, 0.0268752274519418, 0.238146999654007, 0.424175803745118,      0.021181456637386, 0.0337364380380858, 0.030136111418826, 0.0335130031085068,      0.148793370546863, 0.096737574741185, 0.114746050390462, 0.0993760102817675,      0.0985453244761994, 0.100223001228893, 0.0924775724221547), meanflow =       c(0.657410742496051,      0.608271363339857, 0.663241108786611, 0.538259450171821, 0.666299529534762,      0.507156583629893, 0.762448863636364, 37.6559178370787, 50.8557196935557,      31.6601587837838, 225.109658536585, 230.452020833333, 106.450571428571,      106.51658356546, 26.6998648648649, 22.6653709677419, 0.667544935064935,      0.803912464319696, 0.660643403441683, 58.1045714285714, 164.84368,      63.5970920920921, 50.7789043570669, 47.3644579971523, 45.6948401055409,      0.808431216931217, 0.612270689655172, 46.4646057866184, 0.605024778761062,      0.6091, 0.601483928571429, 0.616563426688633), length = c(136,      157, 132, 140, 135, 134, 144, 149, 139, 165, 155, 143, 144, 143,      156, 147, 133, 131, 129, 143, 140, 125, 140, 140, 147, 120, 115,      139, 123, 120, 125, 136), river = structure(c(2L, 2L, 2L, 2L,      2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L,      1L, 1L, 1L, 1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 2L), .Label = c("c",      "f"), class = "factor")), .Names = c("code", "speed", "meanflow",      "length", "river"), class = "data.frame", row.names = c(2L, 4L,      6L, 8L, 10L, 12L, 14L, 16L, 18L, 20L, 22L, 24L, 26L, 28L, 30L,      32L, 34L, 36L, 38L, 40L, 42L, 44L, 48L, 50L, 52L, 54L, 56L, 58L,      60L, 62L, 64L, 68L))       I would like to look at the effect of specific variables on travel speed, particularly to see if travel speed is different between rivers. Fixed effects: river, length Interaction: An interaction between river and length is considered and added as a fixed effect due to fish swimming in slower currents, and body length having a direct impact on swimspeed. currents differ between river and so length will have greater impact on rivers with slower currents So:               model1 <-lm(speed ~ 1 + river + length + river: length, data = data4)      However 'meanflow' has a direct impact on travel speed as a higher flow will enable fish to swim quicker downstream. To account for this I believe this should be added as a random effect?! And be like:               library(lme4)     model2<- lmer(speed ~ river + length + river:length + (1|meanflow), data = data4)      Yet this does not run as 'meanflow' is a continuous variable. would I be better to dis regard this variable, or maybe use a different model all together Any help would be great


title: Different factor levels in R

text: I am running a model using `lm` in R to look at the effects of two independent factors. These independent factors are breeding group (`bgrp`: 2 levels with `bgrp1`=maternal lines; `bgrp2`=terminal lines) and sex (`psex`: 2 levels with `psex1`=males; `psex2`=females); and using number of total born (`tb`) as a linear covariate. I have around 600 records altogether. The model is:               mybgrp <- lm(formula = prtemp ~ bgrp+psex+tb-1, data = apirt)  #removing the intercept     mybgrp          Call:     lm(formula = prtemp ~ bgrp + psex + tb - 1, data = apirt)          Coefficients:        bgrp1     bgrp2     psex2        tb       37.72592  37.98272   0.14467   0.01982          > summary(mybgrp)          Call:     lm(formula = prtemp ~ bgrp + psex + tb - 1, data = apirt)          Residuals:         Min      1Q  Median      3Q     Max      -5.6877 -0.2454  0.0768  0.3916  1.6561           Coefficients:            Estimate Std. Error t value Pr(>|t|)         bgrp1 37.725922   0.135486 278.449  < 2e-16 ***     bgrp2 37.982716   0.129558 293.171  < 2e-16 ***     psex2  0.144669   0.055140   2.624  0.00891 **      tb     0.019818   0.009342   2.121  0.03429 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.6888 on 621 degrees of freedom       (5 observations deleted due to missingness)     Multiple R-squared: 0.9997,     Adjusted R-squared: 0.9997      F-statistic: 4.788e+05 on 4 and 621 DF,  p-value: < 2.2e-16      Now, up to this stage I could see the values of the two levels of the `bgrp` (`bgrp2=37.7`; `bgrp2=38`), but not too certain how to use the value of `psex2` (0.144669) or any other values in the table of coefficients to calculate the values of `psex1` and `psex2`. After asking around, I changed the order of the independent from the first model to this the one below, hoping that it would give me the values of `psex1` and `psex2`:               mypsex <- lm(formula = prtemp ~ psex+bgrp+tb-1, data = apirt)     mypsex     Call:     lm(formula = prtemp ~ psex + bgrp + tb - 1, data = apirt)          Coefficients:        psex1     psex2     bgrp2        tb       37.72592  37.87059   0.25679   0.01982           > summary(mypsex)          Call:     lm(formula = prtemp ~ psex + bgrp + tb - 1, data = apirt)          Residuals:         Min      1Q  Median      3Q     Max      -5.6877 -0.2454  0.0768  0.3916  1.6561           Coefficients:            Estimate Std. Error t value Pr(>|t|)         psex1 37.725922   0.135486 278.449  < 2e-16 ***     psex2 37.870591   0.135908 278.649  < 2e-16 ***     bgrp2  0.256794   0.066167   3.881 0.000115 ***     tb     0.019818   0.009342   2.121 0.034287 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.6888 on 621 degrees of freedom       (5 observations deleted due to missingness)     Multiple R-squared: 0.9997,     Adjusted R-squared: 0.9997      F-statistic: 4.788e+05 on 4 and 621 DF,  p-value: < 2.2e-16       As at now, I realise that `psex1` has the same coefficient value as `bgrp1`. **My questions are** :   1. Am I doing the right thing (in terms of the package and models used)?   2. Am I interpreting the output correctly?   3. Is there a correct way of obtaining the values of these factor levels in R?


title: Calculate ratings

text: Say I have a dataset that contains Users & all the sites that each user visited & how often each user visited it. For example:               User1 -> (Site1, 5 times), (Site5, 3 times)     User2 -> (Site2, 10 times) & so on..       Now based on how may times a user visited each site, I need to come up with a rating in the scale of 1 thru 5. For example: User1's ratings:               Site1 ->  5 (rating)  very interested     Site2 ->  1 (rating)  least interested     & so on...      What's the best statistical technique to compute these ratings? One easy way I can think of is to calculate the range (max - min). Divide the range in 5 buckets & see which bucket each Site falls under. Is this a good way to calculate ratings? Thanks for your help.


title: Transformation by subtraction (for continuous data arbitrarily limited in the experiment)

text: Consider a study looking at the relationship of chili-pepper plant height to number of chilies produced. Data was collected on all plants taller than 10cm. Initial look at the distribution of heights shows that most plants are on the short end, with a few reaching a max height of 150 cm as seen in a histogram. Taking the `log(plantheight)` normalizes the distribution somewhat. However, looking at `log(plantheight - 10)` gives a very nice looking normal distribution to the histogram.   * First, is this a valid transformation? Shifting the values to zero, then taking the log?   * Second, if not, what is the best way to deal with data that had been arbitrarily cut off in the experiment (no values below 10)?


title: Intervention Analysis - Pulse over several periods

text: I have a couple weekly time series and an intervention occurred over several weeks and then for some, after a period of no intervention, began again. So, the pattern is off for a period of weeks, then on for a period of weeks, then off for a period of weeks and perhaps on again. The image below tries to show the possible patterns: ![enter image description here](http://i.stack.imgur.com/DMYJg.jpg) When there is a single pulse, then as discussed here it is pretty simple to code the intervention as a pulse (1 or that period and zero elsewhere) and code the transfer function as one of several shapes - e.g. only effecting a single period (the pulse), a permanent level shift (a step) or perhaps with a lasting (but decaying) effect (a pulse but with decay transfer function). In the situation described here though, is there the same flexibility or do we really only have the option to code the interventions as 10 separate interventions (for the top series) or 5 separate ones for the bottom series and only allow the intervention to effect the single period? I am not sure how to have any of the "pulses" have a decaying effect because where does one start and the other stop? In my reading thus far for this type of analysis I have seen the pulse and step but nothing like this, where there is an "on" for a period of time, then "off" and then on again (maybe).


title: How to combine 2 variables each be strongly correlated with a 3rd variable?

text: I have 2 sets of variables that are weakly correlated to each other but highly correlated with third variable. Is there any method to combine these 2 variables to achieve a much more stronger correlation to third variable? for example: set a (correlation with b:0.2, correlation with c:0.8) set b (correlation with a:0.3, correlation with c:0.75) combination set: m*a + n*b (correlation with c:0.95)??(is it any mathematical approach to make a combination set and find weights?)


title: Which is better ROC curve or Confusion matrix?

text: anyone know about the predictive model evaluation? I'm confused about the ROC curve and confusion matrix. The area under the curve for ROC is represent about the accuracy of the classifier. But what about the confusion matrix, we can also find the accuracy. So which method is better for me to find the accuracy of the model?


title: why is adaboost predicting probabilities with so little standard deviation?

text: I'm using several algorithms to predict a binary target. So far I tried Gradient Boosting, Random Forest, Extra Random Trees and adaboost from scikit learn. All of these algorithms appear to predict probabilities ranging from close to zero to close to 1 with a very similar standard deviation. adaboost is the only one whose predictions are mostly compressed in the 0.4 to 0.6 range with only a minority falling outside of that range. This is not the first time I notice this behavior from this algorithm. Why is that the case? Secondly, if I wanted to blend these models (i.e average the probabilities), how would I account for the fact that adaboost probabilities's standard deviation is so different from any other algorithm? Should I rescale all predicted values from each algorithm to have similar mean and variation?


title: shuffle my data to investigate differences between 3 groups

text: I am sure this already exists but I just don't know the terminology to look for. I have three sets of 10 measurements. Each set corresponds to a different geographic region. So in total I have 30 measurements of my variable, and I have the factor "region" with 3 levels (west region, middle region, east region). Let's say I do a simple ANOVA and I get differences between the 3 regions. But I want to play a little with the possibility of this differences being "by chance". Or, in another scenario, let's say I can't use ANOVA for some reason (eg. strongly inhomogeneous variances) and I use a non-parametric test and I don't find differences I want to know if it's possible to do the following (or if the idea is appropriate): If there is really no difference between the 3 regions, then I can assume that any test (eg ANOVA or a non-parametric equivalent) will find approximate the same results even if I randomly mix all data once and again. So I thought I could simulate this, using my own data but just in different grouping. for example: 1- take all the 30 values from my own measurements 2- shuffle them into 3 groups, ie. randomly choose 10 values and assign them to a randomly chosen group; repeat with the next 10 data and then you have again 3 groups of 10 measurements. 3- run the test (eg. ANOVA) Now I go back to 1, and repeat this eg 1000 times, and see if there is a convergence towards a "stable" pattern. If there is, then there are actually no differences. If the convergence deviates a lot from the results I found with my "real" dataset, then I may think there are actually differences between the 3 regions. Is my reasoning correct/sound? I know there is something like this, I just don't remember the name.. I thought it was related to permutations but I'm not sure...


title: Comparing reaction tests

text: I have an experiment where every person needs to take reaction tests when he/she is sober, then when lightly and then heavily intoxicated. I need to perform a test that tells if the intoxication had effect on the reaction tests on a specific person. (I am NOT testing if the intoxication have effect on the reaction test in general). The reaction test protocol contains 60 hits in each state. I think these numbers do not follow a normal distribution, so I was planning to use Kruskal- Wallis non-parametric ANOVA. However, two of the assumptions are these:   * all individuals must be selected at random from the population   * all individuals must have equal chance of being selected Now this is not something that would be true in my case, right? I am not measuring 1000 hits and select 60 of them for each reaction test, but this 60 hits are all my population. Can I use Kruskal-Wallis?


title: non-collapsed data - repeated measures mixed ANOVA

text: I have to analyze results from an experiment and I’m not sure about how could I do it… The experiment is about the effect of some particularly distorted image typology on response time. The design is as follows: 30 participants watched two movies (A/B), each of them in a different image quality (distorted/no-distorted). Relationship between movie and IQ, order of presentation of IQ, and order of presentation of movie are counter-balanced between participants. Therefore, I understand it is a mixed design with 2 within-subjects IV (movie and IQ) and 3 between-subjects VI (relationship movie-IQ; order of movies; order of IQ). Dependent variable is reaction time (RT) to beeps inserted on movies. There are 10 beeps in each movie. I think it is not a completely mixed design, as all the participants watched both movies AND both IQ, but both movies IN both IQ. This is because we are mainly interested on the effect of IQ on RT, and not in movie effect, and we cannot use the same movie twice within a participant to avoid a “repetition effect”. To analyze results I have calculated mean RT of each participant in each IQ condition, and conducted a repeated measures mixed ANOVA: 1 within- subjects factor (IQ) and 3 between-subjects factors (relationship movie-IQ, order of IQ, order of movie). Then I have conducted another similar ANOVA, but with Movie as within-subjects factor. Results are marginally significant, so I would like to carry out a more powerful analysis. I think that I could do that if I do not collapse means of 20 beeps in each IQ condition, and compare Beep 1 in no-distorted condition to beep 1 in distorted condition, and so for each beep, and finally compare the differences between the two groups of beeps (the 2 IQ levels). But I am not sure about how could I do that in SPSS: I guess I have to create a column for each beep (1, 2,…10) in each condition, but then I don’t know how to create the both groups of beeps (one for each IQ condition), that is my actual interest. I hope someone could help me! Thank you!


title: Comparing coefficients in multilevel models

text: Is it meaningful to compare the coefficients of two different predictors in multilevel model when the two are at different levels? Specifically I have two variables which measure the same construct at two different levels (individual and group). I am using R package `nlme` to fit a multilevel model and want to test the hypothesis that the individual level predictor is significantly better than the group level predictor. Does it make sense to fit the model and then examine whether the confidence intervals of the individual predictor do not overlap with the coefficient of the group level predictor? Is the fact that they have different degrees of freedom a problem? Is there a better way to test this hypothesis?


title: Finding the generator matrix for a Markov jump process

text: Vehicles in a certain country are required to be assessed every year for road- worthiness. At one vehicle assessment center, drivers wait for an average of 15 minutes before the road-worthiness assessment of their vehicle commences. The assessment takes on average 20 minutes to complete. Following the assessment, 80% of vehicles are passed as road-worthy allowing the driver to drive home. A further 15% of vehicles are categorized as a “minor fail”; these vehicles require on average 30 minutes of repair work before the driver is allowed to drive home. The remaining 5% of vehicles are categorized as a “significant fail”; these vehicles require on average three hours of repair work before the driver can go home. A continuous-time Markov model is to be used to model the operation of the vehicle assessment centre, with states W (waiting for assessment), A (assessment taking place), M (minor repair taking place), S (significant repair taking place) and H (travelling home). generator matrix i got: $\space\space\space\space\space\space\space $ $\begin{matrix} W & \space\space A & M & \space\space S &\space\space\space H \end{matrix}$ $\begin{matrix} W \\ A \\ M \\ S \\ H \end{matrix} \left[ \begin{matrix} \frac{-1}{15} & \frac{1}{15} & 0 & 0 & 0 \\ 0 & \frac{-1}{20} & & & \frac{1}{25} \\ 0& 0 & \frac{-1}{30} & 0 & \frac{1}{30} \\ 0 & 0 & 0 & \frac{-1}{180} & \frac{1}{180} \\ 0 & 0 & 0 & 0 & 0 \end{matrix} \right]$ But I don't know how to find the 2 empty spaces in second row. Any help would be appreciated. W, A, M, S, H are labels of rows and columns.


title: What is the exact definition of profile likelihood?

text: Does anyone here know the exact definition of Profile Likelihood? Or does it have one?


title: Modeling time: Probability distribution over time?

text: I'm trying to model users' posting behavior during a day. Say we have a bunch of users, with the time they post tweets. Now, for each user, I would like to estimate the likelihood of he post a new tweet at 9:00am according to his historical posting behaviors. I'm curious what distribution I could pose here. In the literature I saw people using Gaussian, but I'm not sure if that's suitable since it's single peak. (Mixture model would be too complex for this task) Thus I'm wondering is there any distribution over time that I could use? The data I'm having is 2 month worth of tweets. Each tweet contains a timestamp, and the author's id. What I'm trying to model is the user's daily activity in terms of posting a tweet. E.G. the data for a single user would looks like [9:00, 9:12, 17:00, 17:01, 22:22, 22:37, 22:45, 22:47, 22:48...]. So this user post more tweets during night (around 10:30 pm) but very rare during work hours. I wish to model the probability P(user post a tweet|time). Would really appreciate the answers!


title: Anomaly detection: what algorithm to use?

text: Context: I'm developing a system that analyzes clinical data to filter out implausible data that might be typos. **What I did so far:** To quantify the plausibility, my attempt so far was to normalize the data and then compute a plausibility value for point p based on its distance to known data points in set D (= the training set): $$\text{plausibility}(p)=\sum_{q\in D}\text{Gauss}(\text{distance}(p,q))$$ With that quantification, I can then select a threshold that separates the plausible data from the implausible data. I'm using python/numpy. **My problems:**   1. This algorithm cannot detect independent dimensions. Ideally, I could put anything I know about the record into the algorithm and let it find out by itself that dimension X does not influence the plausibility of the record.   2. The algorithm doesn't really work for discrete values like booleans or select inputs. They could be mapped on continuous values, but it is counter-intuitive that Select 1 is closer to Select 2 than to Select 3. **Question:** What sort of algorithms should I look into for this task? There seems to be a ton of options including nearest neighbour based, clustering based and statistical approaches. Also, I have trouble finding papers that deal with anomaly detection of this complexity. Any advice is highly appreciated. **[Edit] Example:** Suppose the data consisted of the Height of a Person, Weight of a Person and Timestamp - so it's 3D-Data. Weight and Height are correlated, but the timestamp is completely independent. If I just consider the euclidean distances, I would have to choose a small threshold to fit most of my cross validation data. Ideally, the algorithm would just ignore the timestamp dimension, because it is irrelevant to determine whether a record is plausible, because the timestamp does not correlate with the other dimensions in any way. Any timestamp is plausible. On the other hand, one could make up examples where the timestamp does matter. For example it could be that value Y for feature X is plausible when measured before a certain date, but not after a certain date.


title: Is a constrained optimization problem equalivant to its Lagrangian form?

text: For the following problem: $\text{min:}\ f(x)\\ s.t. \ g(x)\leq t$ Is the above problem equivalent to the following problem? $\text{min:}\ f(x) + \lambda g(x) \\ s.t. \ \lambda\geq0$ where $t$ and $\lambda$ are variables. It seems equivalent, because if we increase $\lambda$, $t$ will tend to decrease, if we decrease $\lambda$, then $t$ will tend to increase.