title: Algorithm for rating books: Relative perception

text: So I am developing this application for rating books (think like IMDB for books) using relational database. **Problem statement :** Let's say book " **A** " deserves 8.5 in absolute sense. In case if A is the best book I have ever seen, I'll most probably rate it > 9.5 whereas for someone else, it might be just an average book, so he/they will rate it less (say around 8). Let's assume 4 such guys rate it 8. If there are 10 guys who are like me (who haven't ever read great literature) and they all rate it 9.5-10. This will effectively make it's cumulative rating greater than 9 (9.5*10 + 8*4) / 14 = 9.1 whereas we needed the result to be 8.5 ... How can I take care of(normalize) this bias due to incorrect perception of individuals. **MyProposedSolution :** Here's one of the ways how I think it could be solved. We can have a variable **Lit_coefficient** which tells us how much knowledge a user has about literature. If I rate " **A** "(the book) 9.5 and person " **X** " rates it 8, then he must have read books much better than " **A** " and thus his Lit_coefficient should be higher. And then we can normalize the ratings according to the Lit_coefficient of user. Could there be a better algorithm/solution for the same?


title: Clustered standard errors and robust standard errors

text: I was wondering if, when running a regression on panel data, clustered standard errors are already correcting for heteroskedasticity. Actually, I have run such a regression and detected heteroskedasticity. Since I used the pooled OLS model I have to cluster the standard errors anyway. Hence, I was hoping that I can address both issues simultaneously. Is that right?


title: Distribution for semi open interval

text: I have an application where I get back the time to perform an operation, the times are in the interval `]0:inf[` and I have this crazy idea. Lets say I know the mean of my dataset. woule it makes sense to test for normal distribution on the interval `]ln_n(0):mean[` if below mean and `[mean:inf]` above the mean, and if such a distribution existed to calculate std-dev's on it? If this makes sense to anyone but me I'm sure someone has allready done it and could tell me what its called.


title: How do I prove the significance of a non-linear model?

text: For example, for linear models, you take the p-value of the regression, and you deduce if the regression is or not significative. But with non linear models, in R, there isn't shown a p-value associated to a F-statistic, that say if the regression is or not significative. How can I know if a non linear model is significative? What tools can I use in R for that?


title: Bias for kernel density estimator (periodic case)

text: Kernel density estimator is given by $$\hat{f}(x,h)=\frac{1}{nh}\sum_{i=1}^{n}K(\frac{x-X_{i}}{h})$$ where $X_1,...X_n$ i.i.d with some unknown density $f$, $h$ - bandwith, $K$ - kernel function ( $\int_{-\infty}^{\infty}K(x)dx=1$, $\int_{-\infty}^{\infty}K(x)xdx=0$, $\int_{-\infty}^{\infty}K(x)x^2dx<\infty$). The bias can be computed using Taylor expansion: $$\int_{-\infty}^{\infty}\frac{1}{h}K(\frac{x-y}{h})f(y)dy-f(x)=\int_{-\infty}^{\infty}K(y)\left(f(x-hy)-f(x)\right)dy$$ $$=\int_{-\infty}^{\infty}K(y)\left(f'(x)hy+\frac{1}{2}f''(x)(hy)^{2}+o(h^{2})\right)dy=\frac{1}{2}f''(x)h^{2}+o(h^{2})$$ How to deal with periodic kernel and $f$ ($\int_{0}^{1}K(x)dx=1$,$\int_{0}^{1}K(x)xdx=0$,$\int_{0}^{1}K(x)x^2dx<\infty$)? How can I use taylor expansion?($\int_{0}^{1}\frac{1}{h}K(\frac{y-x}{h})f(y)dy=\int_{-\frac{x}{h}}^{1-\frac{x}{h}}K(y)f(x-yh)dy\neq\int_{0}^{1}K(y)f(x-yh)dy$ -I can't use kernel properties) Could you recommend a good book about kernel smoothing for circular data?


title: Quantile-Quantile Plot with Unknown Distribution?

text: ![N=2762](http://i.stack.imgur.com/bBDqg.png) **N=2762** I've been exploring a data set that seems to give rise to this kind of plot rather frequently. Would you say this is one population with a different than normal population? Or are two populations confounding the normal distribution? It used matplotlib and scipy.stats:                     (x,y), (slope, inter, cor) = stats.probplot(data, dist='norm')                osmf = x.take([0, -1])  # endpoints           osrf = slope * osmf + inter           self.ax.plot(x, y, ',', osmf, osrf, '-', **self.kwargs)


title: What are the advantages / disadvantages of using splines, smoothed splines, and gaussian process emulators?

text: I am interested in learning (and implementing) an alternative to polynomial interpolation. However, I am having trouble finding a good description of how these methods work, how they relate, and how they compare. I would appreciate your input on the pros/cons/conditions under which these methods or alternatives would be useful, but some good references to texts, slides, or podcasts would be sufficient.


title: Central Limit Theorem Tails

text: I'm trying to compute some p-values for samples from a distribution of sums of ~1000 random variables. The exact distribution of these random variables isn't known, but I have empirical estimates that I think are pretty accurate. So far I've been using the central limit theorem to produce a normal approximation for this sum, which does OK since n is relatively large, but not great. Computing the exact distribution of the sum via convolutions is too slow, but I only really care about having high accuracy in the tails; everywhere else a rough approximation is fine. Are there any methods that will allow me to improve my estimates in the tails (or just one tail) without having to compute the entire convolutions? I'm not sure if this would be a variant on the CLT implementation or something completely different. My gut feeling is that there isn't any way to do this, so I'd be open to any kind of solution at all!


title: Intraclass Correlation Coefficient Question

text: I have two datasets in which I probably need to use the intraclass correlation coefficient, and I want to make sure I choose the right type. **Dataset 1** : I had 20 randomly selected human subjects rate 50 faces based on their attractiveness on an ordinal scale from 1-5. Each subject rated the same faces. I want to know how similar their ratings are to each other. Which method should I use? 1-k? **Dataset 2** : I have 20 4x4 confusion matrices, one for each human subject. I've been asked by a reviewer to measure how similar the confusion matrices are to each other. ICC seems like a good tool to use here, and I could do this by treating each confusion matrix as a 16-dim vector and then feeding the 16x20 array of data into an ICC algorithm. Again, which ICC method is most appropriate? 1-k? In my niche, I don't see people use ICC much, so I don't have any experience with the method and wanted to double check with experts to prevent myself from doing something stupid. Here are the types in my analysis software: '1-1': The degree of absolute agreement among measurements made on randomly selected objects. It estimates the correlation of any two measurements. '1-k': The degree of absolute agreement of measurements that are averages of k independent measurements on randomly selected objects. 'C-1': The degree of consistency among measurements. Also known as norm- referenced reliability and as Winer's adjustment for anchor points. 'C-k': The degree of consistency for measurements that are averages of k independent measurements on randomly selected objects. Known as Cronbach's alpha in psychometrics. The degree of consistency for averages of k independent measures made under the fixed levels of column factor. 'A-1': The degree of absolute agreement among measurements. Also known as criterion-referenced reliability. The absolute agreement of measurements made under the fixed levels of the column factor. 'A-k': The degree of absolute agreement for measurements that are averages of k independent measurements on randomly selected objects. The degree of absolute agreement for measurements that are based on k independent measurements made under the fixed levels of the column factor. Thanks!!


title: Is it better to use data imputation for missing data or an analysis that is not affected by missing data (e.g., HLM/mixed effects modelling)?

text: I have treated two groups of 100 people with different treatments. I have pre- treatment and post-treatment data for most participants (as well as 1-month follow-up. I also have weekly data for some variables, but may or may not include those). About 10% of cases dropped out before post-treatment in one group and 30% of the other group dropped out at follow-up. I want to compare the efficacy of the treatments in each group (i.e., how group membership is associated with symptoms over time [continuous]). Note that the % of drop-outs is not part of the efficacy analysis (this will be considered in another analysis).   * I have been told to use multiple imputation to help with the missing data problem. This option would involve doing a repeated-measures ANOVA following the data imputation and an attempt to pool the results from the repeated analyses using the different simulated databases.    * I have also been told to use hierarchical linear modelling (i.e., mixed effects modelling), since this analysis is quite robust to missing data. Which option is better? Should I use a method that simulates missing values (e.g., multiple imputation) or use what data I have with a method that is robust to missing data? This question fascinates me and also has important implications for my line of work! Any guidance would be appreciated. Lets assume I can do either analysis (mostly) properly. Edit: Also I should note that another option would be to use BOTH methods (i.e., imputation followed by HLM). I appreciate this is an option. If this is your recommendation, please go ahead and explain why, but please also explain which of the two options is better if they are mutually exclusive.


title: Interpreting the coefficients of ARCH Lagrange Multiplier Test

text: I am new to econometrics and I am building my first econometric model. I ran the LM test on a univariate time series data of 12000 observations and got the following stats:               Chi-squared = 458.5313, df = 12, p-value < 2.2e-16       Can I reject the null? Is there presence of ARCH effects?


title: Predicted probabilities for multinomial logistic regression

text: I have three variables, a factor (`c`) as the dependent variable and two ordinal independent variables (`a, b`). Each variable has five categories (`1,2,3,4,5`). Thus, I fitted a multinomial logistic regression (`testus`, see below) with the `car` package. Now I would like to get the predicted probabilities.               > a <- sample(5,100,TRUE)     > b <- sample(5,100,TRUE)     > c <- sample(5,100,TRUE)     > required(car)     > a <- as.ordered(a)     > b <- as.ordered(b)     > c <- as.factor(c)     > testus        <- multinom(c~a+b)     > predictors    <- expand.grid(b=c("1","2","3","4","5"), a=c("1","2","3","4","5"))     > p.fit         <- predict(testus, predictors, type='probs')     > probabilities <- data.frame(predictors,p.fit)      Now I got the predicted probabilities for a under b and c.               > head(probabilities)       b a         X1         X2         X3         X4        X5     1 1 1 0.10609054 0.22599152 0.20107167 0.21953158 0.2473147     2 2 1 0.20886614 0.27207108 0.08613633 0.18276394 0.2501625     3 3 1 0.17041268 0.24995975 0.16234240 0.13111518 0.2861700     4 4 1 0.23704078 0.21179521 0.08493274 0.03135092 0.4348804     5 5 1 0.09494071 0.09659144 0.24162612 0.21812449 0.3487172     6 1 2 0.14059489 0.17793438 0.29272452 0.26104833 0.1276979`      The first two columns show the categories of the independent variables `a` and `b`. The next five columns show the conditional probabilities (e.g., $P(c=1|b=1~\&~a=1) = 0.10609$. But now I would like to know only the predicted probabilities for `c` under `a` or the predicted probabilities for `c` under `b`. Is this possible?


title: What is the minimum mean squared error estimator of the mean of two normally-distributed variables?

text: Let $Z=(X+Y)/2$, where $X$ and $Y$ are independent normally-distributed random variables with known variances $\sigma^2_X$ and $\sigma^2_Y$ and unknown (and possibly different) means. Given a sample $x_1$ from $X$ and $y_1$ from $Y$, what is the minimum mean squared error estimator of the mean of $Z$? Is there a biased estimator with an MSE improved over the maximum likelihood estimator $\frac{x_1+y_1}{2}$? Can we generalize to $n$ mutually-independent random variables?


title: Intractable posterior distributions

text: In Bayesian statistics, it is often mentioned that the posterior distribution is intractable and thus approximate inference must be applied. What are the factors that cause this intractability?


title: Data transformation to fit gamma distribution in R

text: I'm having trouble to fit and simulate a `gamma distribution` using the `fitdistr` function from the `MASS` package: My data is `daily rainfall`, so i 'm adding `1` to all points to get rid the 0's:               require(MASS)     rainfall = scan('daily_rainfall_january.txt', '')          rainfall = rainfall + 1     # sometimes i get warnings here (NA produced)     fit = fitdistr(rainfall, 'gamma')      Now i want so simulate a January rainfall: (plot shows the accumulated rainfall over N simulations) ![enter image description here](http://i.stack.imgur.com/XDT4B.png) `1 - Blue and red: output from 1000 rgamma simulations` `2 - Dark line: 30 year climatology mean (monthly accumulated rainfall)`               shape = fit$estimate['shape']     rate = fit$estimate['rate']     rainfall.simul = rgamma(31, shape = shape, rate = rate)      Now, how do i transform back `rainfall.simul`? Is there a better way to handle the 0's problem? Thanks!


title: Can Chauvenet's criterion be used with non-normal data?

text: Can I use Chauvenet's criterion on set of observations where a normal distribution cannot be assumed?


title: How to combine annual rainfall data with specific forecast of known accuracy to predict weather for a particular day?

text: > Alice is getting married tomorrow at an outdoor ceremony in a beautiful > garden. > > In recent years, it has rained, on average, only 5 days in a year. Alice is > anxious about the weather for tomorrow and wishes to get the weather > prediction from a weather forecaster. > > From the past records, the weather forecaster is correct 90% of the time. > Unfortunately, the weather forecaster predicts rain for tomorrow. > > Should alice make wet weather plans for tomorrow? Found, _e.g._ , at http://stattrek.com/probability/bayes-theorem.aspx.


title: A simple classical statistics problem that I have trouble with

text: This question was from my CAs. Nothing has been changed. A local weatherman forecasts that there is a 65% chance of rain on the coming Saturday and a 65% chance of rain on the Sunday immediately following that Saturday. What can you conclude from the above weather forecast?   Choose one of the following options. A. It will certainly rain during the weekend (that is, it will rain on at least one of the two weekend days). B. The chance of rain on Saturday and Sunday is 42.25%. C. The weather forecast does not sound right, as P(rain on Saturday)+ P(rain on Sunday) > 100%. D. The chance of rain on Saturday and Sunday is greater than 50%, but less than 85%. E. The chance of rain on Saturday and Sunday is between 30% and 65%. F. None of the above. Chance of raining on Saturday = 0.65 Sunday = 0.65 Therefore, chance of raining on Saturday and Sunday is 0.65*0.65 = 0.4225 Doesn't this means there are two answers, B and E? Which is more correct?


title: Rule of thumb reliability of scales for using manifest variables in path analysis?

text: I am using path analysis to test hypotheses in my study. I have been told by a few people in passing that I can use manifest variables rather than latent variables as long as my reliabilities are high. The Cronbach's alpha reliabilities of my constructs are 0.81, 0.88, 0.87, 0.90 and 0.92.   I know in general, above 0.7 is considered acceptable reliability for scales, but I don't know if that's different for path analysis. Coffman and MacCallum (2005) write that low reliabilities bias the parameter estimates, and that those parameters will in general be lower than if reliabilities were perfect. However, they didn't provide a cut-off... I.e. what level of reliability is generally deemed acceptable before necessitating the use of latent variables? Should I assume that 0.7 is the cutoff? And if so, is there a reference I can cite for this that specifically relates to path analysis? Thank you in advance. EDIT: In my analyses I have used manifest variables. What I really want to know is whether I am justified in using manifest variables, or whether it's necessary for me to use latent variables.


title: Should two equal pseudo-random numbers in a sample count as a "run down" or a "run up" in a runs test?

text: Taking into consideration the runs test proposed by Knuth given a sample of pseudo-random numbers to test independence, and looking at this example:               X = (5,4,1,7,2,3,6), which yields     S = (0,0,1,0,1,1)      I'm a bit confused if:               X = (5,4,1,**1**,2,3,6), would yield      S = (0,0,1,**1**,1,1) ... call it option #1 or yield     S = (0,0,1,**0**,1,1) ... call it option #2      ...and for the case of:               X = (5,**5**,1,7,2,3,6), should it yield     S = (**0**,0,1,0,1,1) ...call it option #1 or     S = (**1**,0,1,0,1,1) ...call it option #2      So what is it in the first and second cases and why?


title: How to report Cronbach's alpha for referent-shift scale (before or after aggregation?)

text: I have a study variable that is measured at the individual level and then aggregated and to the group level. My hypothesis testing is conducted at the group level. I used an 8-item referent-shift measure of team psychological capital. When calculating and reporting the coefficient alpha, should I be looking at it across individual level scores, or averaging (aggregating) each item across each team, and then calculating the coefficient alphas across aggregated item scores. That is, I can find the average value of item 1 within team 1, then the average value of item 2 within team 1 etc. I will then have one value for each scale item per team. I will then calculate the alpha based on n=94 teams. Alternatively, I could leave everything at the individual level, and calculate it that way across all individuals within my sample (n=456 individual team members). I would think I should measure it at the group level, since I am conducting my analyses as the group level. That is, I should be most concerned about the level of measurement error in the aggregated team-level scale, right? If I calculate it at the team level, I get alpha=0.92, and if I calculate it at the individual level, I get alpha=0.86. At the end of the day, they're both good alphas, but I still want to report the correct number. Also, in case this is important in answering the question,   I have calculated ICC(1), ICC(2) and rwg for the scale means (rather than on an item by item basis) and they are acceptable for aggregating the scale **means** to the group level. Thank you in advance! Note: Referent-shift measure means that the referent of the scale is the team rather than the individual. So the individual-referent psychological capital item might be **I** am optimistic about my work. The team-referent psychological capital item would be **Members of my team** are optimistic about their work. See Functional relations among constructs in the same content domain at different levels of analysis: A typology of composition models. Chan, David Journal of Applied Psychology, Vol 83(2), Apr 1998, 234-246. Edit: While the scale is measured at the individual level, analyses are conducted at the group level (after aggregation). Hence my confusion over whether I should report individual level alpha or group level alpha... Group level alpha would be found by aggregating each individual item (find average within-group value for each item and then calculate alpha based on that).


title: Denominator of dependent variable changes from year-to-year

text: I have a model with a ratio (percentage) as dependent variable. The denominator of this DV also changes from year-to-year. I remember that we should be more careful with ratios as DV but I'm a little rusty on the specific steps that should be taken in this situation. Some suggest to transfer the ration into a logic (Papke, L. E., & Wooldridge, J. M. 1996. Econometric methods for fractional response variables with an application to 401(k) plan participation rates. Journal of Applied Econometrics, 11(6): 619–632.). Would this also apply if the ratio exceeds 1? Thanks! Ben


title: Pre/Post analysis without control

text: Quite a newbie to this, any info or pointers are appreciated, but please be as precise as possible. If I have a feature in general that was launched mainly to improve the user experience, but interested to understand if it also helps to boost the movie sales on the website pre and post launch (the feature was launched to the entire site, hence no control group was designed to facilitate statistical testing), what type of analysis and what techniques I can use to account for   1. Movie sales usually decays after release date, e.x. 1st week crazy, 2nd week ok, 3rd no sales (a bit extreme). Down in sales count in post launch does not mean the feature affected the sales. Conversion rate might be more smooth, but it turned out it also has similar trend, but just not so rigid. (1st week tends to have higher conversion rate)    2. Movies are not independent, related movies tend to sale together, new movie will boost sales for old related movies. e.x. Transformers 4 boosts sales for transformers 1~3. So an increase of sales in transformers 1~3 does not mean it's brought up by the feature. How do I isolate those affects so that I know the difference is caused by the feature alone assume that was the only feature that was launched and no other features? Presumably average movie decaying rate over time is available, but just don't know how to use that to weight the sales if possible?


title: Can decision tree be used for fraud detection in this way?

text: A large dataset with more than 100 variables including a target variable. A small portion of target = 1 cases are fraud or due to other errors. I want to identify these target = 1 cases, i.e. fraud or error. Assuming most cases are good, I intend to use decision tree to classify cases. In one leaf with a high percent of target = 0, e.g. in 95% cases target = 0, and in 5% cases target = 1, then I think the 5% cases with target = 1 in this leaf are classification errors that include random error and the error caused by fraudulent behavior, so they can be consider as suspected fraud cases. Does this make sense? As for the decision tree leaves with close percentages of target = 1 and target = 0, e.g. 43% vs. 57%, what can I do? Based on @Zhubarb's suggestion, I made the following changes to the model and its interpretation:   1. Set leaf size, e.g. at least 100 observations in a leaf   2. My goal is to identify fraud at the organization level instead of the individual observation level. Considering the prior probability problem as Zhubarb mentioned, I interpret the decision tree model in a new way: take a leaf with 30% target = 1 and 70% target = 0 for example, and assume the expected positive percentage should be 30% for the individuals in this leaf. Then look at each organization; if in this leaf one organization's positive percentage is much higher than 30%, it can be considered as a suspect organization. In this way, all leaves regardless of whatever positive percentages are high or low can be used. Any comments are appreciated!


title: repeated measure ANOVA or GLMM?

text: I would like to ask for some advice on an analysis that I am doing at the moment. The experiment: I ran an experiment with animals in which I subjected 15 pairs of individuals (one male and one female in each pair) to three olfactory tests: treatment 1, treatment 2, treatment 3. Each test procedure involved a setup with two scent locations; one for the treatment and one for a control scent. Each pair was only subjected to one test per day, with the three tests done on different days. The order of the tests was random. I recorded 7 different response variables per individual for each treatment and each control scent on each test day. Initial analysis showed that other than differences in response to the different treatments, there are also differences in response for individuals, pairs, between sex and also depending on the number of the trial – as the number of trials that the animals are exposed to increased, the overall response decreases (some habituation?). Now I want to test for differences in response between the treatments for each response variable. I am thinking repeated measures ANOVA might be suitable, however, this does not account for the added variation due to the differences in things like sex, and trial number. If I use generalised mixed models, things like trial number can be added as random variables. However, with treatment being a categorical variable, it does not give very clear results as to (significant) differences between means of responses to the different treatments. What would be the best approach here?


title: If $X$~$Exp(1)$ and $Y$~$Exp(2)$, what is the pdf of $Z=X+Y$?

text: I am trying to find the following: If $X\sim\mathrm{Exp}(1)$ and $Y\sim\mathrm{Exp}(2)$, what is the pdf of $Z=X+Y$? I tried to use the convolution formula but am not sure what the limits of the integral are: $$g(z) = \int_{0}^{z} 1 \cdot 2\cdot\exp(-1\cdot(y-x))\cdot\exp(-2x) \mathrm dx.$$ I ended up getting: $g(z) = 2\cdot(1-e^{-z})\cdot e^{-y}$. Would this be correct? Thanks!


title: Does filtering of data by effect size violate some assumption of P value adjustment methods?

text: I have pre- and post-treatment continuous data for a large number of variables that I am analyzing for treatment effect. Normally I would obtain the P values and then adjust them for multiple testing with a method such as Benajmini- Hochberg. However, the statistical test that I am using is somewhat computationally intensive. To reduce this load, I am thinking of first filtering the data by removing variables for which the average effect size is less than 2-fold in either direction; i.e., `mean(absolute(log2(post/pre))) <1`. Will such filtering of data by effect size violate some assumption of P value adjustment methods? It seems that the filtering is likely to enrich for variables for which a truly positive treatment effect exists.


title: An unbiased estimate of the median

text: Suppose we have a random variable $X$ supported on $[0,1]$ from which we can draw samples. How can we come up with an unbiased estimate of the median of $X$? We can, of course, generate some samples and take the sample median, but I understand this will not in general be unbiased. Note: this question is related, but not identical, to my last question, in which case $X$ could only be sampled approximately.


title: Philosophy of sampling for experiments: finite versus infinite

text: What is the difference between a finite population and an infinite one - when you are designing an experiment (sample/power and interpretation of the results)? Say a company has a database of 20,000 customers. Given that a response to some stimulus is relatively small, and a meaningful minimum detectable difference is also small, if you run a power analysis for a 2 sample proportion, you may find that you need 2 groups of 15,000 for the experiment. Do you quit and say you cant experiment on this population? Or do you (somehow) instead treat the population as finite and run a power analysis that way? What are the implications? **ADD:** What I would like to know, and wanted to add this detail in case the last part of the question wasn't completely clear - is what the difference is between   1. The inference with assuming an infinite population, say the inference from a logistic regression model with `glm()` in R.   2. The inference with assuming a finite population, say inference from a logistic regression model with svyglm in R? Will #1 allow inference about the wider population / data generating process and #2 only allow inference about that particular population (which is assumed fixed)?


title: How to compute distances within & between clusters, if clusters contain sequences?

text: I applied number of methods of clustering, and I want to evaluate these different methods using Dunn index, in this method I have to calculate the distance among clusters and among points in clusters. My question is: if my algorithm is clustering users depending on their sequences, i.e., each user has sequences, then the algorithm measures the similarity among sequences using seq. alignment technique, then cluster them. Sequence alignment tool considers two users most similar when they have max score. If my clusters are id's of users or sequences of users, how do I compute the distance among clusters and among sequences within a cluster?


title: Forward Filtering Backwards Sampling (FFBS) and Look-Ahead Bias

text: **Assumptions / Context:** Let's assume that I have data that can be modeled as a dynamic linear model. To estimate the parameters (e.g., covariance matrix of the state/system equation), I use a Gibbs sampling or other Markov Chain Monte Carlo procedure. **Problem:** The problem is that to use the Gibbs sampler I use the FFBS algorithm which uses the entire history of the data series. Thus, if I want to test the accuracy of my model's predictions in the past I would be using a model fitted with data from the future. Therefore, there is look-ahead bias in the performance measures that I may calculate (e.g., mean squared error of the residuals). **My Initial Though Surely Inferior Thoughts:** Assuming I have "point-in- time" data (e.g., I know when the data became available to modeler in the past), a "brute-force" way would be to estimate the model's parameters using data from period $1$ to $t-1$ and forecast the observation variable $y_{t}$ and repeat this procedure for all $t \in T$. So, if I had a univariate time series with 80 observations I would need to repeat this process 79 times (first for the data in period $t=1$, second for data from periods $t=1:2$, third for data from periods $t=1:3$, and so on). I would have 79 different models for this time series. Then, I would calculate some performance measure (e.g., _Theil's U_ ) using the residuals from my forecasts of $y_{t}$. **Questions:** 1) Am I incorrect about the look-ahead bias created by the FFBS algorithm?, and 2) Does anybody know a way to efficiently estimate a State Space time series model using Bayesian techniques?


title: Do you need to change instances to rates for OLS regressions

text: I am interested in performing a regression on data on a population. This population causes events X and Y. I have monthly data for population, event X, and event Y. Do I need to change my variables X and Y into rates (X divided by population) before running the regression? Some details about the model: the population are students at a university. That's one IV. Another IV will be a yet to be determined demographic characteristic. As output variables, I have cheating instances and plagiarism instances. I'm going to combine these into one DV "academic dishonesty" and run three models, one with cheating as a DV, the other with plagiarism, and the last with academic dishonesty as the DV. In all cases university population and demographics will be the two IVs. Should I change these cheating and plagiarism instances into rates? Thanks for the advice.


title: Controlling type I error in permutation tests

text: I'm implementing a permutation test for a test statistic with unknown distribution. Apart from my concrete situation (see below), I'd like to know first   1. If and under what circumstances permutation tests might inflate the type I error, that is, in what kind of situations they may produce anti-conservative results? In my concrete example, I'm trying a modification of the 'optimal cutpoint'-approach (see for instance "Practical p-value adjustment for optimally selected cutpoint" by S. Hilsenbeck and G. Clark from 1996 in Stat. Med. (15): 103-112). But instead of dichotomizing the whole continuous scale, I start with a fixed group at the lower end of the continuous scale and compare this group with a group from the higher end of the continuous scale (the data in the middle are ignored). In a systematic fashion, only the upper cutpoint is moved towards lower values of the continuous predictor searching for a maximized effect. In order to adjust for the maximum p-value selected from this search, the data is resampled in a permutation approach. However, in a simulation study with 1000 repetitions with each 500 permutations, the type I error seems to be inflated (i.e., under a postulated null effect, the proportion of significant test results exceeds the nominal alpha niveau). So my second, more concrete question is   1. Are the simulation results reasonable and what might be an explanation for the results? Are 500 permutations too few and not precise enough (however, the deviances from the nominal alpha niveau are almost exclusively above the nominal alpha niveau)? What could be done to fix this problem? I'd be glad for any kind of help!


title: Prejudice in blind test

text: There is a blind test that asks participants to distinguish Coke and Pepsi. A participant will test 6 cups of drink and tell whether it was Coke or Pepsi. Assuming that the participant can tell the difference between them, though not perfect, if he judged that the first three were all Coke, then he would think it is much more likely that the remaining three would be Pepsi rather than Coke. Then the last three experiments cannot give precise information about the participant's ability to distinguish. Such a problem can happen when the participant knows a priori that there would be equal numbers of Coke and Pepsi. Then the experimenter might think of flipping a coin six times to decide which cola is given to the participant. Then it is possible that all or most of 6 cups are Cokes. It may be problematic when the participant is good at telling Coke a Coke but not as good at telling Pepsi a Pepsi. Is there any methods to solve this problem?


title: The standard normal distribution vs the t-distribution

text: Given an IID normally distributed sample $X_1,...,X_n$ for $n$ small with mean $\mu$, standard deviation $\sigma$, sample mean $\overline{X}$ and sample standard deviation $s$ (the unbiased estimator form). I understand that $$\frac{\overline{X} - \mu}{\frac{\sigma}{\sqrt{n}}} \sim N(0,1),$$ but I'm having trouble reconciling this with the fact that $$\frac{\overline{X} - \mu}{\frac{s}{\sqrt{n}}} \sim t_{n-1}.$$ Since the $t$-distribution is like the standard normal distribution but with a higher variance (smaller peak and fatter tails), this would seem to suggest that the sample standard deviation $s$ systematically underestimates the population standard deviation, making it in fact a _biased_ estimator.


title: Good CFA fit indices, but poor fit indices at structural equation modeling

text: I am currently working on my research and I have the following problem. I have conducted a CFA on all model variables and after deleting the troublesome items, the CFA yielded good fit indices. However, when i conduct a structural equation modeling (with item-averaged composites of the resulting items), I have poor fit indices (e.g. RMSEA 0.145, CFI: 0.80 etc.). In my model I have two independent variables, two mediators and one dependent variable. I added the error terms on the dependent variables and the mediators in AMOS. In addition, I covaried the independent varuables with each other and I also covaried the control variables with the independent variables. Finally, I also fixed a path from the control variables to the dependent variable. What are the methods to improve the fit indices in such situations?


title: What is the difference between GARCH and ARMA

text: I am confused. I don't understand the difference a ARMA and a GARCH process.. to me there are the same no ? Here is the (G)ARCH process $\sigma_t^2 = \underbrace{ \underbrace{ \alpha_0 \+ \sum_{i=1}^q \alpha_ir_{t-i}^2} _{ARCH} \+ \sum_{i=1}^p\beta_i\sigma_{t-i}^2} _{GARCH}$ And here the ARMA $ X_t = c + \varepsilon_t + \sum_{i=1}^p \varphi_i X_{t-i} + \sum_{i=1}^q \theta_i \varepsilon_{t-i}.\,$ Is the ARMA simply an extension of the GARCH; GARCH being used only for returns and with the assumption $\sigma = r\epsilon$ where $\epsilon$ follows a strong white process ... Can someone please help me out ..? Thanks


title: How to evaluate $\int_{0}^{\infty} x^4 e^{-x^2/{\beta}^2}dx$?

text: Can someone show me how to evaluate $\int_{0}^{\infty} x^4 e^{-x^2/{\beta}^2}dx$? I am working on verifying that a function is a pdf, and finding its expectation and variance. I want to be able to do it step-by-step, without just getting the answer from a CAS. Thanks in advance!


title: Interpretation of PCA Report from Machine Learning Toolkit

text: I am using a Machine Learning package called Orange. This software provides a Python interface and a visual programming GUI. I am using the GUI at present. There were 48 original variables. I have been using the PCA widget for dimensionality reduction and I want to investigate the relationship of PCA components to the original variables. I have done a little reading but that seem to be confusing me. My naive understanding is that PCA transforms the data set by a linear combination of the original variables. Is it possible find out the contribution of the original variables to the PC? The PCA widget produces a report, a snippet is below (didn't bother to include all PC).                               Comp.1  Comp.2  Comp.3  Comp.4 ....     Std. deviation  0.593   0.388   0.296   0.238 ....     Proportion Var  35.160  15.021  8.739   5.672 ....     Cumulative Var  35.160  50.181  58.920  64.592 ...      So if PC1 has a Proportion of Var 35.16, then essentially this is 35.16% of each of the original variables + proportion of PC2 etc? The widget also provides EigenVector (2x48). Do I need to do some magic with the EigenVectors to assess the influence of the original variables? Finally, if some PCA components have proportion of variance of zero or close to zero does this mean I have some poor/bad/redundant features? If so, is there a simple way to determine which they are form the PCA results?


title: Random effect significance in linear mixed model

text: I have performed LRT:               > (lL1 <- logLik(mixedModel))          'log Lik.' -65435.99 (df=4)     > (lL2 <- logLik(linearModel<-lm(y~x)))          'log Lik.' -78234.86 (df=3)     > ((Delta <- as.numeric(lL2-lL1)))          [1] -12798.86     > pchisq(- 2 * Delta, 1, lower.tail = FALSE)          [1] 0      How can I understand that 0 from chi-square? Are random effects are statistically important?


title: Negative upper bound value in genetic model estimate

text: I'm running twin models in which I estimate additive genetic and environmental influences in a bivariate model, mixing continuous and ordinal measures. In my model I estimated an overlap between E estimates (non shared environmental influences), which brings out the following output:               ACE.Re[2,1]   3.020748e-01         6.575223e-01        -7.953737e-01      What I don't understand is how the upper value can be negative? Has anybody else come across this problem?


title: Why should random numbers be drawn from a distribution?

text: What does it mean to draw a number at random from a normal distribution? Why can't you just "draw" numbers at random? Why do you have to draw random numbers from a distribution?


title: How to impute or predict a characteristic when one of the IVs in the prediction is other household members having that same characteristic?

text: In some data set A, we have: household id, person id, age, sex, and then a simple binary `likes donuts / does not like donuts` variable. In some other data set B, we have just household id, person id, age, and sex. We have set up a regression in data set A with independent variables: age, sex, and `does anyone in the household like donuts?` in order to predict who likes donuts in data set B. People living in households where _anyone_ likes donuts are more likely to like donuts themselves, so we don't want to throw out that _intra_ -household probability when imputing donut likeage. However, I'm unclear how to implement something computationally that maintains these within-household donut affinities when we're imputing `do you like donuts` in data set B? Do we impute at the household-level first, and then impute a subset of the people within sampled households? I feel like I should take a page out of the cluster sampling playbook, but I'm not really sure where to start looking.


title: How to pick a trained model instance when its performance fluctuates a lot?

text: I have a dataset that I have to perform a regression task. I split the dataset into 80% for training and 20% for validation.               train.index <- sample(N.rows, N.rows * 4/5, F)     train.data <- data[train.index,]     val.data <- data[-train.index,]      I used the GBM in R which has the parameters specified (and other take defaults). However, this applies to any kind of algorithms with a grain of randomness embedded.               train.model <- function(train.data, val.data)     {         # 1) model <- gbm(formula = /* ... */, data = train.data, n.trees = 1000, interaction.depth = 3, n.minobsinnode = 5, shrinkage =  0.01, distribution = "gaussian")         # 2) predict via model on val.data         # 3) calculate error     }      Suppose I executed train.model with the same train.data and val.data every time, and I do it about 50 times. The error (Mean absolute error) percentile is about 5-10%, which is rather high. I could select the best model out of the 50 runs, but this assumes I have a validation set which I do not use as training set. My question is this: **Suppose I am using all the data for training (normal cross validation practice), how could I tell if the model is performing poorly and avoid these major fluctuations?** Thank you.


title: How to qualitatively describe distribution shapes

text: I have a distribution which looks a bit binomial, though with a cutoff, i.e.: $P(x) = 0$ if $x \le 0$, rises to a bell-shaped peak at around $x=0.2%$, appears to tail off exponentially until an abrubt cut such that $P(x)=0 $ if $x>1$. I don't want to make any claim as to what kind of distribution this is, as the mechanism behind it is too complex to analyze. Are there some adjectives that formally describe its shape though? A word that means "a bit binomial- looking"?


title: Selecting PCA components which separate groups

text: I frequently used to diagnose my multivariate data using PCA (omics data with hundreds of thousands of variables and dozens or hundreds of samples). The data often come from experiments with several categorical independent variables defining some groups, and I often have to go through a few components before I can find the ones that show a separation between the groups of interest. I have came up with a rather primitive way of finding such discriminating components, and I wonder   1. to what extent this is reasonable / justifiable, and   2. whether there are better ways of achieving the same. Note that this is exploratory. Before I convince anyone else, I want to convince myself. If I see that there are components that clearly distinguish the groups of interest (e.g. control vs treatment), even if they are responsible for a minor portion of the variance of the responses, I trust it more than a result from, say, supervised machine learning. Here is my approach. I will use the "metabo" example data set from pca3d in R. The idea is to assess how much variance of each of the component can be explained by the independent variable. For this, I calculate a simple model for each component and use $R^2$ as a metric to order the components from "most interesting" to "least interesting".               require( pca3d )     # data on metabolic profiles of TB patients and controls     data( metabo )     # first column is the independent variable     pca <- prcomp( metabo[,-1], scale.= T )           # create a model for each component     lm.m <- lm( pca$x ~ metabo[,1] )         lm.s <- summary( lm.m )         lm.r2 <- sapply( lm.s, function( x ) x$r.squared )     plot( lm.r2, type= "l" )     text( 1:length( lm.r2 ), lm.r2, 1:length( lm.r2 ), pos= 3 )      Here is the result. The plot shows the percentage in variance of each component explained by the independent variable in `metabo[,1]`. ![enter image description here](http://i.stack.imgur.com/RMb1L.png) We can sort the components by $r^2$ to find out which ones to display with `order( lm.r2, decreasing= TRUE )`; the first three components are 2, 1 and 7.               pca3d( pca, components= c( 1, 2, 7 ), group= metabo[,1] )      Here is the plot: ![enter image description here](http://i.stack.imgur.com/EWluB.png) (The red and green categories are two groups of subjects who are not patients, and it is to be expected that they cannot be distinguished.) To reformulate my questions,   1. Does this approach make sense to you? My problem is that it looks too much like data dredging. Also, intuitively I think maybe I should turn the table and ask what part of the variance in the independent variable is explained by each variable? Finally, I'm (almost) sure that I am reinventing the wheel, poorly, so my second question is   2. Is there anything better? Note that I do not want to switch to partial least squares or anything similar at this stage; I just want to diagnose the PCA in the context of my classification.


title: How to compare rates of occurence in consecutive time series count data?

text: My data consists of occurrences of words in time windows. E.g.:               Day; Word; Frequency     1; "dog"; 45     1; "cat"; 2     ...     2; "dog"; 90     2; "cat"; 4     ...      I would like to estimate the ratios of all day-to-day differences (i.e., for dog day 1->2: 90-45/45 = 100%). For cat the increase is also 100%, but due to the small sample size I would like to somehow quantify that it is "less trustworthy". Something similar (for binomial data) is proposed here: > http://www.evanmiller.org/how-not-to-sort-by-average-rating.html But with count data it's not quite the same... Any ideas are most welcome.


title: In importance sampling, why should the importance density have heavier tails?

text: Why should the importance density (or biasing density) $g$ have heavier tails than the original distribution $f$? Equivalently, why should $$\frac{f(x)}{g(x)}<M , \forall x$$ for some $M>0$?


title: Statistical analysis for relationships between three variables at two locations?

text: I have measured total organic carbon (TOC), total nitrogen (TN) and chlorophyll content in seawater at two depths i.e. upper and lower. Now I want to do statistical analysis of my data to compare the relation between the three - whether they are related or unrelated, whether they affect one another significantly or not. Please tell me what statistics test I should use.


title: Examples of Hetero- and Homoscedasticity with the Same Dependent Variable?

text: What are good examples to show that, even if a regression of $Y$ on $X$ is heteroscedastic, a regression of $Y$ on a different independent variable $Z$ could be homoscedastic? More formally, what are populations of triples $(Y_i, X_i, Z_i)$ for which all of the following hold: a) $E[Y_i | X_i]$ is an increasing (or decreasing) function of $X_i$; b) $Var[Y_i | X_i]$ is an increasing (or decreasing) function of $X_i$; c) $E[Y_i | Z_i]$ is an increasing (or decreasing) function of $Z_i$; d) $Var[Y_i | Z_i]$ is constant? I am interested in examples where it is known or plausible that these properties apply to a _population_ (ie they are not just quirks of particular samples that may not be representative). Although this scenario seems quite possible in theory, it is (to me at least) quite hard to find convincing examples. Eventually I thought of the following. $Y_i$ is the number of years individual i is in full-time education. $X_i$ is i's average family income during the period of i's education. $Z_i$ is i's month of birth (1,...,12) in terms of the school year. Suppose the individuals i are in a country in which: 1) children are required to start school on reaching a certain birthday (so that they start at different times of the academic year, depending when their birthday falls); 2) compulsory education ends at the end of the academic year in which children reach a certain age, and most education beyond the compulsory period also ends at the end of the academic year; 3) children from families with higher incomes are more likely to stay in education beyond the compulsory period. Then children from higher-income families are likely to spend longer in education on average, but the time they spend in education is likely to vary according to aptitudes and preferences, while children from lower-income families are more likely to spend in education only the compulsory period. Children born towards the beginning of the academic year are likely to spend longer in education on average, but there is unlikely to be much variation with month of birth in how much longer they spend (those born in different months being equally likely to be from higher-income families and to stay in education beyond the compulsory period). Are there better (simpler, more obvious) examples? Addendum: whuber has given an excellent theoretical answer showing how many examples can be generated. Answers involving specific examples with familiar variables (eg $Y$ = Income) would also be of interest.


title: machine learning model with 50% chance to be deterministic

text: I am trying to predict auction prices from a self gathered data set. I have titles available, and the price it was offered for. Now I noticed that whenever there is a 3 digit string of numbers in the title, there is about 48% that this is the actual price offered (I constrained to only gather data of the price between 100-800). How can I use this information in a further model, where I will make features of other text attributes (where 1 will mean the word exists in a training sample, and 0 will mean it does not exist in a training sample)? Is the best way to split the data up into where there exists a string in the title with numerical value between 100-800, versus where it does not exist, or is there a way to combine it in a single regression? i think what makes it different is that there is a 50% chance to be deterministic.


title: Neural network weights explode in linear unit

text: I am currently implementing a simple neural network and the backprop algorithm in **Python** with **numpy**. I have already tested my backprop method using central differences and the resulting gradient is equal. However, the network fails to approximate a simple sine curve. The network has one hidden layer (100 neurons) with $\tanh$ activation functions and a output layer with a linear activation function. Each unit also has a bias input. The training is done by simple gradient descent with a learning rate of 0.2. The problem arises from the gradient, which gets larger with every epoch, but I don't know why. Further, the problem is unchanged if I decrease the learning rate. EDIT: I have uploaded the code to pastebin: http://pastebin.com/R7tviZUJ


title: Normalize the ratio predictor variable by logit or not?

text: I have a Poisson GLM model               m1 <- glm(count ~ x, data = data, family = "quasipoisson")      The predictor variable x is a cumulative sum of ratios, i.e. real numbers between 0 and 1. (Not important for this question, but if you are interested why the sum is needed, it is because it is a transformation of a simpler but recursive model). ## Is it desirable to normalize the predictor variable by logit? i.e. instead of using cumulative sum of ratios, use cumulative sum of logit of the ratios (`x2`)?               m2 <- glm(count ~ x2, data = data, family = "quasipoisson")      I tried to look at the q-q plots (`plot(model)`) but they seem to be the same. However I tried to do another plot which is not present in the `plot(model)` output - I don't know if it's actually important or not. So see comparison of the two models on this plot: ![enter image description here](http://i.stack.imgur.com/oOWhx.png) Is this plot actually important for this decision? Does the second plot actually show heteroscedasticity? If yes, then is the original model (without transformation) possibly better according to these plots?


title: General approaches to model car traffic in a parking garage

text: a friend of mine has asked me to help him with predictive modelling of car traffic in a medium sized parking garage. The garage has its busy and easy days, its peak hours, dead hours opening hours (it is opened during 12 hours during weekdays and during 8 hours during weekends). The goal is to predict how many cars will enter the garage during a given day (say, tomorrow) and how are these cars supposed to be distributed over the day. Please point me to general references (preferably, publicly available) to strategies and techniques. Thank you


title: Change of variable with non bijective function

text: I have a logistic regression model: $softmax(WX)$ where $W$ is my parameter matrix and $X$ is my input. I want a density function over the outputs of that model. Say I know that my $X$ are distributed according to some density $p$. From the change of variables I know that the density of the outputs will be something like $p'(x) = p(x) det|J^{-1}|$ where $J = \frac{\partial{softmax(WX)}}{\partial{X}}$. However, this holds only for a bijective function $softmax(WX)$. What can I do if this is not the case - especially if $W$ is not square?


title: Non-repeated measures ANOVA of time series with different subjects per treatment in R?

text: I have the following experimental layout:   1. five different treatments - harvesting rates, ranging from 0 to 1, indicating proportion of branches per plant harvested   2. 75 plants, randomly assigned plants to each treatment, resulting in 15 plants per treatment   3. followed over 5 years, resulting in five data points per individual plant Now I want to analyse the data and see if the treatment (harvesting rate) has an impact on different measured variables (e.g., height of plant, cumulative number of stems harvested), and preferably also in which years they differ. Initially I thought using a repeated-measures ANOVA, but the subjects per treatment are different, but they are followed over time, which violates independence assumptions of normal ANOVAs. So which statistical test can I use here? I use R, so an example in R would be nice.


title: Choose best model between logit, probit and nls

text: I'm analyzing a certain dataset, and I need to understand how to choose the best model that fits my data. I'm using R. An example of data I have is the following:               corr <- c(0, 0, 10, 50, 70, 100, 100, 100, 90, 100, 100)      These numbers correspond to the percentage of correct answers, under 11 different conditions (`cnt`):               cnt <- c(0, 82, 163, 242, 318, 390, 458, 521, 578, 628, 673)      Firstly I tried to fit a probit model, and a logit model. Just now I found in the literature another equation to fit data similar to mine, so I tried to fit my data, using the `nls` function, according to that equation (but I don't agree with that, and the author does not explain why he used that equation). Here is the code for the three models I get:               resp.mat <- as.matrix(cbind(corr/10, (100-corr)/10))     ddprob.glm1 <- glm(resp.mat ~ cnt, family = binomial(link = "logit"))     ddprob.glm2 <- glm(resp.mat ~ cnt, family = binomial(link = "probit"))          ddprob.nls <- nls(corr ~ 100 / (1 + exp(k*(AMP-cnt))), start=list(k=0.01, AMP=5))      Now I plotted data and the three fitted curves:               pcnt <- seq(min(cnt), max(cnt), len = max(cnt)-min(cnt))      pred.glm1 <- predict(ddprob.glm1, data.frame(cnt = pcnt), type = "response", se.fit=T)      pred.glm2 <- predict(ddprob.glm2, data.frame(cnt = pcnt), type = "response", se.fit=T)      pred.nls <- predict(ddprob.nls, data.frame(cnt = pcnt), type = "response", se.fit=T)          plot(cnt, corr, xlim=c(0,673), ylim = c(0, 100), cex=1.5)     lines(pcnt, pred.nls, lwd = 2, lty=1, col="red", xlim=c(0,673))     lines(pcnt, pred.glm2$fit*100, lwd = 2, lty=1, col="black", xlim=c(0,673)) #$     lines(pcnt, pred.glm1$fit*100, lwd = 2, lty=1, col="green", xlim=c(0,673))      ![](http://i.stack.imgur.com/HbFm6.png) Now, I would like to know: what is the best model for my data?   * probit   * logit   * nls The logLik for the three models are:               > logLik(ddprob.nls)     'log Lik.' -33.15399 (df=3)     > logLik(ddprob.glm1)     'log Lik.' -9.193351 (df=2)     > logLik(ddprob.glm2)     'log Lik.' -10.32332 (df=2)      Is the logLik sufficient to choose the best model? (It would be the logit- model, right?) Or is there something else I need to calculate?


title: Predicting Data.

text: I have this data:               DAY  # of units Manuf     %      1          800           20     2          400           10     3         1200           30     4          800           20     5          400           10     6          200           5     7          200           5        * First column is day of week.   * Second column is # of units that we need to manufacture on that day.   * Third column is the % of units manufactured that day out of the total manufactured during the week. For example-total units manuf during these 7 days are 4000. On day 1 800 were manuf out of 400 i.e. 20%. My objective is to calculate the maximum total # of hours the machines can be rested during a specific day. I know that the total # of hours the machines will remain rested for these 7 days is 400. I just need to distribute these hours into 7 days. Now it's a safe assumption for us that the more units we manufacture the lesser we rest we give to the machines. Based on that this is how I solved this.               DAY      # of units Manuf %     Rest Hrs Var    Rest Hrs/Day            1          800           20       1.5x              27.3         2          400           10        3x               54.6         3         1200           30        x                18.2         4          800           20       1.5x              27.3         5          400           10        3x               54.6         6          200           5         6x               109.2         7          200           5         6x               109.2      So I assumed that the rest hrs for the day with max # of units be x. Based on that we can calculate the remaining values rest hrs var of each day. The sum of these hrs is 22x which is equal to 400. So 22x = 400 gives me x=18.2. Using this I calculated rest hrs for ech day. Please let me know if this the best way to do this. Alternative solution:               DAY      # of units Manuf  %      Inv of %    Min   Factor            SUM   %RestHrs  RestH       1          800           20         5         5  Sqrt(5-(5-1)=1.63   18.21    .089   27.2      2          400           10        10         5  Sqrt(10-(5-1)=2.77  18.21    .15    54.5      3         1200           30        3.33       5  Sqrt(3.33-(5-1)=1   18.21    .05    18.2      4          800           20         5         5  Sqrt(5-(5-1)=1.63   18.21    .089   27.2      5          400           10        10         5  Sqrt(10-(5-1)=2.77  18.21    .15    54.5      6          200           5         20         5  Sqrt(20-(5-1)=4.2   18.21    .23    109.1      7          200           5         20         5  Sqrt(20-(5-1)=4.2   18.21    .23    109.1      Please let me know whats the best solution for this kind of problem. Thank You


title: Smoothing constant in single exponential smoothing

text: I have some SKUs and I'd like to do a forecast using single exponential smoothing as a forecasting method, when should we go for small value of alpha (.05,.1,...) and when for bigger values(.8,.9,...)? Does it depend on the characteristics of the series?


title: AIC for non-nested models: normalizing constant

text: The AIC is defined as $AIC=-2 \log(L(\hat\theta))+2p$, where $\hat\theta$ is the maximum likelihood estimator and $p$ is the dimension of the parameter space. For the estimation of $\theta$, one usually neglects the constant factor of the density. This is, the factor that does not depend on the parameters, in order to simplify the likelihood. On the other hand, this factor is very important for the calculation of the AIC, given that when comparing non-nested models this factor is not common and then the order of the corresponding AICs might be different if it is not considered. **My question is** , do we need to compute $\log(L(\hat\theta))$ including all the terms of the density when comparing non-nested models?


title: Is Joel Spolsky's "Hunting of the Snark" post valid statistical content analysis?

text: If you've been reading the community bulletins lately, you've likely seen The Hunting of the Snark, a post on the official StackExchange blog by Joel Spolsky, the CEO of the StackExchange network. He discusses a statistical analysis conducted on a sample of SE comments to evaluate their "friendliness" from an outside user's perspective. The comments were randomly sampled from StackOverflow and the content analysts were members of Amazon's Mechanical Turk community, a market for work that connects companies to workers who do small, short tasks for affordable fees. Not so long ago, I was a graduate student in political science and one of the classes I took was _Statistical Content Analysis_. The class's final project, in fact its entire purpose, was to conduct a detailed analysis of the New York Times' war reporting, to test whether or not many assumptions Americans make about news coverage during wars were accurate (spoiler: evidence suggests they're not). The project was huge and quite fun, but by far its most painful section was the 'training and reliability testing phase', which occurred before we could conduct full analysis. It had two purposes (see page 9 of the linked paper for a detailed description, as well as references to intercoder reliability standards in the content analysis statistical literature):   1. Confirm all coders, i.e., readers of the content, were trained on the same qualitative definitions. In Joel's analysis, this meant everyone would know exactly how the project defined "friendly" and "unfriendly."   2. Confirm all coders interpreted these rules reliably, i.e. we sampled our sample, analyzed the subset, and then statistically demonstrated our pairwise correlations on qualitative evaluations were quite similar. Reliability testing hurt because we had to do it three or four times. Until -1- was locked down and -2- showed high enough pairwise correlations, our results for the full analysis were suspect. They couldn't be demonstrated valid or invalid. Most importantly, _we had to do pilot tests of reliability before the final sample set._ My question is this: Joel's statistical analysis lacked a pilot reliability test and didn't establish any operational definitions of "friendliness". Was the final data reliable enough to say anything about the statistical validity of his results? For one perspective, consider this primer on the value of intercoder reliability and consistent operational definitions. From deeper in the same source, you can read about pilot reliability tests (item 5 in the list). Per Andy W.'s suggestion in his answer, I'm attempting to calculate a variety of reliability statistics on the dataset, which is available here, using this command series in R (updated as I calculate new statistics). Descriptive statistics are here Percentage agreement (with tolerance = 0): 0.0143 Percentage agreement (with tolerance = 1): 11.8 Krippendorff's alpha: 0.1529467 I also attempted an item-response model for this data in another question.


title: How to understand / visualize the error surface in online learning algorithms

text: I have a question about the shape of the error surface for online gradient descent algorithms. Take into account that I am trying to translate my specific question into a more general and idealized situation so that you can understand it. This means that the situation and constraints you will read here may look a bit odd and simplistic. Imagine the following situation:   * I have only one parameter to optimize. This is really like having an extremely simple connectionist network with 2 units (input and output).    * The data that are provided to the system are also extremely simple: Both the input and output units can be thought as binary elements (they are either active, value=1, or inactive, value=0). On each trial, I present one exemplar consisting of a pair of numbers, corresponding to the activation of either unit.    * Moreover, let's make it simpler. Imagine that the input unit is always active (value=1) and I just need to find the best parameter (the weight) to predict the current outcome (target = 1 or 0). So each exemplar actually consists of a 1 and a target to be predicted (either 1 or 0).    * My learning algorithm is a variant of the delta rule that performs steepest descent on the error function.    * Note that the parameter is updated online, after presenting each exemplar (as opposite to the more usual batch update method). This is the main complication of my idealized problem. I need my learning algorithm to adapt to unexpected changes in the input stream.  Now, the question: Given this situation, I have a problem "visualizing" and understanding the real task the algorithm is carrying out. In particular, I am wondering what exactly the error function looks like, even in this extremely simplified situation. My intuition tells me that the curvature of the surface will depend on the consistency in the input stream. Say, if many consecutive trials consist of the same target (e.g., "1-1, 1-1, 1-1..."), then the surface would be flatter, while if they point to inconsistent directions (e.g., "1-0, 1-1, 1-0...") it would look more curvy / steep.   * **Is this intuition right?**   * **Would you suggest any metaphor or image to help me "visualize" the problem the online algorithm is actually solving?**


title: Maximam r distance for Ripley's K-function

text: I am using R's package spatstat to study the locational pattern of conflict events in Africa (around 8.000 points) using point pattern analysis techniques. I was able to obtain the plot of g(r), the pair correlation function, using the contour of the African continent as my study region. The size of the enclosing rectangle is around 8.546km by 10.423km. Spatstat does a nice job and gives me the plot of g(r), but I noticed that the computation is not performed for r > 625km. Looking at this reply and in the function's source code, I understood that spatstat uses a rule of thumb stating that the biggest r for which to compute Ripley's K-function and its variations should not excede 25% of the shortest side of the rectangle enclosing the study region. **My questions are thus the following:**   1. Is this rule of thumb also valid for complex polygons as it is for rectangles?   2. What happens to my K-function if I compute it for a larger r than what is advised?   3. Is there any literature addressing this specific problem (size of r, polygonal windows)? Thank you in advance for your help!


title: libsvm "reaching max number of iterations" warning and cross-validation

text: I'm using libsvm in C-SVC mode with a polynomial kernel of degree 2 and I'm required to train multiple SVMs. Each training set has 10 features and 5000 vectors. During training, I am getting this warning for most of the SVMs that I train:               WARNING: reaching max number of iterations     optimization finished, #iter = 10000000      Could someone please explain what does this warning implies and, perhaps, how to avoid it? I also want to apply cross-validation for my models in order to determine the best choices for gamma and C (regularization). My plan is to just try every combinations of these 10 values: 0.00001, 0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000 for both parameters and see which combination produces the best accuracy during cross-validation. Is this enough? Should I use more values in this interval, or should I choose a wider interval?


title: predictors and varImp functions in caret packages, should results match for a glmnet model?

text: I am using the train function in caret package. I have 55 predictors and a continuous outcome. The model that I selected is glmnet. Looking the results of using the "predictors" function-just to check what variables are involved in the final model- it seems that 29 out off 55 variables are taking part in the model. Right. The problem comes up when after that, I used varImp to also list the variable importance (useModel="lm", but in fact if I leave the default and ROC analysis is performed, I had exactly the same values for importance). In the importance values list, there is just n=1 variables with weight equal to 0. And one of the variables that previously the "predictors" function chose as one of the predictors used in the model (smoke), now appeared as one of the less weigthed in varImp list. So essentially, I have different results from "predictors" and "varImp". Does this usually happen or am I wrong? I guess the most correct is the results I had from "predictors' function, as varImp is an external process.. Thanks guys for your help! Raquel


title: How to find a conditional probability using copula-based Markov process?

text: I have a monthly time series of a water quality parameter. I used copula-based Markov process of C(Y(t), Y(t-1) and I forecasted the mean behavior of Yt by following equation: ![enter image description here](http://i.stack.imgur.com/pGjBe.jpg) Now, I need to find the probability that for a given month, my water quality (Yt) will be lower than a predefined value (a critical value that has been set up by regulations). Can you please help me how to find this probability? If there are methods other than copula-based Markov process that I can use please share it with me. I used this paper as my reference: Chen, X., Fan, Y., 2006, Estimation of copula-based semiparametric time series models


title: How to check if my data fits log normal distribution?

text: I'd like to check in R if my data fits log-normal or Pareto distributions. How could I do that? Perhaps ks.test could help me do that, but how could I get the alpha and k parameters for Pareto distribution for my data? Any help will be appreciated.


title: Bootstrap for estimating parameters with only one sample

text: Let say we have only one sample in which we are interested only on one statistics c of this sample. This is usually true in practice when we usually have only one data set. So, most probably the usual estimating method like MLE cant be applied. And let say we want to estimate 2 parameters (i.e. a & b) of a certain dist which depend on c. Is bootstrap technique will work for this? Is there any theoretical proof for such practice? Any help will be appreciated. tsukiko


title: Chi-square analysis

text: I am trying to determine association between my independent variables (categorical, nominal) and dependent variables (categorical, nominal). What other statistical analysis I can use aside from chi-square? Almost all my variables have two levels except two independent variables that had more than 2 levels but chi-square analysis didn't show any significance association for those variables. Also, can I determine the magnitude of the significance by using chi-square? I have a few variables that chi-square showed significant associations, $p<.05$, but I don't know if I need to run further tests to determine the magnitude of association or chi-square should be able to tell me that without further analysis.


title: What analysis should I perform?

text: What analysis should I perform if I am trying to find if there is an association between socio-demographic variables and my dependent variables. All my variables are categorical nominal. Also I am trying to see if there is a change in my dependent variables in my control group at baseline and follow up. Again, all variables used here are categorical nominal. Please let me know what analysis is best at answering my research questions fully.


title: repeated measures ANOVA with gender as between subject factor

text: I have conducted an intervention study on preschool children on behavior problems, social-emotional competence with control and experimental groups. I have analysed the results using ANCOVA to find out main intervention effects on measures. Now I want to see if there are any gender differences in intervention group on the measures. on one variable I have got significant gender*time effects. My question is that how would I interpret which group has performed better on intervention, males or females?


title: Is it acceptable to run two linear models on the same data set?

text: For a linear regression with multiple groups (natural groups defined a priori) is it acceptable to run two different models on the same data set to answer the following two questions?   1. Does each group have a non-zero slope and non-zero intercept and what are the parameters for each within group regression?   2. Is there, regardless of group membership, a non-zero trend and non-zero intercept and what are the parameters for this across groups regression? In R, the first model would be `lm(y ~ group + x:group - 1)`, so that the estimated coefficients could be directly interpreted as the intercept and slope for each group.The second model would be `lm(y ~ x + 1)`. The alternative would be `lm(y ~ x + group + x:group + 1)`, which results in a complicated summary table of coefficients, with within group slopes and intercepts having to be calculated from the differences in slopes and intercepts from some reference. Also you have to reorder the groups and run the model a second time anyway in order to get a p-value for the last group difference (sometimes). Does this using two separate models negatively affect inference in any way or this this standard practice? To put this into context, consider x to be a drug dosage and the groups to be different races. It may be interesting to know the dose-response relationship for a particular race for a doctor, or which races the drug works for at all, but it may also be interesting sometimes to know the dose-response relationship for the entire (human) population regardless of race for a public health official. This is just an example of how one might be interested in both within group and across group regressions separately. Whether a dose- response relationship should be linear isn't important.


title: Effect of Age-Period-Cohort on prediction and interactions

text: I am performing a longitudinal analysis and I am curious if the predictors we are including in the model will introduce any unexpected effects. We have subjects with multiple points of follow-up and they are entered into the study at various times and ages. The main hypothesis we want to test is if there is an effect of what age you entered the study and how long you've been in the study (e.g. at a given age, if you've been in the study for 10 years or 20, is there a difference in your measures). It is also possible that different criteria were used to allow subjects to enter the study depending on when they entered. Thus, we are adjusting for their age at entry, age at measurement, and decade of entry. To test the main hypothesis, we are looking at the age-at-entry * age-at-measurement interaction. I think people usually do age at entry and time since entry, but this should give identical predictions, and time itself is not directly interesting to us, but both age variables are. My understanding of APC models is that prediction is unaffected. But the problem is in the interpretation of individual effects -- each coefficient for age at entry and age at measurement is confounded with the duration of study. Similarly, I think that age at entry and decade of entry will be confounded with birth year, though we assume birth year has no effect. Are any new problems introduced with when including interactions? Is including all 3 of these variables problematic? I feel like we might be doing something profoundly stupid with this model, yet I'm not quite sure how else to test the hypotheses and control for selection criteria at entry.


title: Comparing multiclass classification algorithms for a particular application

text: I am assessing a bunch of classification algorithms for a specific application with multiple classes. The classification algorithms that I am considering are:   1. Multinomial Logistic Regression (Matlab's 'mnrfit')   2. Multiclass SVM (K. Crammer and Y. Singer. On the Algorithmic Implementation of Multi-class SVMs, JMLR, 2001.). I will use the code provided by the authors since Matlab's 'svmtrain' only does binary classification.   3. Neural Networks (Matlab's 'nprtool')   4. Decision Trees (C4.5 and CART from Matlab's 'classregtree')   5. k-Nearest Neighbors (Matlab's 'ClassificationKNN')   6. Naive Bayes Classifier (Matlab's 'naiveBayes')   7. Discriminant Analysis (Matlab's 'ClassificationDiscriminant')   8. Random Forests (Matlab's 'TreeBagger') I have the following questions:   1. Have I omitted any "obvious" multiclass classification algorithm that's a must-try? Or, are there any binary classifiers that can easily be used for multiclass with one-vs-all method.    2. Which of these are linear classifiers and which are non-linear classifiers? I know 3, 4 and 5 are non-linear by nature and 2 can be non-linear with the kernel trick.   3. Which of these are discrete classifiers and which are probabilistic? For example, logistic regression gives a probability for each class, while decision trees give exactly one class. On a related note, if I want to use ROC curves for comparison, how do I compare a discrete classifier with a probabilistic one, as the former gives only a point on the ROC plot?   4. Which of these are deterministic classifiers and which are stochastic? In other words, which classifiers will yield exactly the same results for multiple runs? This is important for me: For example in logistic regression I can be done with just one run, whereas for neural networks, I will have to train the net multiple times to avoid biased results. Finally, are there any good papers which show how to compare various classification algorithms. I found this one which is quite good.


title: Choosing sample size to achieve pre-specified margin-of-error

text: If I want to achieve a margin-of-error of <= 5 % for a representative population sample, how large a sample do I need when: The interviewees are picked from _X_ regions and from _Y_ age groups? That is, how many samples from each age group & region? I know how many people live in each region and how they are distributed in the _Y_ age groups in each region.


title: generalized linear model help in SPSS (and in general)

text: My dependent variable is not normally distributed. I chose the Gamma log link, and I hope this is correct. I wasn't entirely sure with this because the dependent variable is actually a whole lot of percentages. I know this isn't technically continuous, but if I choose Poisson, it excludes most of my data. so question 1. if this is the wrong distribution, what is going wrong with the others? Next, when I run the model, I have a few factors which are sex, stage (5 categories) and number of offspring (1 or 2 only). My covariates are weight and the other dependent variable, which are blood cell counts. I then selected all of these for the model, with an interaction between sex and stage, sex stage and offspring, sex and offspring. I have also played around by taking these in and out, and not having interactions at all. However, I have done single analysis on this before, and am interested in sex and stage related to my dependent variable. I also checked the box for pairwise comparisons, and Bonferroni. So these are the questions I have: 1) if my interactions are either approaching significance or not significant, can I still look at the pairwise Bonferroni figures and quote the significant values from that? 2) If I take all my interactions out, but want to look at how my dependent variable is affected by stage in females only, can I exclude male cases, and run two separate models (then exclude females). If I explain this in a paper, saying that I ran three different (looking at the original also) models and use the results for all, is this valid? Or would this open up the possibility of errors. 3) I have been told by one person that if I have an interaction (significant or not), I can no longer look at the single effects. Is this true? 4) If my interaction between sex/stage and my dependent variable is not significant, should I exclude it from the model? The AIC doesn't change very much between them. If I am allowed to do 1 from above, there are a couple of interesting things happening. 5) why are my tests of model effects always significant but my parameter estimates aren't? is it because of the various levels in some of my factors?


title: Does SPSS rescale dendrograms?

text: A colleague and I have been clustering some data in SPSS (v19) and R (2.15), respectively. Using the same distance metric and agglomeration method, we get identical merge orders/agglomeration schedules in both programs, and the dendrograms have very similar shapes, but the actual height values are quite different. On this particular data set, the R dendrogram is about 150 units tall, but the SPSS dendrogram is only 25 units tall. This is somewhat $\ldots$ unsettling, obviously. While thumbing through some manuscripts and webpages, I noticed that all but one of the SPSS-derived dendrograms were also exactly 25 units tall (example 1, 2, 3; we also have a pretty hefty pile of papers that are unfortunately all pay-walled, including the single counter-example). The caption on the SPSS output says something about rescaling, but the documentation is oddly silent about if, how, and why SPSS might be rescaling the dendrograms. Could someone please confirm that SPSS does rescales dendrograms (and rescales them onto [0,25])? For extra credit, is there a way to turn this rescaling off? It seems to cause SPSS to cut our dendrogram a few levels above the leaves.


title: Testing hierarchical data

text: It is easy to compare _paired_ data, using paired t-test. But suppose this pairing is hierarchical. Below is an example: The stem cells of two mouse strains S1 and S2 are cultured in two different culture conditions C1 and C2, having 3 replicates for each Strain/Condition (2 x 2 x 3 = 12 samples total). We have examined the expression level of N genes in all of these samples. We want to identify the genes which are significantly (p < 0.02) higher expressed in culture condition C1 in any strain. The case was easy if we had no replicates: just pairing the identical samples in two conditions C1 and C2 and running a paired t-test. If we had replicates but no strains the problem was easier: just running a normal t-test on two population of cells cultured in C1 and C2. But how to deal the case while our data is hierarchical?


title: Clustering time series

text: I want to create forecasting for a large quantity of time series. Since they are too many, I am thinking on reducing my data by clustering it into to similar groups. However, I am using SPSS modeler and it is not possible to cluster time series (only static data). Do you think it makes sense to apply clustering on static data and fit the forecasting model on its centroide? Which software do you suggest to cluster the time series directly? Thanks!


title: M-PLUS and HLM 7 for Hierarchical Linear Modeling

text: Many of my colleagues are considering using the M-PLUS package to carry out Hierarchical Linear Modeling, rather than the HLM 7 package. Both packages are designed by very well respected figures in the field of Hierachical Linear Modeling. Have any users seen distinct statistical advantages in using one package or the other? Thank you in advance, Stephen Brand


title: Understanding the random effect in linear mixed models (lme4, R)

text: I am familiar with linear regression models but the random section of linear mixed models just melts my mind. I did find an excellent guide that could have helped me but the languageR package is not compatible with newer versions of lme4 so I've been unable to implement it in my work. For me the fixed effects are very understandable (below lactation and a higher yr2 value both contribute to a higher weight but the lactation effect is more consistent which results in a higher t-value). The first problem is to understand what I am actually putting in. To a certain extent I understand that `(1|P$grupp)` means that the mixed model add to the base line (intercept) while `(P$grupp|P$lweek)` mean that belong to a group is expected to affect the average weight increase (or decrease) while `P$lweek` adds to the baseline value. But why does all tutorials seem to favor write ups like `(1+P$fgrupp|P$lweek)` rather than `(P$grupp|P$lweek)`? Now on to the actual output (see below for full output). I've used the following models (sorry for the Swenglish but the sample is the weight of cows `P$vikt` is the weight at certain time points and `P$lweek` is the time since a calf was born, `P$fgrupp` is a factor telling if the cow belongs to feed group 1,2 or 3):               Formula: P$vikt ~ P$lweek + P$laktation + P$yr + (1 | P$fgrupp)            #$     Formula: P$vikt ~ P$lweek + P$laktation + P$yr + (1 + P$fgrupp | P$lweek)      Where I understand it as the first one being rather useless (essentially it tells us that the average weight of the cows isn't affected of which feed group it belongs too). This is reflected by fgrupp having variance 0 in the first formula below. The second is more interesting as the `P$fgrupp|P$lweek` as I understand it should show if different feed groups affect the weight increase of cows as function of the time. But I am really not competent enough to understand the input. I understand that variance somehow mean that belonging to group2 or 3 explain some of the variation in the growth curves but I really don't understand how to interpret this.               Random effects:      Groups   Name        Variance Std.Dev. Corr              P$lweek  (Intercept)   13.068  3.6149                                     #$               P$fgrupp2     77.230  8.7881  1.000                              #$               P$fgrupp3     81.188  9.0104  1.000 1.000                        #$      Residual             4031.831 63.4967                   Number of obs: 1048, groups: P$lweek, 84      **Full output**               #First model#     Linear mixed model fit by REML      Formula: P$vikt ~ P$lweek + P$laktation + P$yr + (1 | P$fgrupp)            #$        AIC   BIC logLik deviance REMLdev      11703 11732  -5845    11698   11691     Random effects:      Groups   Name        Variance Std.Dev.      P$fgrupp (Intercept)    0.0    0.000                                      #$      Residual             4139.9   64.342       Number of obs: 1048, groups: P$fgrupp, 3                                   #$          Fixed effects:                 Estimate Std. Error t value     (Intercept)  509.593      4.683  108.82     P$lweek        1.028      0.105    9.79                                    #$     P$laktation   22.789      1.454   15.67                                    #$     P$yr2         35.294      4.093    8.62                                    #$          Correlation of Fixed Effects:                 (Intr) P$lwek P$lktt                                           #$     P$lweek     -0.560                                                         #$     P$laktation -0.636  0.030                                                  #$     P$yr2       -0.240 -0.034 -0.141                                           #$               #Second model#          Linear mixed model fit by REML      Formula: P$vikt ~ P$lweek + P$laktation + P$yr + (1 + P$fgrupp | P$lweek)  #$            AIC   BIC logLik deviance REMLdev          11707 11761  -5842    11693   11685         Random effects:          Groups   Name        Variance Std.Dev. Corr                  P$lweek  (Intercept)   13.068  3.6149                                     #$                   P$fgrupp2     77.230  8.7881  1.000                              #$                   P$fgrupp3     81.188  9.0104  1.000 1.000                        #$          Residual             4031.831 63.4967                       Number of obs: 1048, groups: P$lweek, 84                                   #$          Fixed effects:                 Estimate Std. Error t value     (Intercept) 508.2291     5.1770   98.17     P$lweek       1.0662     0.1192    8.94                                    #$     P$laktation  22.6525     1.4459   15.67                                    #$     P$yr2        35.6343     4.0848    8.72                                    #$          Correlation of Fixed Effects:                 (Intr) P$lwek P$lktt                                           #$$     P$lweek     -0.627                                                         #$     P$laktation -0.570  0.025                                                  #$     P$yr2       -0.224 -0.018 -0.136                                           #$


title: What's the tractable data size for Sparse PCA or LASSO?

text: I'm up to perform certain kinds of sparse decomposition methods on my dataset. However, I'm not sure: what's the tractable data size for the Sparse Decomposition methods? The dataset is a $10^3\times10^5$ binary matrix, and the expected methods are Sparse PCA, which aims to find a decomposition with minimized coordinate L1-norm and positive constraint. I found these things useful: penalized package for R language. nsprcomp pacakge for R


title: How to measure student skill level based on a set of questions answered at different points in time and of different difficulties levels?

text: I am working on an e-learning system with a friend for our final year (Computer Science) project which is part of the under-graduate programs mandatory 'courses'. I have a question about making inferences and therefore gauging the skill level of the 'subject (student) at hand' We have gathered the following data thus far for every attempt made at answering a question:   1. Unique ID of the question   2. Attempt number   3. Outcome of the attempt (Correct/Incorrect)   4. Difficulty of the question (assigned a nominal score between 1-5)   5. Date and Time at which the attempt was made **How can I use the entire data set regardless of what point in time the answers were given and judge the skill level of the student from that particular model?** I know I could just look at the last X number of attempts for a section of a particular difficulty level and pick up a trend from that, but is there a better way to do this?


title: Statistical test for zero valued variable

text: I am testing a product (liquid bag) for leakage and my primary variable is the amount of leakage in ml. I have a specification that states, on average, leakage must be below XXml. In total I tested 20 bags. Now my problem is this: none of the bags leaked so essentially my data consists of 20 values all equal to zero. Since I have no variability in the data, a t-test offers no utility (average CI for leakage 0 +/- 0 ml????). Aside from treating the data as attribute (leakage/no leakage) and getting a confidence interval for the proportion leaked, does anyone have any opinion on what type of test I perform if I want to keep the data variable (and not as attribute)?


title: ART neural network disambiguation

text: I have an assignment to implement the adaptive resonance theory (ART) type network (as part of a bigger project). I have red a lot of Internet resources on the topic and I think I've got the essence of it, but I am not sure. So far I have the following findings and questions:   * The ART1 network has two layers (input and output) which are fully connected in both directions.   * The input passes through the bottom-up connections, and the output neuron with highest value is the winner. Than the input is somehow compared against the prototype stored in the top-down connection weights (not sure how this is done, does the out-value of the winning output neuron pass through the top-down connections or not, it looks to me that the top-down connections are not really proper connections by definition)   * Does the ART1 network have a learning phase, and operating phase (like a multi-layer perceptron) or is it adjusting its weights constantly and it is up to the user to decide when is the "learning" phase over?


title: Pre- and post- presentation data analysis

text: I have carried out some research on whether awareness and knowledge of HIV had been increased on university students by giving a questionnaire to students before and after a presentation on HIV. Only 140 students completed the pre questionnaire and only 129 completed the post questionnaire after the presentation. So how is it possible to compare whether knowledge in students increased following the presentation as the numbers are unequal before and after? Could I compare the mean or mode score before and after even though the number of participants was different?


title: Estimating variance from sequence of random variable

text: Given $X\sim N(0,\Omega)$. Suppose that we can construct a sequence $\\{X_n\\}$ based on the observation such that $\\{X_n\\}\to X$ in distribution. My problem is to estimate $\Omega$ consistently using this sequence. If $\operatorname{var} X_n$ converges to a "finite" matrix, then of course $\operatorname{var}X_n$ is the estimator that we want using the dominated convergence theorem. But, as we know $\operatorname{var} X_n$ can go to infinity, so in this case, what should I do? Any help will be appreciated. Thanks


title: Sufficient statistics of posterior (with Poisson data)

text: Suppose that, for year $t$, the data $y$ is Poisson with mean $a + bt$. Assume also a uniform prior on $(a,b)$. If we have $n$ years of data then I think the posterior for $(a,b)$ will be \begin{align*} p(a,b | y) &\propto p(y|a,b)p(a,b) \\\ &\propto p(y | a, a) \\\ &=\prod_{i=1}^n {\rm Poisson}(y_i | a+ bt_i) \\\ &\propto \prod_{i=1}^n (a+ bt_i)^{y_i} e^{-(a+ bt_i)} \\\ &= e^{-(na + b\sum_{i=1}^n t_i)}\prod_{i=1}^n (a + bt_i)^{y_i} \\\ \end{align*} What I'm struggling with is how to find the sufficient statistics. At first I thought maybe continue like this \begin{align*} e^{-(na + b\sum_{i=1}^n t_i)}\prod_{i=1}^n (a + bt_i)^{y_i} &=e^{-(na + b\sum_{i=1}^n t_i)}e^{\sum_{i=1}^n y_i \log(a+bt_i)} \end{align*} which seems promising but I'm confused about how to deal with the argument of the log since it also depends on $i$. Any suggestions?


title: Covariance between two sample means of correlated data

text: I have two sets of random data $X=\\{x_1,...,x_N\\}$ and $Y\\{y_1,...,y_N\\}$ both of length $N$. The sets are autocorrelated such that the correlation between $x_i$ and $x_j$ depends only on $|i-j|$. From both of these I can find the sample mean, $$ \bar{X} = \frac{1}{N}\sum_{i=1}^N x_i $$ and similarly for $Y$. I believe I can find the variance of each mean as follows, $$ \text{var}(\bar{X}) = \frac{S^2}{N}\frac{N-1}{\frac{N}{\gamma_2}-1} $$ where $S^2 = \frac{1}{N}\sum_{i=1}^N(x_i-\bar{X})^2$ and $\gamma_2 = 1+2\sum_{j=1}^{N-1}(1-j/N)\rho_j$, with $\rho_j$ being the autocorrelation function at a lag $j$. My question is this. What is $\text{cov}(\bar{X},\bar{Y})$ if the data sets are also correlated? The correlation between the sets should be the same for all pairs i.e. $\text{cov}(x_i,y_i) = \text{cov}(x_j,y_j)$. Will it be some generalisation of $\text{var}(\bar{X})$ or can I calculate it from the variances of the means and the correlation alone?


title: mixed model with both random effects and specified covariance structure

text: A linear mixed effects model allows you to specify random subject-level effects, or a (non-zero) covariance structure for within-subject residuals, or both. I've read a number of separate explanations of subject-level random effects and residual covariance structures, but I don't understand how these relate to each other or whether it's ever appropriate to use both in the same model. For example, if you specify a random intercept by subject and also a compound symmetry covariance structure, those both predict that the observations for a given participant will tend to be higher or lower than the relevant fixed- effects mean. Are they completely redundant with each other? Is it just an artifact of the software interface that I'm even able to include both? On the other hand, you could specify a subject-specific intercept and an autoregressive covariance pattern, which don't entirely overlap. Does this mean that an autoregressive pattern will be used to model each data point's deviations from the value predicted by the fixed and random effects?


title: False alarm probability: general vs. element by element

text: Suppose I make measurements of variable $P$ and a sequence of variables $D^1, D^2, \ldots$. There is correlation with $P$ and $D^i$, but not with $D^i$ and $D^j$. I can make a detection metric that looks like: $$\boldsymbol{z} = [P, D^1, D^2, \ldots, D^N]^T$$ or of the form $$\boldsymbol{z}^1 = [P, D^1]^T; \quad \boldsymbol{z}^2 = [P, D^2]^T \ldots ; \quad \boldsymbol{z}^N = [P, D^N]^T$$ In a hypothesis testing framework, the probability of false alarm is given by $$ P_F = \int p_{\boldsymbol{z}|H_0}(\boldsymbol{\xi}|H_0) d\boldsymbol{\xi}$$ **Question** : Is $P_F$ for the complete $z$ equal to the element wise false alarm probability for the $z^i$? That is does the following hold: $$P_F \stackrel{?}{=} 1 - (1-P_{F1})(1-P_{F2})\cdots(1-P_{FN})$$ **Thoughts** : I was trying to use the [VanTrees] general Gaussian problem framework where I model everything as Gaussian, but I found the analysis intractable. With the special form of the covariance matrix (zeros everwhere but along diagonal, first row, and first colum), I'm wondering if there is some property to exploit. edit: typo