title: Differences between groups using proportions

text: I am doing two analyses. First, I want to compare performance of two groups (p<.05), with the independent variable being proportional data (total number of categorical responses/total number of responses). I need to covary age. I am told that normal MANCOVA cannot be done because of homogeneity assumption violations. Arc-sin transformations were recommended, but one site suggested I multiple the transformation times 2. Why? Next, I want to identify the best predictor(s) of the proportional categories to a continuous variable (or binary if I transform to Hi/Lo, but I'd rather not), again, all independent variables all being proportions. Thank you in advance for your assistance.


title: Mixed-effects for repeated measures in Matlab -- how to choose a model?

text: This post follows on from this one. I am trying to analyse a dataset: $S$ subjects were tested for $A$ measures in $B\times C$ conditions. I have worked out that I can perform a repeated measures analysis by using a mixed-effects design where the subjects are random factors. Then I am left with a '3 way' $A\times B \times C$ ANOVA, which I perform with MATLAB's `anovan` function. Where I am stuck is: **what model structure should I be using, and what do I report? Specifically, what, if not all, (random) interactions, i.e. $?\times S$, should I be including in the model?** (If I don't include all $?\times S$ interactions, MATLAB reverts to using the MSE of the error degrees of freedom as the denominator in the F-statistic calculations. I don't know whether that's good or bad.) **Specific questions:**   1. Nomenclature: Am I performing an $A \times B \times C$ repeated measures ANOVA, or an $S \times A \times B \times C$ mixed-effects ANOVA, or maybe an$A \times B \times C$ repeated-measures mixed effects ANOVA?   2. Should I be including all the $?\times S$ interaction terms in the analysis, or is this a choice-of-model question?   3. Should I report the effects of $S$ and/or the $?\times S$ interactions? Do I want some kind of $\eta^2$ effect measure to communicate how much variance is due the subjects?


title: How to compare two implementations of web service by parameters of distributions of their response times?

text: We are migrating web services to new platform. Current WS will be wrapped behind new platform's layer. We would like to define criteria for response times of wrapped WS that should be met so we can say that wrapped WS is at least as good as original. My idea is that these criteria should be based on parameters of distribution of response times of unwrapped original. These can be obtained from access logs. At first I thought that define median as only criterion would be enough. Then I added SD . But now I'm confused because I'm not sure how to compare these two distributions. If I got distribution 1 with parameters m1 (median) and SD1 and distribution 2 with parameters m2 and SD2 how can I say what distribution is "better". My idea is that its somehow similar to two sets of measurements of length of two objects. First set contains measurements of object 1 and has parameters m1 and sd1 and other contains measurements of object 2 and has parameters m2 and sd2. I think that one can't say that object 1 is longer than object 2 just because m1 > m2 or can one? I'm not even sure if median and SD are right parameters to describe distribution of response times of WS. According to this answer http://stats.stackexchange.com/a/46374/44256 response times follow Log-Normal distribution. What parameters of distribution should I use for comparison of Log-Normal distributions? What about outliers, should I filter them out before comparison? Are there some general rules that can be followed in comparison of two distributions that are not dependent on type of distribution?


title: Question on the usage of regularization in applied statistics/science

text: I was reading the paper _``A significance test for the lasso''_ by Lockhart, Tibshirani _et al_ and was considering the issue of applying regularization in the applied sciences (for example, behavioral sciences). When do statisticians consider it acceptable to apply regularization methods like LASSO to shrink coefficients? There seems to be a clash between statistical convenience and external scientific theory here. Or does one use regularization only when it fits a priori knowledge of what we want to model? I realize this answer might vary between statisticians and I hope the question is not too soft/vague. Any insights are appreciated.


title: Is a Gaussian-Gaussian RBM just a linear model?

text: The 'conventional' configuration of RBMs are Binary-Binary and Gaussian-Binary (and sometimes Binary-Gaussian) units. Although it is possible for both the visible and hidden units to be gaussian, wouldn't a Gaussian-Gaussian RBM just resemble a linear model, since there is no non-linearity in the networks units anymore? Thus, stacking them would not not have the same benefit as for, say, Binary-Binary RBMs. And when using them for dimensionality reduction, a simple PCA would achieve better results? Am I missing any significant points in the training of RBMs or are Gaussian- Gaussian RBMs just that limited?


title: Queries on the Bayesian method

text: Currently I am working on a bayesian model framework and have questions related to the philosophy of using such techniques of modeling.   1. How do I know that the prior which I have captured from the experts is valid. There are parameters of the model which has captured a very wide range - say, 10% to 90%. It does not give me comfort, rather, it may show that the expert panel inputs missed out on the clear range. Is there any method out there which may allow me to check this?    2. We do not have enough data to work on, thus, the Bayesian framework. When can we say that Bayesian analysis is not required and the whole analysis/ model can be done using data. Is there any threshold on data availability/ philosophy where it indicates the transition from Bayesian to classical? (I understand that they are techniques from two different school of thought, so Bayesian is used in our case for lack of data)   3. The conjugate prior method for a beta-beta-beta model gives alpha estimate for posterior = a1+a2-1 beta estimate for posterior = b1+b2-1. the question which haunts me is, how will the effect of a larger data set be captured in the posterior parameters, if this model is being used post the required amount of data is available? It would be great if someone can answer my questions. If further clarification is required on my thoughts/ questions... please do let me know. Thanks!


title: What causes standardised regression coefficients in a confirmatory factor analysis to be greater than one?

text: I performed a CFA to reduce my items into factors. However, due to the violations of non-normality, the results were really weird, and transformations didn't help. There were some standard regression coefficients that were negative and ones that were over 1 despite the factors itself having good cronbach alphas.   1. Where could I have gone wrong in my CFA? What's the meaning behind getting path coefficients over 1?   2. Would it be suitable to even include such a bad-fit model into the write-up of my results?


title: Why are all the spearman's rho correlations with one variable low when the other variables all have moderate correlations?

text: I was wondering what could be the reasons for a very low Spearman correlation, as low as 0.01-0.06. Basically no correlation even exists! All my other variables have moderate correlations except any correlations with this one particular variable.


title: Regression by the model $y_t = β + u_t$

text: I have 30 data points and should calculate the residuals $u_t$ from them. My regression model is: $$y_t = β + u_t$$               Data Points     -15,6     -21,6     -19,5     -19,1     -20,9     -20,7     -19,3     -18,3     -15,1     -14,1     -14,9     -26,4     -26,87     -23,11     -25,38     -20,60     -14,96     -10,54     -4,72     -0,36     -4,17     4,14     1,09     6,81     17,27     19,24     24,35     25,35     30,02      I calculated by the ordinary least square method that β is $β=-y_t$ so from all the data points $β=26,8722467$ and then $u_t = y_t-β$ Is this correct?


title: Interpretation of incidence-rate ratios

text: So, I want to fit a random effects negativ-binomial model. For such a model STATA can produce exponentiated coefficients. According to the help file such coefficients can be interpreted as incidence-rate ratios. Unfortunately I am not a native English speaker and I do not really understand what incidence- rate ratios are or how I could translate them. So my question is, how can I interpret incidence-rate ratios. E.g.: If the model gives me an incidence-rate ratios of .7 for one var. would that mean that the number of expected observations (counts) on the dependent var. changes by .7 if the independent var changes by one unit? Can anyone help?


title: resettest in lmtest: unexpected behavior on perfect fit

text: The function `resettest` in the `R` package `lmtest` runs the Ramsey RESET test. A student today brought to my attention that this function behaves oddly when it is handed a regression which displays perfect fit:               library(lmtest)          x <- seq(1,100)     y <- x          reg1 <- lm(y~x)     summary(reg1)          # Should match but don't     summary(lm(y~x+I(x^2)))     resettest(reg1,power=2,type="regressor")          set.seed(12344321)     x <- rnorm(1000)     y <- x          reg2 <- lm(y~x)     summary(reg2)          # Should match but don't     summary(lm(y~x+I(x^2)))     resettest(reg2,power=2,type="regressor")          set.seed(12344321)     x <- seq(1,100)     y <- x + 0.001*rnorm(100)          reg3 <- lm(y~x)     summary(reg3)          # Should match and do, even though almost the same as reg1 example     summary(lm(y~x+I(x^2)))     resettest(reg3,power=2,type="regressor")      In this code sample, the p-value for the x^2 term should match the p-value for the `resttest` in each of the three pairs in the code. What actually happens is that they fail to match in examples reg1 and reg2---two cases where the regression fits perfectly. In example reg3---which is almost identical to example reg1, differing only by the addition of a tiny bit of noise---they do match. This does not have a lot of practical significance since real world regressions basically never fit perfectly, but it is disconcerting if you are just playing around with the function to convince yourself that you understand what it does. Does anyone understand why this seemingly undesirable behavior occurs or if I am doing something wrong in my use of `resttest`? Maybe the test statistic ends up being calculated as (roundoff error)/(roundoff error)?


title: Probability that a random number will be the largest from sets

text: Let's say that I have three ranges of numbers: $$ A: [0.15, 0.26]\\\ B: [0.25, 0.34]\\\ C: [0.20, 0.35]$$ If a number was randomly picked from the range of each set, what is the probability that the number from set $A$ will be the largest? From set $B$? etc?


title: How to use extracted data from factor analysis in subsequent 2x2 factorial design?

text: I’m conducting a 2x2 between subjects factorial design for my dissertation. **IV’s** :   1. advertising appeal (two levels: rational appeal vs. emotional appeal)    2. media type (two levels: TV vs. Internet).  **DV** : cognitive (thinking) and affective (emotional) customer response **Covariates** :   1. subjects’ product knowledge    2. involvement with the product    3. yearly spending on high-involvement products    4. general attitudes toward the advertisement    5. general attitude toward the brand.  **Hypotheses** :   1. A combination of different advertising media and appeals (TV/emotional, TV/rational, Internet/emotional, Internet/rational) differently affect consumers' cognitive and emotional responses.   2. Rational Internet ad for high-involvement product is more effective in generating cognitive response compared to rational TV ad for the same category of product.   3. Emotional TV ad for high-involvement product is more effective in generating affective response compared to emotional Internet ad for the same category of product. Factor analyzed data that has been collected from four treatment conditions (TV / emotional, TV / rational, Internet / emotional, Internet / rational) need to be used in further causal analysis to find out interaction effects (or testing hypothesis 2 and 3). **Question** : How can the extracted factors (4 emotional and 4 cognitive) be used in 2x2 MANCOVA analysis?


title: Finite output in logistic regression with divergent `glm` algorithm

text: let `m` be my matrix of data                     x_i y_i      [1,] 0.0   0      [2,] 0.0   0      [3,] 0.0   0      [4,] 0.0   0      [5,] 0.1   0      [6,] 0.2   0      [7,] 0.3   0      [8,] 0.4   0      [9,] 0.5   0     [10,] 0.6   0     [11,] 0.0   1     [12,] 0.0   1     [13,] 0.0   1     [14,] 0.9   1     [15,] 1.0   1      My aim is to study the logistic regression `y~x`, where the covariate `x` has observations `m[,1]` and similarly for `y`. Please note that we have no complete separation in the data _but_ the "anomalous" entries in rows `m[11,], m[12,]` and `m[13,]` all correspond to observations with `x_i=0`. I expect `glm` to diverge as the likelihood function reaches no maximum in the ray $k\beta$, for $k\rightarrow \infty$ and $\beta=(-0.7,1)$. Using `glm` with 1 iteration I get the output                 Coefficients:                 Estimate Std. Error z value Pr(>|z|)     (Intercept)  -1.2552     0.7648  -1.641    0.101     x             1.6671     1.7961   0.928    0.353          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 19.095  on 14  degrees of freedom     Residual deviance: 18.275  on 13  degrees of freedom     AIC: 22.275          Number of Fisher Scoring iterations: 1      with an error message (the algorithm does not converge). Moreover, with the default number of iterations `(=25)` the output is               Coefficients:                 Estimate Std. Error z value Pr(>|z|)     (Intercept)  -1.1257     0.7552  -1.491    0.136     x             1.4990     1.6486   0.909    0.363          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 19.095  on 14  degrees of freedom     Residual deviance: 18.246  on 13  degrees of freedom     AIC: 22.246          Number of Fisher Scoring iterations: 4      and no error warning. I see a contradiction; even in presence of 1 iteration the algorithm does not converge but the output is "finite" (I have not explicitly computed the inverse of the Hessian of the likelihood function, unfortunately). Moreover, with 25 iterations the warning message disappears and the output is still finite. What do you think about this situation? Is it possible that `glm` stops automatically after the first iteration? Thank you, Avitus


title: ANOVA sum of squares between groups

text: I'm trying to get an intuitive understanding of why the sum of squares between groups needs to be multiplied by the number of observations within each group. Using the `iris` dataset in `R` as an example, the sum of squares needs to be multiplied by 50:               meanOverall = mean(iris$Sepal.Length)          meanSetosa = mean(iris$Sepal.Length[iris$Species=='setosa'])          meanVersicolor = mean(iris$Sepal.Length[iris$Species=='versicolor'])          meanVirginica = mean(iris$Sepal.Length[iris$Species=='virginica'])           sumsquaresBetweenGroups = sum((meanSetosa - meanOverall) * (meanSetosa - meanOverall) + (meanVersicolor - meanOverall) * (meanVersicolor - meanOverall) +      (meanVirginica - meanOverall) * (meanVirginica - meanOverall)) * 50      Why do you need to multiply the between groups sum of squares by the number of observations in each group? In addition, what if $n$ is different in each group. What would you do if there were 45, 55 and 60 observations in each of the three groups of `Species` in the `iris` dataset?


title: Evaluate statistical significance of difference between outcomes of tests

text: I have 400 patients and I apply 4 tests to all of them. This gives me a number of positive and negative for each test. I want to compare the tests 2 by 2 to find out if there is any statistically significant difference between them.               Test  Pos   Neg     1     100   300     2     150   250     3     160   240     4     200   200      What statistical test I can use to determine this? I'm very new to the field so I appreciate any pointers.


title: How to display error bars for cross-over (paired) experiments

text: The following scenario has become the Most-FAQ in the trio of investigator (I), reviewer/editor (R, not related to CRAN) and me (M) as plot creator. We can assume that (R) is the typical medical big boss reviewer, who only knows that each plot must have error bar, otherwise it is wrong. When a statistical reviewer is involved, problems are much less critical. **Scenario** In a typical pharmacological cross-over study, two drugs A and B are tested for their effect on glucose level. Each patient is tested twice in random order and under the assumption of no carry-over. The primary endpoint is the difference between glucose (B-A), and we assume that a paired t-test is adequate. (I) wants a plot that shows the absolute glucose levels in both cases. He fears (R)'s desire for error bars, and asks for standard errors in bar graphs. Let's not start the bar graph war here ._) ![bar graphs and SE glucose for two treatments](http://i.stack.imgur.com/FU7Rz.png) (I): That cannot be true. The bars overlap, and we have p=0.03? That's not what I have learned in high school. (M): We have a paired design here. The requested error bars are totally irrelevant, what counts is the SE/CI of the paired differences, which are not shown in the plot. If I had a choice and there were not too many data, I would prefer the following plot ![Joined lines show the pairing, dots the original values](http://i.stack.imgur.com/iK9bI.png) **Added 1:** This is the parallel coordinate plot mentioned in several responses (M): The lines show the pairing, and most lines go up, and that's the right impression, because the slope is what counts (ok, this is categorical, but nevertheless). (I): That picture is confusing. Nobody understands it, and it has no error bars (R is lurking). (M): We could also add another plot that shows the relevant confidence interval of the difference. The distance from the zero-line gives an impression of the effect size. (I): Nobody does it (R): And it wastes precious trees (M): (As a good German): Yes, point on the trees is taken. But I nevertheless use this (and never get it published) when we have multiple treatments and multiple contrasts. ![confidence interval of the difference](http://i.stack.imgur.com/bziqw.png) **Any suggestions**? The R-Code is below, if you want to create a plot.               # Graphics for Crossover experiments     library(ggplot2)     library(plyr)     theme_set(theme_bw()+theme(panel.margin=grid::unit(0,"lines")))     n = 20     effect = 5      set.seed(4711)     glu0 = rnorm(n,120,30)     glu1 = glu0 + rnorm(n,effect,7)     dt = data.frame(patient = rep(paste0("P",10:(9+n))),                                   treatment = rep(c("A","B"), each=n),glucose = c(glu0,glu1))          dt1 = ddply(dt,.(treatment), function(x){       data.frame(glucose = mean(x$glucose), se = sqrt(var(x$glucose)/nrow(x)) )})          tt = t.test(glucose~treatment,paired=TRUE,data=dt,conf.int=TRUE)     dt2 = data.frame(diff = -tt$estimate,low=-tt$conf.int[2], up=-tt$conf.int[1])         p = paste("p =",signif(tt$p.value,2))          png(height=300,width=300)     ggplot(dt1, aes(x=treatment, y=glucose, fill=treatment))+             geom_bar(stat="identity")+         geom_errorbar(aes(ymin=glucose-se, ymax=glucose+se),size=1., width=0.3)+       geom_text(aes(1.5,150),label=p,size=6)          ggplot(dt,aes(x=treatment,y=glucose, group=patient))+ylim(0,190)+       geom_line()+geom_point(size=4.5)+       geom_text(aes(1.5,60),label=p,size=6)          ggplot(dt2,aes(x="",y=diff))+       geom_errorbar(aes(ymin=low,ymax=up),size=1.5,width=0.2)+        geom_text(aes(1,-0.8),label=p,size=6)+       ylab("95% CI of difference glucose B-A")+  ylim(-10,10)+       theme(panel.border=element_blank(), panel.grid.major.x=element_blank(),              panel.grid.major.y=element_line(size=1,colour="grey88"))          dev.off()


title: Cell chi square test

text: I need to do a cell chi square test, but googling availeth me not (I suspect because of all the pages about the cells in chi square discussions generally). How does one calculate this? Any good sources on it, or on it in R specifically?


title: Generator and/or interpreter of probability models

text: Even simple models in probability theory can be quite confusing (e.g. drawing something with or without returning it, conditional probabilities, one and only one outcome, at least one outcome asoasf...) In other areas of mathematical modelling (e.g. differential equations) you often have powerful and easy to use modelling tools - I have never seen this for stochastic models. Do you know of any   (a) generator with which you could in a (e.g.) point-and-click style setup certain probability models and the generator gives you the respective formulas.   (b) an interpreter which gives you possible interpretations of probability formulas in basic (e.g.) urn or coin-toss-terms. Thank you!


title: Large n small p regression - Machine Learning

text: In the area of machine learning, most of the algorithms are intended for small n large p problems. I am familiar with the statistical techniques of PCA, etc but was wondering what algorithms are available in the machine learning area for similar analysis work. The aim is to determine the effect of a limited number of predictor variables on the outcome using robust methods such as pca in the ML domain, Thanks in advance.


title: A measure of the "flatness" of log likelihood at the MLE

text: Are there meaningful ways to quantify how "flat" the log likelihood function is around the MLE when the parameter has more than one dimension? In particular is the determinant of the Hessian a reasonable measure?


title: Interpreting Intercept when doing logistic regression with categorical data in R

text: In `R` I have a categorical variable that I performed logistic regression on and got the following result:               glm(formula = mortality ~ SMOKE, family = binomial, data = c.data)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -0.2155  -0.2155  -0.2155  -0.1860   2.8515            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -4.0483     0.3189 -12.694   <2e-16 ***     SMOKEN        0.2968     0.3559   0.834    0.404         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for binomial family taken to be 1)              Null deviance: 492.45  on 2369  degrees of freedom     Residual deviance: 491.72  on 2368  degrees of freedom     AIC: 495.72          Number of Fisher Scoring iterations: 6      Is the value for the intercept the same as `SMOKEY` (has a history of smoking)?


title: Plotting Group-Differentiated Coefficients

text: Using lme package in R to simulate bivariate regression by interacting all IVs with a dummy-coded groupid variable, I've used coefplot2, but would like plots grouped by IV to visually represent the differential association of x and y by group. I can do the grouping by hand using ggplot2 if I have to, but I'm not 100% sure about which points to plot, especially for the min and max of the confidence bands. With the example Y = b0 + b1(X) + b2(Z) + b3(XZ) + e, plots for Group 1 (Z=0) will be centered at b1 with min and max for the confidence band formulated as b1+-1.96*SE. Plots for Group 2 (Z=1) will be centered at b1+b3, correct? If that's right, how do I formulate the min and max for the confidence band for Group 2? Would I just use the SE associated with b3?


title: Why are lower p-values not more evidence against the null?

text: Johansson (2011) in "Hail the impossible: p-values, evidence, and likelihood" (here is also link to the journal) states that lower p-values are often considered as stronger evidence against the null. Johansson implies that people would consider evidence against the null to be stronger if their statistical test outputted a p-value of 0.01, than if their statistical test outputted a p-value of 0.45. Johansson lists four reasons why the p-value cannot be used as evidence against the null: >   1. p is uniformly distributed under the null hypothesis and can therefore > never indicate evidence for the null. >   2. p is conditioned solely on the null hypothesis and is therefore > unsuited to quantify evidence, because evidence is always relative in the > sense of being evidence for or against a hypothesis relative to another > hypothesis. >   3. p designates probability of obtaining evidence (given the null), rather > than strength of evidence. >   4. p depends on unobserved data and subjective intentions and therefore > implies, given the evidential interpretation, that the evidential strength > of observed data depends on things that did not happen and subjective > intentions. > Unfortunately I cannot get an intuitive understanding from Johansson's article. To me a p-value of 0.01 indicates there is less chance the null is true, than a p-value of 0.45. Why are lower p-values not stronger evidence against null?


title: How to carry out multiple post-hoc chi-square tests on a 2 X 3 table?

text: My data set is comprised of either total mortality or survival of an organism at three site types, inshore, midchannel and offshore. The numbers in the table below represent the number of sites.                             100% Mortality            100% Survival     Inshore             30                       31      Midchannel          10                       20      Offshore             1                       10      I would like to know if the # of sites where 100% mortality occurred is significant based on site type. If I run a 2 x 3 chi-square, I get a significant result. Is there a post-hoc pairwise comparison that I can run or should I actually be using a logistical ANOVA or regression with binomial distribution? Thanks!


title: The estimation of variance in simulation

text: I am confused about the calculation of variance in a regression simulation. I can get the sample of the estimator from the simulation. So there are two ways to get the variance of this estimator: one is just the sample standard error while the other is the estimator variance that we infer before (e.g. the variance formula of OLS estimator). Sometimes the results above are not the same. What is the difference between these two calculations of variance? And which one shall I choose if I do a simulation?


title: Multiple Level Regression FDI

text: I am trying to work out an econometric model that has multiple layers, and I am not sure if it is feasible. I have 8 home countries that are investing in 3 host countries, and I am looking at the evolution of investment over a period of 10 years (2002-2011). As my independent variable I have the FDI inflow into the the 3 host countries from the 8 source countries for each year. As my independent variables, I have for example GDP/capita of host country i over the 10 year period. Now the same GDP/capita, will apply to the 8 source countries in my panel. How do I deal with these repeated observations in my panel. I am in dire need of help. I need to know know if such a panel is possible and what type of regression I need to apply. Thank you for your time!


title: Repair flawed experiment design

text: I seem to have uncovered a methodological flaw in a published large study dataset that originally took about 30 years to collect and collate. The experiment had a sequence of stimuli in a globally unique randomized order — i.e. each interview was conducted in the same order, but the order was randomized once as part of the experiment design. By analyzing the data, I have discovered that the first stimulus behaves differently from the rest: adjacent stimuli correlate pretty well (Absolute value of Pearson r between 0.3 and 0.85, with p-values far under 0.05); except for the first stimulus in the randomized order, which correlates with its neighbors far less (Pearson r at around -0.07, p ~ 1e-25). Before writing these findings up and publishing them I would like to be able to suggest methods to still use the collected data for further analysis, which leads me to my question: What strategies might be worth pursuing for extracting reliable conclusions from a dataset suffering from this type of sequential “burn-in” issues with the data collection method? Is there a good source I can read on handling data from flawed experimental design?


title: When we are dealing with instrumental variables is this correct?

text: If I'm to find an instrumental variable for an equation, am I getting this idea right? I have this regression: $$\text{stndfnl}=\beta_0+\beta_1\text{atndrt}e+\beta_2\text{priGPA}+\beta_3\text{ACT}+\beta_4\text{priGPA}*\text{atndrte}+u$$ where $\text{stndfnl}$ is the standardized outcome on a final exam, $\text{atndrte}$ is attendance rate, $\text{priGPA}$ is prior college GPA, and $\text{ACT}$ (if not obvious given the context) is the score of the standardized test, the ACT. I'm told that if $\text{atndrte}$ is correlated with $u$, then, in general, so is $\text{priGPA*atndrte}$. So what might make a good instrumental variable for $\text{priGPA*atndrte}$? Do you think, say, the mother's education level ($\text{motheduc}$) would be a good IV? It could have a significant effect on $\text{priGPA*atndrte}$ and it could also be uncorrelated with the error term in the main regression. Is this the correct way to approach this?


title: Jaccard Indexes and PCA

text: Does it make sense to use a matrix that is made up of jaccard indexes instead of a covariance matrix and perform principal component analysis on that?


title: Learning material about time series

text: > **Possible Duplicate:**   >  Books for self-studying time series analysis? I am new to time series modelling altogether. But I am aware about regression modelling and some data mining algorithms like decision trees. I want to learn time series from scratch. My background is Mathematics. Please suggest me any good book/material/web resource for step by step learning with case studies. Thanks


title: Graph Theory - Creating an Index of Familiarity, Given Trade Frequency Counts

text: **Set Up** I'm hoping to create an "index of familiarity" between traders on a barter market. I have data from a peer-to-peer barter market (i.e. people come in with their wares, and can trade with other individuals). The dataset records how many times each trader traded with each of all other traders who on the market (i.e. a symmetric n-by-n matrix of trade frequency counts). Here's a sample dataset (.CVS & .RData). I also have access to a similar symmetric n-by-n matrix with _VALUES_ of items exchanged (which I suspect will be more telling). **Middlemen** : I also need to be careful about 'middlemen'. These traders frequently serve as middlemen (duh), but who are not really 'familiar' with their counter-parties. I have a list of who these traders are (in my sample dataset, trader 5 is the only middlemen), and so hope to control for this. At least I hope their existence doesn't throw off index values for other traders. I suspect I should exclude middlemen from my matrix. **Example Clusters** : If you look at my example dataset (see CSV or RData links above), it should be clear that traders t1, t2, t3, and t4 form a familiarity cluster (I'll call Cluster 1) and so they should all have values that suggest they are close to one-another, and distant from the t6, t7 and t8. Also, traders t6, t7, and t8 cluster into Cluster 2. **Also:** As an added layer of nuance I hope to capture, note that t2 only trades with the middleman and with one counterparty in her cluster (Cluster 1). Thus I'd like to familiarity index calculated to suggest t2 is far "more familiar" with t1, t3 and t4 (folk in Cluster 1) than with t6, t7 and t8 (traders in Cluster 2). The output I am aiming for is an n-by-n symmetric matrix of familiarity index values between each traders on the market. I have no experience with graph theory or existing tools, but am happy to do what it takes to learn. I have this weird feeling that my question is perhaps a slight tweak on a well understood problem - one with lots of tools available to me. Thus I'm hoping there's a method you can recommend with an existing implementation. **My Background** : I mostly program in R (I'd prefer to work with cran packages), but can work with Python and Java if that's going to be preferable. Thanks for any advice you might have! **My Best Stab at a Familiarity Function** (in R code)               FamiliarityIndex <- function(data){       # Takes input dataset of symmetric n-by-n matrix of trade frequency counts       # return a 'familiarity index'       # How to interpret: how familair Trader COLUMN is with this Trader number ROW       #    e.g., if col 4 row 1 is large, that's how familiar trader 4 is with trader 1        #            (it might differ for col 1 row 4!)            #rescale counts       newData = data       for (i in 1:nrow(data)){         colSum = sum(data[,i], na.rm=T)         newData[,i] = data[,i] / colSum       }            #        data <- newData       for (col in 1:ncol(data)){         for (row in 1:ncol(data)){           if (is.na(newData[row,col])){} else {             newData[row,col] = newData[row,col] +                                 sum(data[row,1:(ncol(data))] * data[,col], na.rm=T)                }              }       }            return(newData)          }


title: Procedure to select transformation on linear regression

text: I need to automate the transformation on some linear regression models. There is only one predictor in this case. Sometimes i get a good model with the original variables, sometimes i need to log the predictor, and in some cases log both sides. I'm using R, so what kind of tests/packages can i use to automate this? I'm using Pearson correlation now, but i'm not sure if it makes sense. thanks! PS: This may look a duplicate question, but i couldn't find yet the methodology to apply.


title: How can I check if two distributions are independent

text: A telephone salesman calls 3 people in a month and sells a product with probability 30%. In each call he sells either 0 or 1 products. He calls three people this month. If he manages to sell the product in his first or second call he gets 100 euro per successful call and 200 if he sells a product in his last call. Let $X$ be number of products he sells and $Y$ the sum of euros he sells. Are $X$ and $Y$ independent? I managed to compute total distributions and marginal distributions of $X$ and $Y$ but how can I compute $\mathbb{P}(X=x\cap Y=y)?$


title: Finding choice probabilities by using utility with logit and probit models

text: I am using a formula to calculate the utility, which is as follows:               v_{ij} = 1 - x*beta + delta_i + e_{ij}          delta_i ~ N(0,phi^2)     e_ij ~ N(0,sigma^2)      v_{ij} is the utility of state ij, x is the vector of dummy variables, beta is the vector of regression coefficients, delta is the i-th random error and e_{ij} is the usual error term. So, let's say there are two alternatives to make a choice from and their utility functions are                   v_{i1} = 1 - x_{i1}*beta + delta_i + e_{i1}       and                   v_{i2} = 1 - x_{i2}*beta + delta_i + e_{i2} .      I am wondering how I could calculate the probability of choosing the first option. I know I need to use this               P(Choosing the first) = P( v_{i1} >  v_{i2})      = ( 1 -x_{i1}*beta + delta_i + e_{i1}  > 1 - x_{i2}*beta+ delta_i + e_{i2} )      Can anyone please help me find this probability by using both logit and probit models.


title: Multidimensional d prime

text: If you have a one dimensional space, like you have one variable, and two sets of observations, and want a measure of separation between the two sets of observations you can calculate a d prime like $d=(\mu1-\mu2)/\sqrt{(s_1^2+s_2^2)/2}$ basically distance between means over the square root of sum of variances over two,   I am not sure, but I guess there is an assumption of normality and positive means here, My question is, what if you are in a multi dimensional (or at least two dimensional) space, so you have more than one variable that defines each observation? How do you calculate the 2D or ND d prime?


title: How to combine sample means and sample variances?

text: Given a set of data, assumed normally distributed from 1 to 19. And two subsets of that data, one with the numbers below 10, and one with the numbers above:               subset1 = [1 2 3 4 5 6 7 8 9 9 8 7 6 5 4 3 2 1];     subset2 = [15 16 17 18 19 19 18 17 16 15 15 16 17 18 19 18 17 16];     data    = [subset1 subset2];      We can calculate the means ($\mu_1,\mu_2,\mu_{tot}$) and SD's ($\sigma_1,\sigma_2,\sigma_{tot}$) for all three vectors above and if we assume the distributions to be Gaussian, we get the following graph: ![Gaussian distributions of the data sets](http://i.imgur.com/FhI7fq0.png) But now let's assume the data set gets lost before I was able to calculate the total mean and variance (I know, it's ridiculous). Only being left with $\mu_1,\mu_2$ , $\sigma_1,\sigma_2$ and the lengths of the subsets $n_1,n_2$. Could I reconstruct the total pdf with the parameters of the subsets? The mean is rather easy: $\hat{\mu}=\frac{n_1\mu_1+n_2\mu_2}{n_1+n_2}=\mu_{tot}$ But what about the variance $\sigma_{tot}$? I read about pooled variance, but this is not depending on the means of the subsets and can't therefore possibly be true, since the spacing of the two "subgaussians" is obviously of influence. [EDIT] As a matter of fact, there is! How to calculate pooled variance of two groups given known group variances, means, and sample sizes? I apologize. Follow-up question: Does the same hold for multivariate pdf's?


title: How to compute the density of the product of a discrete and a continuous variable?

text: Let's say $X$ and $Y$ are continuous and discrete random variables, respectively, with $f(x)$ and $g(y)$ being a probability density for $X$ and probability mass for $Y$, respectively. Can I say that $Z=XY$ has density function equal to $$\sum_{y} f(z/y)g(y)$$ and can I generalize this for all pairs of random variables?


title: Properties of AR(p) model to be invertible

text: Suppose we are working with an $AR(p)$ model. We assume that $E( \epsilon^2 ) < \infty$. We know that an $AR(P)$ model is invertible $\text{iff } \phi(z) \neq 0$ for $|z| \le 1$. But why do we have to make this assumption? Why can we not have a zero inside the unit circle, but outside the unit circle it's fine?


title: What statistical model for price optimisation

text: I'm not an expert so please be understanding if I don't formalize the problem in a professional/correct way. I have a data set with more than 160.000 offers that were made to the users of a website. Each offer has a binary outcome, offer accepted or offer rejected. The data is:   * Y purchase status, binary outcome variable which can be 0 or 1 if the user accepted or not the offer   * X set of attributes available for each user (such as gender, age group, income group, relationship status, number of kids etc)   * DR discount rate of the offer, which is positively correlated with the outcome variable as you can see in the image below ![enter image description here](http://i.stack.imgur.com/gY38f.jpg) What I'd like to do: I'd like to design a model to define an optimum offer for each user based on the user's attributes. Basically being able to offer a lower discount rate to a user with a certain set of demographics. The assumption is that different users are more or less sensible to the discount rate. The DR can variates in an interval which we can consider given, for simplicity. How would you tackle this problem? An initial idea was to proceed this way: doing a regression on the entire population with Y as the output variable and X,DR as dependent variables. After that apply the coefficients on single cases to get an idea of the single user respect to the population. After that calculate some sort of score that will then be used to calculate the DR to apply to the user. **Am I completely going in the wrong direction?** Could you please give me a feedback and eventually point me to some relevant papers on the argument? Thanks for your help, Giacomo


title: Alternative name for dispersion or spread

text: The determinant of the covariance matrix gives us the square of the volume of a hyperrectangle aligned with the distribution whose side lengths are determined by the variances of the distribution in the eigenspace. In the univariate case this boils down to variance. But what should we call such a measure? What does this measure mean? On the one hand we could say it measures the spreadedness of the distribution, or simply its spread. Synonymical names used in literature are 'dispersion', 'variability' and 'scatter'. My problem with these names is that they are too ambiguous. 'Spread' could also refer to uniformity. I need to disambiguate between uniformity and some other concept for which I am trying to find the right word. What word would be appropriate?


title: Maximize Sortino Ratio in R

text: I'm trying to find portfolio weights, which will maximize my Sortino Ratio. I have weekly log-returns on 16 stocks. I've used all day yesterday to install and implement the following example on my own data: https://r-forge.r-project.org/scm/viewvc.php/pkg/PortfolioAnalytics/demo/sortino.R?view=markup&root=returnanalytics&pathrev=1765 The code runs fine, but my results are way off. So I figured I could write my own function to calculate the Sortino Ratio and then use some sort of optimization function in R such as optim. So I've written a function to calculate the Sortino Ratio and I'm then using optim to maximize it:               sorsol <- optim(rep(0,16), method="BFGS",      fn = Sortino, dataset = afkastweekly,     control=list(fnscale=-1))      I need to impose a constraint which says the weights much sum to one, but having used many hours to try this in ConstrOpt witout a solution, I figured I could rescale my weights after the maximation, but is this approach incorrect?               sorsol$par/sum(sorsol$par)      And if so why? I believe it is, since I have another portfolio (The tangency portfolio), which yields a higher Sortino Ratio, so it must be incorrect, but I can't see why. If it is incorrect, can anyone point me in a helpfull direction as to maximize the Sortino Ratio? Below you will find a dput of (some) of my data:               > dput(head(afkastweekly))     structure(c(0.0496716097972385, 0.0196298560448849, 0.0150056406178702,      -0.0390505148718727, 0.0326470778366659, 0.0460332318059749,      -0.0125472160520887, 0.0396091380950461, -0.0146701897477941,      0, 0.0195128142235816, 0.00961545869944214, 0.10227924687079,      -0.00464064410673171, -0.0187799192715268, 0.0417697208491958,      -0.0137309788880824, 0.0450532353179636, 0.023311078868447, 0.0212773984472854,      0.0030030052597696, 0.01340303012734, -0.0300322870988752, -0.0122700925918142,      0.0747800151897677, 0.0108578699720492, 0.0297894368021074, 0.0668942348300301,      -0.0440953198652241, 0.00612872194137326, -0.0528105809038451,      0.0299167610379927, -0.00845671001822357, 0.0126583968719234,      0.105360515657826, 0.0149815536156166, -0.0504708516508101, 0.00235118644773902,      0.0277798072067919, -0.0230968977290251, -0.00703409592550575,      -0.0214010936503484, 0.00775197680431772, 0.0415938972988372,      -0.0168071183163812, 0.046006273008643, 0.00716848947861326,      0, -0.0294666592968058, 0.0258298703198072, -0.0371094812328154,      -0.001889903852625, 0.0131695147656332, -0.018870660553068, 0.0514031784599656,      0.016563525671673, 0.0203259030140401, -0.0473827044015529, 0.0331980694095968,      0.0497615095590636, 0.0250639686632166, 0.110052461541528, -0.017897569457543,      -0.0274616799413878, 0.0229367853430982, 0.0378243086673389,      0.0616935690053397, -0.00858374369139181, 0.0504308536268923,      0.123232640423948, -0.0989399478549036, 0.0314986670593704, -0.0451204352804702,      0.0741079721537226, 0, -0.0588405000229342, 0.00754720563538314,      0.00749067172915741, 0.0112571545246345, 0.00186393343806257,      0.0184507079131166, -0.0353525187157198, 0.0150378773645405,      -0.00186741417479563, -0.0141843740698722, 0.010656907162943,      0.00352746690692918, -0.00352746690692918, 0.0278766136011539,      -0.0243491466942247, 0, 0.0137301928119018, 0.00227014853453911,      0.0721624462482309, -0.00634922767865875, 0.0375043954584542), .Dim = c(6L,      16L), .Dimnames = list(NULL, c("MRSKA", "WDH", "VWS", "TOP",      "TDC", "NOVOB", "NRDEA", "COLOB", "CARLB", "MAERSKB", "GN", "GENMAB",      "FLS", "DSV", "DANSKE", "NVZY")), index = structure(c(12426,      12433, 12440, 12447, 12454, 12461), class = "Date"), class = "zoo")      Here's the function that calculates my Sortino Ratio:               Sortino <- function(weights, dataset, MAR = 0.0319, scale = 52)                 {                     pafkast <- NULL                     for (i in 1:dim(dataset)[1])                         {                             pafkast <- rbind(pafkast, as.vector(weights) %*% as.vector(dataset[i,]))                         }                     Er <- ExpectedReturns(pafkast) * scale                     DD <- DownsideDeviation(pafkast, MAR = MAR) * sqrt(scale)                     SR <- (Er-MAR) / DD                     return(SR)                 }      The function DownsideDeviation is from the package PerformanceAnalytics. If anything else is needed, please let me know.


title: Multi-class Confusion Matrix to Binary confusion matrix

text: i know the main concepts of data/text mining but i used them mainly in binary classification problems (just two classes). i am now dealing with a problem with 8 classes and i am atruggling how to perform evaluation measurement, like precision, recall and others. My question is, can i convert a multi class matrix to a binary conf matrix with TP, FP, TN FN values and then calculare the measures i want? if not, can you give me some help on how i perform evaluation os a classifier when i got a confusison matrix like: https://docs.google.com/spreadsheets/d/1yhGqO5KK_B9U6HMmxIQyincE3NNht- zUgdMQiPOLeu8/edit?usp=sharing Feel free to edit document.


title: Categorical Binary Variable

text: I'm attempting to draw an inference of accuracy be it by confidence interval or response rate, except typical methods don't really suite my problem for the reason that the entire population was surveyed for a yes or no vote. The 'sample' of data which is categorical binary variables, then is simply the proportion of the population that responded, which was 63% n=47. Accordingly, the non response isn't represented as it is in typical literature, whereby non response is framed in the proportion of response within a sample taken from the population and not where the entire population is surveyed. Further, central limit theorem seems to be empirically more robust regarding continuous data and particular parameters such as mean. Which is not to say that the convergence and attractor functions will not be relevant here. I'm assuming both ways of framing the problem are empirically descriptive, but I am struggling to frame the issue with technical precision. Any input would be great.


title: Using Wilcoxon and Cliff's delta for comparing two sets of monitored variables?

text: I want to compare the following example vectors:               a: [0, 1, 2, 3, 4]     b: [4, 3, 2, 1, 0]      `a[i]` represents a monitored variable `i` during execution `a` (likewise for `b`) The things I would like to find out are:   1. Are the pairs significantly different?   2. How different are they?   3. Which ones are different from each other (pairwise comparisons)? For (1) and (2) I am thinking about the Wilcoxon signed rank test and Cliff's delta effect size measurement. (3) can be done using simple programming. So, my questions are:   1. Are these the right test and method for this problem?   2. When I perform Wilcoxon in Python for vectors `a` and `b`, I get `(5.0, 1.0)` as output. I do not understand the p-value of 1.0 as I consider the sets to be very different. Also, what is the 5.0 value? Is it the W value? When I manually calculate this, I get 0.0.... **Edit** To clarify what 'different' means to me: I am working on a method to compare usage profiles monitored during test execution of software. So, an example profile may be a:               a: [func1 = 0sec, func2 = 1sec, func3 = 2sec, func4 = 3sec, func5 = 4sec]     b: [func1 = 0sec, func2 = 3sec, func3 = 2sec, func4 = 3sec, func5 = 4sec]      with sec being total seconds execution time in that function during run a or b of a specific test, but with b ran on (slightly) changed code (= a new version). I would like to compare this profile to another profile to verify whether execution time has changed significantly. From the answers I understand that it is difficult to say something about significance for only two executions. Would it be better if I have 10 executions of a and 10 of b, calculate mean and stdev for all functions and compare those using Wilcoxon?


title: Label on the y-axis in a normalised histogram

text: If you have a histogram with frequency on the y-axis and bins for different ranges of values on the x-axis, then it is reasonable that the label on the y-axis should be frequency. But if these frequencies are normalised, what is the correct y-label? Normalised frequency? Quota? I use Python's matplotlib:               import matplotlib.pyplot as plt     l = [3,3,3,2,1,4,4,5,5,5,5,5]     plt.hist(l,normed=True)     plt.show()      It seems that a bin can reach above 1.0 which I don't really want. I want it to be more like proportion or fraction; the height of the bins should sum up to one.


title: Fourier/trigonometric interpolation

text: ## Background In a paper from Epstein (1991): On obtaining daily climatological values from monthly means, the formulation and an algorithm for calculating Fourier interpolation for periodical and even-spaced values are given. In the paper, the goal is to **obtain daily values from monthly means** by interpolation. In short, it is assumed that unknown daily values can be represented by the sum of harmonic components: $$ y(t) = a_{0} + \sum_{j}\left[a_{j}\,\cos(2\pi jt/12)+b_{j}\,\sin(2\pi jt/12)\right] $$ In the paper $t$ (time) is expressed in months. After some derviation, it is shown that the terms can be calculated by: $$ \begin{align} a_{0} &= \sum_{T}Y_{T}/12 \\\ a_{j} &= \left[ (\pi j/12)/\sin(\pi j/12)\right] \times \sum_{T}\left[Y_{T}\,\cos(2\pi jT/12)/6 \right]~~~~~~~j=1,\ldots, 5 \\\ b_{j} &= \left[ (\pi j/12)/\sin(\pi j/12)\right] \times \sum_{T}\left[Y_{T}\,\sin(2\pi jT/12)/6 \right]~~~~~~~j=1,\ldots, 5 \\\ a_{6} &= \left[ (\pi j/12)/\sin(\pi j/12)\right]\times \sum_{T}\left[Y_{T}\cos(\pi T)/12\right] \\\ b_{6} &= 0 \end{align} $$ Where $Y_{T}$ denote the monthly means and $T$ the month. Harzallah (1995) summarizes this aproach as follows: "The interpolation is carried out by adding zeros to the spectral coefficients of data and by performing an inverse Fourier transform to the resulting extended coefficients. The method is equivalent to applying a rectangular filter to Fourier coefficients." * * * ## Questions My goal is to use the above methodology for interpolation of **weekly means to obtain daily data** (see my previous question). In summary, I have 835 weekly means of count data (see the example dataset at the bottom of the question). There are quite a few things that I don't understand before I can apply the approach outlined above:   1. How would the formulas have to be changed for my situation (weekly instead of monthly values)?   2. How could the time $t$ be expressed? I assumed $t/835$ (or $t/n$ with $n$ data points in general), is that correct?   3. Why does the author calculate 7 terms (i.e. $0\leq j \leq 6$)? How many terms would I have to consider?   4. I understand that the question can probably be solved by using a regression approach and using the predictions for interpolation (thanks to Nick). Still, some things are unclear to me: How many terms of harmonics should be included in the regression? And what period should I take? How can the regression be done to ensure that the weekly means are preserved (as I don't want an exact harmonic fit to the data)? Using the regression approach (which is also explained in this paper), I managed to get an exact harmonic fit to the data (the $j$ in my example would run through $1, \ldots, 417$, so I fitted 417 terms). **How can this approach be modified -$~$if possible$~$- to achieve the conservation of the weekly means?** Maybe by applying correction factors to each regression term? The plot of the exact harmonic fit is: ![Exact harmonic fit](http://i.stack.imgur.com/7XuxU.png) **EDIT** Using the signal package and the `interp1` function, here's what I've managed to do using the example data set from below (many thanks to @noumenal). I use `q=7` as we have weekly data:               # Set up the time scale          daily.ts <- seq(from=as.Date("1995-01-01"), to=as.Date("2010-12-31"), by="day")          # Set up data frame           ts.frame <- data.frame(daily.ts=daily.ts, wdayno=as.POSIXlt(daily.ts)$wday,                            yearday = 1:5844,                            no.influ.cases=NA)          # Add the data from the example dataset called "my.dat"          ts.frame$no.influ.cases[ts.frame$wdayno==3] <- my.dat$case          # Interpolation          case.interp1 <- interp1(x=ts.frame$yearday[!is.na(ts.frame$no.influ.case)],y=(ts.frame$no.influ.cases[!is.na(ts.frame$no.influ.case)]),xi=ts.frame$yearday, method = c("cubic"))          # Plot subset for better interpretation     par(bg="white", cex=1.2, las=1)     plot((ts.frame$no.influ.cases)~ts.frame$yearday, pch=20,          col=grey(0.4),          cex=1, las=1,xlim=c(0,400), xlab="Day", ylab="Influenza cases")     lines(case.interp1, col="steelblue", lwd=1)      ![Cubicinterpo](http://i.stack.imgur.com/R1FE8.png) There are two issues here:   1. The curve seem to fit "too good": it goes through every point    2. The weekly means are not conserved **Example dataset**               structure(list(date = structure(c(9134, 9141, 9148, 9155, 9162,      9169, 9176, 9183, 9190, 9197, 9204, 9211, 9218, 9225, 9232, 9239,      9246, 9253, 9260, 9267, 9274, 9281, 9288, 9295, 9302, 9309, 9316,      9323, 9330, 9337, 9344, 9351, 9358, 9365, 9372, 9379, 9386, 9393,      9400, 9407, 9414, 9421, 9428, 9435, 9442, 9449, 9456, 9463, 9470,      9477, 9484, 9491, 9498, 9505, 9512, 9519, 9526, 9533, 9540, 9547,      9554, 9561, 9568, 9575, 9582, 9589, 9596, 9603, 9610, 9617, 9624,      9631, 9638, 9645, 9652, 9659, 9666, 9673, 9680, 9687, 9694, 9701,      9708, 9715, 9722, 9729, 9736, 9743, 9750, 9757, 9764, 9771, 9778,      9785, 9792, 9799, 9806, 9813, 9820, 9827, 9834, 9841, 9848, 9855,      9862, 9869, 9876, 9883, 9890, 9897, 9904, 9911, 9918, 9925, 9932,      9939, 9946, 9953, 9960, 9967, 9974, 9981, 9988, 9995, 10002,      10009, 10016, 10023, 10030, 10037, 10044, 10051, 10058, 10065,      10072, 10079, 10086, 10093, 10100, 10107, 10114, 10121, 10128,      10135, 10142, 10149, 10156, 10163, 10170, 10177, 10184, 10191,      10198, 10205, 10212, 10219, 10226, 10233, 10240, 10247, 10254,      10261, 10268, 10275, 10282, 10289, 10296, 10303, 10310, 10317,      10324, 10331, 10338, 10345, 10352, 10359, 10366, 10373, 10380,      10387, 10394, 10401, 10408, 10415, 10422, 10429, 10436, 10443,      10450, 10457, 10464, 10471, 10478, 10485, 10492, 10499, 10506,      10513, 10520, 10527, 10534, 10541, 10548, 10555, 10562, 10569,      10576, 10583, 10590, 10597, 10604, 10611, 10618, 10625, 10632,      10639, 10646, 10653, 10660, 10667, 10674, 10681, 10688, 10695,      10702, 10709, 10716, 10723, 10730, 10737, 10744, 10751, 10758,      10765, 10772, 10779, 10786, 10793, 10800, 10807, 10814, 10821,      10828, 10835, 10842, 10849, 10856, 10863, 10870, 10877, 10884,      10891, 10898, 10905, 10912, 10919, 10926, 10933, 10940, 10947,      10954, 10961, 10968, 10975, 10982, 10989, 10996, 11003, 11010,      11017, 11024, 11031, 11038, 11045, 11052, 11059, 11066, 11073,      11080, 11087, 11094, 11101, 11108, 11115, 11122, 11129, 11136,      11143, 11150, 11157, 11164, 11171, 11178, 11185, 11192, 11199,      11206, 11213, 11220, 11227, 11234, 11241, 11248, 11255, 11262,      11269, 11276, 11283, 11290, 11297, 11304, 11311, 11318, 11325,      11332, 11339, 11346, 11353, 11360, 11367, 11374, 11381, 11388,      11395, 11402, 11409, 11416, 11423, 11430, 11437, 11444, 11451,      11458, 11465, 11472, 11479, 11486, 11493, 11500, 11507, 11514,      11521, 11528, 11535, 11542, 11549, 11556, 11563, 11570, 11577,      11584, 11591, 11598, 11605, 11612, 11619, 11626, 11633, 11640,      11647, 11654, 11661, 11668, 11675, 11682, 11689, 11696, 11703,      11710, 11717, 11724, 11731, 11738, 11745, 11752, 11759, 11766,      11773, 11780, 11787, 11794, 11801, 11808, 11815, 11822, 11829,      11836, 11843, 11850, 11857, 11864, 11871, 11878, 11885, 11892,      11899, 11906, 11913, 11920, 11927, 11934, 11941, 11948, 11955,      11962, 11969, 11976, 11983, 11990, 11997, 12004, 12011, 12018,      12025, 12032, 12039, 12046, 12053, 12060, 12067, 12074, 12081,      12088, 12095, 12102, 12109, 12116, 12123, 12130, 12137, 12144,      12151, 12158, 12165, 12172, 12179, 12186, 12193, 12200, 12207,      12214, 12221, 12228, 12235, 12242, 12249, 12256, 12263, 12270,      12277, 12284, 12291, 12298, 12305, 12312, 12319, 12326, 12333,      12340, 12347, 12354, 12361, 12368, 12375, 12382, 12389, 12396,      12403, 12410, 12417, 12424, 12431, 12438, 12445, 12452, 12459,      12466, 12473, 12480, 12487, 12494, 12501, 12508, 12515, 12522,      12529, 12536, 12543, 12550, 12557, 12564, 12571, 12578, 12585,      12592, 12599, 12606, 12613, 12620, 12627, 12634, 12641, 12648,      12655, 12662, 12669, 12676, 12683, 12690, 12697, 12704, 12711,      12718, 12725, 12732, 12739, 12746, 12753, 12760, 12767, 12774,      12781, 12788, 12795, 12802, 12809, 12816, 12823, 12830, 12837,      12844, 12851, 12858, 12865, 12872, 12879, 12886, 12893, 12900,      12907, 12914, 12921, 12928, 12935, 12942, 12949, 12956, 12963,      12970, 12977, 12984, 12991, 12998, 13005, 13012, 13019, 13026,      13033, 13040, 13047, 13054, 13061, 13068, 13075, 13082, 13089,      13096, 13103, 13110, 13117, 13124, 13131, 13138, 13145, 13152,      13159, 13166, 13173, 13180, 13187, 13194, 13201, 13208, 13215,      13222, 13229, 13236, 13243, 13250, 13257, 13264, 13271, 13278,      13285, 13292, 13299, 13306, 13313, 13320, 13327, 13334, 13341,      13348, 13355, 13362, 13369, 13376, 13383, 13390, 13397, 13404,      13411, 13418, 13425, 13432, 13439, 13446, 13453, 13460, 13467,      13474, 13481, 13488, 13495, 13502, 13509, 13516, 13523, 13530,      13537, 13544, 13551, 13558, 13565, 13572, 13579, 13586, 13593,      13600, 13607, 13614, 13621, 13628, 13635, 13642, 13649, 13656,      13663, 13670, 13677, 13684, 13691, 13698, 13705, 13712, 13719,      13726, 13733, 13740, 13747, 13754, 13761, 13768, 13775, 13782,      13789, 13796, 13803, 13810, 13817, 13824, 13831, 13838, 13845,      13852, 13859, 13866, 13873, 13880, 13887, 13894, 13901, 13908,      13915, 13922, 13929, 13936, 13943, 13950, 13957, 13964, 13971,      13978, 13985, 13992, 13999, 14006, 14013, 14020, 14027, 14034,      14041, 14048, 14055, 14062, 14069, 14076, 14083, 14090, 14097,      14104, 14111, 14118, 14125, 14132, 14139, 14146, 14153, 14160,      14167, 14174, 14181, 14188, 14195, 14202, 14209, 14216, 14223,      14230, 14237, 14244, 14251, 14258, 14265, 14272, 14279, 14286,      14293, 14300, 14307, 14314, 14321, 14328, 14335, 14342, 14349,      14356, 14363, 14370, 14377, 14384, 14391, 14398, 14405, 14412,      14419, 14426, 14433, 14440, 14447, 14454, 14461, 14468, 14475,      14482, 14489, 14496, 14503, 14510, 14517, 14524, 14531, 14538,      14545, 14552, 14559, 14566, 14573, 14580, 14587, 14594, 14601,      14608, 14615, 14622, 14629, 14636, 14643, 14650, 14657, 14664,      14671, 14678, 14685, 14692, 14699, 14706, 14713, 14720, 14727,      14734, 14741, 14748, 14755, 14762, 14769, 14776, 14783, 14790,      14797, 14804, 14811, 14818, 14825, 14832, 14839, 14846, 14853,      14860, 14867, 14874, 14881, 14888, 14895, 14902, 14909, 14916,      14923, 14930, 14937, 14944, 14951, 14958, 14965, 14972), class = "Date"),          cases = c(168L, 199L, 214L, 230L, 267L, 373L, 387L, 443L,          579L, 821L, 1229L, 1014L, 831L, 648L, 257L, 203L, 137L, 78L,          82L, 69L, 45L, 51L, 45L, 63L, 55L, 54L, 52L, 27L, 24L, 12L,          10L, 22L, 42L, 32L, 52L, 82L, 95L, 91L, 104L, 143L, 114L,          100L, 83L, 113L, 145L, 175L, 222L, 258L, 384L, 755L, 976L,          879L, 846L, 1004L, 801L, 799L, 680L, 530L, 410L, 302L, 288L,          234L, 269L, 245L, 240L, 176L, 188L, 128L, 96L, 59L, 63L,          44L, 52L, 39L, 50L, 36L, 40L, 48L, 32L, 39L, 28L, 29L, 16L,          20L, 25L, 25L, 48L, 57L, 76L, 117L, 107L, 91L, 90L, 83L,          76L, 86L, 104L, 101L, 116L, 120L, 185L, 290L, 537L, 485L,          561L, 1142L, 1213L, 1235L, 1085L, 1052L, 987L, 918L, 746L,          620L, 396L, 280L, 214L, 148L, 148L, 94L, 107L, 69L, 55L,          69L, 47L, 43L, 49L, 30L, 42L, 51L, 41L, 39L, 40L, 38L, 22L,          37L, 26L, 40L, 56L, 54L, 74L, 99L, 114L, 114L, 120L, 114L,          123L, 131L, 170L, 147L, 163L, 163L, 160L, 158L, 163L, 124L,          115L, 176L, 171L, 214L, 320L, 507L, 902L, 1190L, 1272L, 1282L,          1146L, 896L, 597L, 434L, 216L, 141L, 101L, 86L, 65L, 55L,          35L, 49L, 29L, 55L, 53L, 57L, 34L, 43L, 42L, 13L, 17L, 20L,          27L, 36L, 47L, 64L, 77L, 82L, 82L, 95L, 107L, 96L, 106L,          93L, 114L, 102L, 116L, 128L, 123L, 212L, 203L, 165L, 267L,          550L, 761L, 998L, 1308L, 1613L, 1704L, 1669L, 1296L, 975L,          600L, 337L, 259L, 145L, 91L, 70L, 79L, 63L, 58L, 51L, 53L,          39L, 49L, 33L, 47L, 56L, 32L, 43L, 47L, 19L, 32L, 18L, 34L,          39L, 63L, 57L, 55L, 69L, 76L, 103L, 99L, 108L, 131L, 113L,          106L, 122L, 138L, 136L, 175L, 207L, 324L, 499L, 985L, 1674L,          1753L, 1419L, 1105L, 821L, 466L, 274L, 180L, 143L, 82L, 101L,          72L, 55L, 71L, 50L, 33L, 26L, 25L, 27L, 21L, 24L, 24L, 20L,          18L, 18L, 25L, 23L, 13L, 10L, 16L, 9L, 12L, 16L, 25L, 31L,          36L, 40L, 36L, 47L, 32L, 46L, 75L, 63L, 49L, 90L, 83L, 101L,          78L, 79L, 98L, 131L, 83L, 122L, 179L, 334L, 544L, 656L, 718L,          570L, 323L, 220L, 194L, 125L, 95L, 77L, 46L, 42L, 29L, 35L,          21L, 29L, 16L, 14L, 19L, 15L, 19L, 18L, 21L, 10L, 14L, 7L,          7L, 5L, 9L, 14L, 11L, 18L, 22L, 39L, 36L, 46L, 44L, 37L,          30L, 39L, 37L, 45L, 71L, 59L, 57L, 80L, 68L, 88L, 72L, 74L,          208L, 357L, 621L, 839L, 964L, 835L, 735L, 651L, 400L, 292L,          198L, 85L, 64L, 41L, 40L, 23L, 18L, 14L, 22L, 9L, 19L, 8L,          14L, 12L, 15L, 14L, 4L, 6L, 7L, 7L, 8L, 13L, 10L, 19L, 17L,          20L, 22L, 40L, 37L, 45L, 34L, 26L, 35L, 67L, 49L, 77L, 82L,          80L, 104L, 88L, 49L, 73L, 113L, 142L, 152L, 206L, 293L, 513L,          657L, 919L, 930L, 793L, 603L, 323L, 202L, 112L, 55L, 31L,          27L, 15L, 15L, 6L, 13L, 21L, 10L, 11L, 9L, 8L, 11L, 7L, 5L,          1L, 4L, 7L, 2L, 6L, 12L, 14L, 21L, 29L, 32L, 26L, 22L, 44L,          39L, 47L, 44L, 93L, 145L, 289L, 456L, 685L, 548L, 687L, 773L,          575L, 355L, 248L, 179L, 129L, 122L, 103L, 72L, 72L, 36L,          26L, 31L, 12L, 14L, 14L, 14L, 7L, 8L, 2L, 7L, 8L, 9L, 26L,          10L, 13L, 13L, 5L, 5L, 3L, 6L, 1L, 10L, 6L, 7L, 17L, 12L,          21L, 32L, 29L, 18L, 22L, 24L, 38L, 52L, 53L, 73L, 49L, 52L,          70L, 77L, 95L, 135L, 163L, 303L, 473L, 823L, 1126L, 1052L,          794L, 459L, 314L, 252L, 111L, 55L, 35L, 14L, 30L, 21L, 16L,          9L, 11L, 6L, 6L, 8L, 9L, 9L, 10L, 15L, 15L, 11L, 6L, 3L,          8L, 4L, 7L, 7L, 13L, 10L, 23L, 24L, 36L, 25L, 34L, 37L, 46L,          39L, 37L, 55L, 65L, 54L, 60L, 82L, 55L, 53L, 61L, 52L, 75L,          92L, 121L, 170L, 199L, 231L, 259L, 331L, 357L, 262L, 154L,          77L, 34L, 41L, 21L, 17L, 16L, 7L, 15L, 11L, 7L, 5L, 6L, 13L,          7L, 6L, 8L, 7L, 1L, 11L, 9L, 3L, 9L, 9L, 8L, 15L, 19L, 16L,          10L, 12L, 26L, 35L, 35L, 41L, 34L, 30L, 36L, 43L, 23L, 55L,          107L, 141L, 217L, 381L, 736L, 782L, 663L, 398L, 182L, 137L,          79L, 28L, 26L, 16L, 14L, 8L, 4L, 4L, 6L, 6L, 11L, 4L, 5L,          7L, 7L, 6L, 8L, 2L, 3L, 3L, 1L, 1L, 3L, 3L, 2L, 8L, 8L, 11L,          10L, 11L, 8L, 24L, 25L, 25L, 33L, 36L, 51L, 61L, 74L, 92L,          89L, 123L, 402L, 602L, 524L, 494L, 406L, 344L, 329L, 225L,          136L, 136L, 84L, 55L, 55L, 42L, 19L, 28L, 8L, 7L, 2L, 7L,          6L, 4L, 3L, 5L, 3L, 3L, 0L, 1L, 2L, 3L, 2L, 1L, 2L, 2L, 9L,          4L, 9L, 10L, 18L, 15L, 13L, 12L, 10L, 19L, 15L, 22L, 23L,          34L, 43L, 53L, 47L, 57L, 328L, 552L, 787L, 736L, 578L, 374L,          228L, 161L, 121L, 96L, 58L, 50L, 37L, 14L, 9L, 6L, 15L, 12L,          9L, 1L, 6L, 4L, 7L, 7L, 3L, 6L, 9L, 15L, 22L, 28L, 34L, 62L,          54L, 75L, 65L, 58L, 57L, 60L, 37L, 47L, 60L, 89L, 90L, 193L,          364L, 553L, 543L, 676L, 550L, 403L, 252L, 140L, 125L, 99L,          63L, 63L, 76L, 85L, 68L, 67L, 38L, 25L, 24L, 11L, 9L, 9L,          4L, 8L, 4L, 6L, 5L, 2L, 6L, 4L, 4L, 1L, 5L, 4L, 1L, 2L, 2L,          2L, 2L, 3L, 4L, 4L, 7L, 5L, 2L, 10L, 11L, 17L, 11L, 16L,          15L, 11L, 12L, 21L, 20L, 25L, 46L, 51L, 90L, 123L)), .Names = c("date",      "cases"), row.names = c(NA, -835L), class = "data.frame")


title: Jack-knife with time series models

text: **Introduction** I am aiming to forecast the annual growth rates for a number of macroeconomic indicators (denote one by $Y_t$). One of the tasks is to test the forecasting performance of rival time series models with and without exogenous variables ($X_t$, a $T\times k$ matrix). The list of rival models include:   1. AR(I)MA model (annual growth rates are unlikely to have "unit Roo", though the latter is either assumed or tested) $$A(L)Y_t =\mu+ B(L)\varepsilon_t$$   2. linear regression model with ARMA errors $$Y_t = X_t\beta + \eta_t, \ \ A(L)\eta_t = B(L)\varepsilon_t $$   3. lagged dependent variable model (autoregressive model with exogenous variables) $$A(L)Y_t = X_t\beta + \varepsilon_t$$    4. linear regression model $$Y_t = X_t\beta + \varepsilon_t$$ Where $\varepsilon_t$ is assumed to be a strong white noise, zero-mean constant variance i.i.d. process; $A(L)$ and $B(L)$ are autoregressive (of order $p$) and moving average (of order $q$) polynomials with $L$ - a back- shift (lag) operator. Note that the primary and the only goal is forecasting performance, thus any "good" properties of parameter estimates are of the secondary concern. All I need is to test for the most parsimonious, robust to starting conditions forecaster. Decision will be made with one of the `accuracy()` options, but first I need to obtain the material for the comparison. Models 1. and 2. are estimated by `auto.arima()` with default `"CSS-ML"` estimation method. Models 3. and 4. are estimated by ordinary least squares (`lm()`). $T$ is about $40$ quarters. **Approaches tried so far** To make the jack-knifed residuals the first approach denoted by "rolling" has been implemented. Starting from feasibly large sub-sample of time series data, parameters are estimated and an $h$ ahead forecast is done by the `predict()` function (EDIT: it is the same suggestion as in the first part Rob's answer to the second question). After that one point is added and the estimation\prediction steps are repeated. A weak point of such experiments is that the number of time ticks (sample size) used to estimate the parameters is different. While I would like to test the robustness to the starting conditions, keeping the sample size for estimation fixed. Bearing this in mind I tried to set the several subsequent values (EDIT: for the interval $k+p+q<t_0<t_1<T-h+1$) in $Y_t$ being missing values (NA). In models 2.-4. this also implies dropping the corresponding subsequent rows in data matrix $X_t$. Prediction for 3. and 4. is straightforward (the same `predict()` with omitted $X_t$ data rows works well). All my concerns are about models 1. and 2. With just the AR($p$) part the predictions are done sequentially $Y_{t+1|t} = \hat A(L)Y_t$. But with the presence of MA($q$) one could not (?) use the estimated parameters directly. From the Brockwell and Davis "Introduction to Time Series and Forecasting" Chapter 3.3 it follows that one needs an innovation algorithm to recursively estimate the $\theta_{n,j}$ from the specific system of equations that involves estimated autoregressive and moving average parameters. EDIT: these $\theta_{n,j}$ parameters are used to make the ARMA prediction, not the originally estimated parameters $\theta_{j}$. However it is remarked in the same chapter that $\theta_{n,j}$ asymptotically approaches $\theta_{j}$ if the process is invertible. It is not evident that 30-40 points is enough for the asymptotic result to be used even if is invertible. Notes: I don't want to restrict $q$ to zero, since I am not doing so in true out-of-sample forecasting. EDIT: also not that it is not missing value imputation problem, but forecasting experiment, that the trajectory is not supposed to bridge two sub-samples by this imputing the missing values. **Questions**   1. Does `auto.arima()` performs correctly with the presence of missing values inside of the sample? [Already answered by Rob.]   2. (The actually crucial part of this post) How to correctly forecast (NOT impute) these missed points from the ARMA model when both $p>0$ and $q>0$? (I hope there are the ways already implemented in R language, but I simply is missing something.) EDIT: since the parameters for ARMA parts are estimated correctly could I legally rearrange the arima object to include the estimated parameters and the data only for the first subsample and then use a predict function? EDIT2: I have tried to modify the estimated `mod` structure - the resulting forecast from `predict.Arima` is identical (double precision differences) to the forecast where I use the estimated MA and AR coefficients predicting $Y_{t+1|t}$ directly as $\hat A(L)(Y_t-X_t\hat \beta)+ X_t\hat \beta+\hat B(L)\hat \varepsilon_t $, without `KalmanForecast()`. This was expected since state space representation is supplied with the same estimated $\theta_j$, not $\theta_{n,j}$. So the only question remains is the difference between $\theta_j$ and $\theta_{n,j}$ significant to influence the point forecasts? I hope the answer is negative.


title: Standardized generalized hyperbolic distribution

text: I am interested in the standardized version (mean zero, variance one) of the generalized hyperbolic and the hyperbolic distribution. I want to include this in my analysis and therefore I need the dervations and the implementations. In R, there is an implementation of the standardized generalized hyperbolic distribution.   1. I need the used formula and how to derive it for this implementation, where can I get it?   2. Where can I see, what the code is doing? So where can I get a 'look behind the scenes' to see what is implemented?   3. I searched via google but I could not find any derviations of how to get a standardized version of the ghyp or hyp. Where can I find this?   4. I could use the R command to also fit a standardized hyperbolic distribution if I do the following: Inser for lambda=1 (this gives hyperbolic) and then use the `dsgh` command along together with the `optim` command to get a good fit of the standardized hyperbolic distribution to my data, is this correct? edit: OK, the answer to 2) is, that I just can load the package and enter the command without any specifications. I then get the code behind it. edit: I also noticed the sghFIT command, where can I find the theory behind this? Because they have to use a standardized version, so they need the formula and the derivation, this is what I need. edit: I just saw, that they also only use the dsgh command for fitting. Ok, my last edit, all comes back to the implementation of the dsgh command, the code for this command is:               dsgh <-       function(x, zeta = 1, rho = 0, lambda = 1, log = FALSE)      {         # A function implemented by Diethelm Wuertz              # Description:         #   Returns density of the sgh distribution              # FUNCTION:                  # Parameters:         if (length(zeta) == 3) {            lambda = zeta[3]            rho = zeta[2]            zeta = zeta[1]         }               # Compute Density:         param = .paramGH(zeta, rho, lambda)         ans = dgh(x, param[1], param[2], param[3], param[4], lambda, log)              # Return Value:         ans     }      He computes the density with               # Compute Density:             param = .paramGH(zeta, rho, lambda)             ans = dgh(x, param[1], param[2], param[3], param[4], lambda, log)      but I am not getting the idea behind the code. dgh is just giving the moments of the generalized hyperbolic distribution, see here? So where is the density? How does he do the modification to get a standardized version? Where is the theory behind this? So there must be a puplished paper or so, where it is described how to get the standardized version? I cannot find this?


title: KPSS test - output interpretation in stata

text: I did KPSS test for some variables in stata to check for stationarity; I want to interpret the the stata outputs, but I don't know how to do that. For instance, in the following case:                   KPSS test for mIlliq1         Maxlag = 13 chosen by Schwert criterion         Autocovariances weighted by Bartlett kernel         Critical values for H0: mIlliq1 is trend stationary         10%: 0.119  5% : 0.146  2.5%: 0.176  1% : 0.216              Lag order    Test statistic         0           .557         1           .309         2           .229         3           .188         4           .162         5           .144         6           .129         7           .118         8           .109         9           .102        10          .0974        11          .0934        12          .0901        13          .0873        Maxlag = 13 chosen by Schwert criterion        Autocovariances weighted by Bartlett kernel      How Can I interpret this result? Does it mean that the TS is stationary or not? If it isn't, how can I can make it a stationary time series? Thank to everyone.


title: A course in experimental design for data miners

text: I'm a computer scientist working in data mining. It's no secret to say that computer scientists are fairly poor at doing systematic experimental design and evaluation - the use of p-values and confidence estimates is considered advanced :). What I'd like to know if there are good courses/material to teach computer scientists about good experimental design. To make this more specific, I'll add the following information:   * The course should be targeted at graduate students who can be assumed to have a reasonable understanding of probability, but limited background in statistics.    * The course should focus on experimental design in "uncontrolled unnatural settings": in other words there is neither an underlying physical ground truth or a way to control the data gathering process (as with human subjects). Of course a good course will focus on fundamentals, but it should deal with this scenario in a significant way.    * A computational element would be a bonus but is not mandatory. We deal with lots of data, but can figure out computational issues ourselves if need be.


title: Efficient way to run multiple ANOVAs on subsets of data

text: I have a data set that is repeated-measures ratings of four devices (two two- level variables) for 6 dimensions (ease, confidence, comfort, control, size and fit), with a two-level between-subjects variable. Like so:               ptcip / grp / fdback / dur / assem / accep / freq / t.vs.r / attrib / meas / d.rating     1    RA   binary    short       <NA>  <NA> <NA> rating       ease    9         2     1    RA   binary    short       <NA>  <NA> <NA> rating confidence    7         1     1    RA   binary    short       <NA>  <NA> <NA> rating    comfort    6         4     1    RA   binary    short       <NA>  <NA> <NA> rating    control    5         3     1    RA   binary    short       <NA>  <NA> <NA> rating       size    7        -1     1    RA   binary    short       <NA>  <NA> <NA> rating        fit    6         0     1    RA   binary      med       <NA>  <NA> <NA> rating       ease    9         6     1    RA   binary      med       <NA>  <NA> <NA> rating confidence    5         6     1    RA   binary      med       <NA>  <NA> <NA> rating    comfort    6         5     1    RA   binary      med       <NA>  <NA> <NA> rating    control    9        -1     1    RA   binary      med       <NA>  <NA> <NA> rating       size    2        -1     1    RA   binary      med       <NA>  <NA> <NA> rating        fit    8        -2      I need to run separate ANOVAs on each of the dimensions (we know from other analyses that the dimensions are rated differently, so including it as a factor wouldn't be very informative). What's the most efficient way to write this--split each of the dimensions off into separate data frames, or specify aov to be run on data iff attrib == "ease"? Seems like the latter would be better, but I'm new to R. The model I am currently specifying is               aov(d.ratings ~ grp*fdback*dur + Error(particip/(fdback*dur))+grp, PDdonly)      but as that collapses over the six dimensions, it's not actually meaningful atm.


title: Formula to calculate a t-distribution

text: I am writing an application that will be dealing with <30 observations in a normal distribution. My understanding is that this point I would need to use t-distribution. The thing is, this is easy enough to look up in a table. However, I need this to be something to calculate it programmatically. I can't seem to find a formula to find the t-distribution score. How can I build this table so to speak? I am writing this in PHP. I should also note that I am by no means a math/stat expert. So if my use of terms/vocabulary seems off just bear with me.


title: Selecting priors for logistic functions

text: I have this confusion related to how to select priors for a logistic regression By Bayes theorem $P(\theta|D) = \frac{P(D|\theta) * P(\theta)}{P(D)}$. Now my likelihood $P(D|\theta)$ is given by logistic function $P(y=1) = \frac{1}{1+exp^{(-\theta'x)}}$. I was thinking what sort of priors I can use for P(theta). Can I use gamma distribution for P(theta). Will that help me to calculate $\int P(\theta)*P(D|\theta) d\theta$


title: Poisson distribution conditional on Poisson-distributed mean

text: I am looking for some details on the following distribution: assume the probability of a random variable $n \sim \text{Pois}(k)$ is conditional on $k \sim \text{Pois}(\mu)$. After some trying I figured out that the conditional probability $\Pr(n\mid\mu)$ should be $$ \Pr(n\mid\mu) = \frac{e^{-\mu/e\,(e-1)}}{n!} \, T_n\left(\mu/e\right) \, \text{,} $$ where the $T_n(x)$ are the Touchard polynomials. Does anyone know the name of this distribution? I've worked out what I need here, but I would much rather just quote this from somewhere. Any references? (I posted the same question to http://math.stackexchange.com/ without reply, then I discovered this here might be a better place to ask. If a cross-post is frowned upon, I will remove one of the question.)


title: Priors for parameters of normal distribution leading to same results as frequentist formula

text: Given a sample vector $x$ of size $N$ from a normally distributed population. With frequentist methods the population mean is estimated as $\hat{\mu}=\frac{\Sigma{}x_i}{N}$, population sigma is estimated as $\hat{\sigma}=\sqrt{\frac{\Sigma{(x_i - \hat{\mu})^2}}{N - 1}}$. Is there a Bayesian prior for the population mean and population sigma that will lead to the same estimates as the above frequentist formulas for any $x$?


title: Cross-validating time-series analysis

text: I've been using the caret package in R to build predictive models for classification and regression. Caret provides a unified interface to tune model hyper-parameters by cross validation or boot strapping. For example, if you are building a simple 'nearest neighbors' model for classification, how many neighbors should you use? 2? 10? 100? Caret helps you answer this question by re-sampling your data, trying different parameters, and then aggregating the results to decide which yield the best predictive accuracy. I like this approach because it is provides a robust methodology for choosing model hyper-parameters, and once you've chosen the final hyper-parameters it provides a cross-validated estimate of how 'good' the model is, using accuracy for classification models and RMSE for regression models. I now have some time-series data that I want to build a regression model for, probably using a random forest. What is a good technique to assess the predictive accuracy of my model, given the nature of the data? If random forests don't really apply to time series data, what's the best way to build an accurate ensemble model for time series analysis?


title: Is a power of an inverse gamma random variable itself inverse gamma?

text: If $X$ is an inverse gamma distributed random variable, then would $X^p$ also be distributed as inverse gamma? I found someone asked about the square root of inverse gamma, but I didn't find a direct answer to that question.


title: Mean absolute percentage error (MAPE) in Scikit-learn

text: How can we calculate the Mean absolute percentage error (MAPE) of our predictions using Python and scikit-learn? From the docs, we have only these 4 metric functions for Regressions:   * metrics.explained_variance_score(y_true, y_pred)   * metrics.mean_absolute_error(y_true, y_pred)    * metrics.mean_squared_error(y_true, y_pred)   * metrics.r2_score(y_true, y_pred)


title: Using incomplete regression output to determine the value of the test statistic

text: I have incomplete regression output subsisting of the sample size, the standard error of the regression coefficient, and the p-value for two variables in a regression model. Can I determine what the regression coefficients are? Can I infer which coefficient indicates a larger value for the regression coefficient?


title: Results on continuity corrections

text: For me, the first thing that comes to mind when I come across the term _continuity correction_ is that when $X\sim\mathrm{Bin}(n,p)$, one approximates $\Pr(X\le x)=\Pr(X<x)$ by $\Pr(Y<x+1/2)$ where $Y$ is normally distributed and has the same expectation and variance as $X$. Of course one may also use continuity corrections with Poisson distributions, and among the automatic list of "Questions that may already have your answer", I find continuity corrections for Pearson and McNemar chi-square tests and Yates' continuity correction for $2\times2$ tables. The fact that one gets reasonably close approximations even for fairly small $n$ is easily verified. The fact that things work well as $n\to\infty$ is easy to see. But what about precisely stated and proved mathematical results? Those might include precise error bounds, things about rates of convergence, and other things. What such results are known and where are they found?


title: ARIMA model sensitivity

text: Just wondering if a model was sensitive to the overall assumptions made when building the ARIMA model, what would be the impact on the model if any of these assumptions were violated?


title: Compare coefficients across groups in panel data set

text: I want to test for differences in (regression) coefficients across two groups. I understand, that this is possible by including interaction terms and performing a Chow-Test. But what do you do if you can't identify the groups in your panel dataset anymore? Have a look at this example. This may clarify my problem: Imagine 2 groups of people. Over 1 month, you record how much everybody spend on ice-cream (4 sorts), how they rated the ice-cream and what size of ice- cream they got. Your dependent variable is rating, your independent variables are cost and size.                  +-------------------------------------------------+         | day | ice-cream | rating | cost  | size  | group|         |-----+-----------+--------+-------+--------------+      1. |  1  |  vanilla  |  2     |  3    |  4    | 1    |      2. |  1  |  choco    |  1     |  4    |  3    | 0    |      3. |  1  |  vanilla  |  2     |  2    |  1    | 1    |      4. |  1  |  hazelnut |  9     |  2    |  1    | 0    |      5. |  2  |  hazelnut |  2     |  1    |  1    | 0    |      6. |  2  |  hazelnut |  5     |  1    |  1    | 0    |      7. |  2  |  vanilla  |  6     |  5    |  2    | 1    |      8. |  3  |  berry    |  1     |  6    |  3    | 0    |      9. |  3  |  berry    |  6     |  4    |  2    | 0    |     10. |  3  |  berry    |  7     |  3    |  4    | 1    |     11. |  3  |  berry    |  5     |  2    |  1    | 1    |     12. |  3  |  vanilla  |  4     |  1    |  2    | 0    |     13. |  3  |  vanilla  |  4     |  5    |  2    | 0    |     etc |  .  |  ......   |  .     |  .    |  ...  | .    |            +-------------------------------------------------+      **Not caring** about groups at all, I could just aggregate the data to a day- icecream level (because I want to consider fixed effects at the ice-cream level). My panel dataset would look like this:                   +------------------------------------------+         | day | ice-cream | rating | cost  | size  |          |-----+-----------+--------+-------+-------+      1. |  1  |  vanilla  |  2     |  2.5  |  2.5  |       2. |  1  |  choco    |  1     |  4    |  3    |       3. |  1  |  hazelnut |  9     |  2    |  1    |       4. |  2  |  hazelnut |  3.5   |  1    |  1    |       5. |  2  |  vanilla  |  6     |  5    |  2    |       6. |  3  |  berry    |  4.75  |  3.7  |  2.5  |       7. |  3  |  vanilla  |  4     |  3    |  2    |      etc |  .  |  ......   |  .     |  .    |  ...  |          +------------------------------------------+      However, this doesn't allow me to calculate a regression model for each group. So, I followed a different path: I split the raw data set according to groups (one set for each group) and aggregated to the daily level afterwards. The two panel datasets looks like this: **Group 1**                   +-------------------------------------------------+         | day | ice-cream | rating | cost  | size  | group|         |-----+-----------+--------+-------+--------------+      1. |  1  |  vanilla  |  1.5   |  2.5  |  2.5  | 1    |      2. |  2  |  vanilla  |  6     |  5    |  2    | 1    |      3. |  3  |  berry    |  7     |  3    |  4    | 1    |      4. |  3  |  vanilla  |  6     |  2.5  |  2.5  | 1    |     etc |  .  |  ......   |  .     |  .    |  ...  | .    |            +-------------------------------------------------+      **Group 2**                   +-------------------------------------------------+         | day | ice-cream | rating | cost  | size  | group|         |-----+-----------+--------+-------+--------------+      1. |  1  |  choco    |  1     |  4    |  3    | 0    |      2. |  1  |  hazelnut |  9     |  2    |  1    | 0    |      3. |  2  |  hazelnut |  3.5   |  1    |  1    | 0    |      4. |  3  |  berry    |  3.5   |  5    |  2.5  | 0    |      5. |  3  |  vanilla  |  4     |  3    |  2    | 0    |     etc |  .  |  ......   |  .     |  .    |  ...  | .    |            +-------------------------------------------------+      For each dataset, I run a fixed-effects regression. This leaves me with two model outputs which I want to compare and see whether the coefficients are different. Any ideas how to test coefficients across models in this case? Thank you very much! Simon


title: Poisson regression to estimate relative risk for binary outcomes

text: **Brief Summary** Why is it more common for logistic regression (with odds ratios) to be used in cohort studies with binary outcomes, as opposed to poisson regression (with relative risks)? **Background** Undergraduate and graduate statistics and epidemiology courses, in my experience, generally teach that logistic regression should be used for modelling data with binary outcomes, with risk estimates reported as odds ratios. However, poisson regression (and related: quasipoisson, negative binomial, etc.) can also be used to model data with binary outcomes and, with appropriate methods (e.g. robust sandwich variance estimator), provides valid risk estimates and confidence levels. E.g.,   * Greenland S., Model-based estimation of relative risks and other epidemiologic measures in studies of common outcomes and in case-control studies, Am J Epidemiol. 2004 Aug 15;160(4):301-5.   * Zou G., A modified poisson regression approach to prospective studies with binary data, Am J Epidemiol. 2004 Apr 1;159(7):702-6.    * Zou G.Y. and Donner A., Extension of the modified Poisson regression model to prospective studies with correlated binary data, Stat Methods Med Res. 2011 Nov 8. From poisson regression, relative risks can be reported, which some have argued are easier to interpret compared to odds ratios, especially for frequent outcomes, and especially by individuals without a strong background in statistics, see Zhang J. and Yu K.F., What's the relative risk? A method of correcting the odds ratio in cohort studies of common outcomes, JAMA. 1998 Nov 18;280(19):1690-1. From reading the medical litterature, among cohort studies with binary outcomes it seems that it is still far more common to report odds ratios from logistic regressions rather than relative risks from poisson regressions. **Questions** For cohort studies with binary outcomes:   1. Is there good reason to report odds ratios from logistic regressions rather than relative risks from poisson regressions?   2. If not, can the infrequency of poisson regressions with relative risks in the medical literature be attributed mostly to a lag between methodological theory and practice among scientists, clinicians, statisticians, and epidemiologists?   3. Should intermediate statistics and epidemiology courses include more discussion of poisson regression for binary outcomes?   4. Should I be encouraging students and colleagues to consider poisson regression over logistic regression when appropriate?


title: What is the distribution of the ratio of two t-distributed random variables?

text: x is t-distributed; y is t-distributed. How is x/y distributed? Does it have a closed-form formula?


title: Wilcoxon one-tailed test

text: It's regarding a one-tailed Wilcoxon sign test. I am comparing (fictional) ratings of liking before and after a woman shaves off the hair. The design is complicated. I will do my best to explain. There are 9 women. Participants are asked to rate how much they like each women. After that they are asked if the women should shave her hair. Those who answered yes, are asked to rate how much they like the woman after she has shaved her hair while those who answer no move to rating the next woman. (please remember this is fictional, my actual research is diff. but I don't want to post it here so using a fictional research which is nearly similar). So basically, I have pre and post shaving ratings. Participant one may have 5 pre and post ratings while part. 2 has 7 pre and post ratings. * * * I expect pre shaving ratings to be more than post shaving ratings, hence one- tailed. This is the output from SPSS: Descriptive Statistics   N Mean Std. Deviation Minimum Maximum Percentiles   25th 50th (Median) 75th rating1 379 4.000 1.08588 1.00 5.00 3.0000 4.0000 5.0000 rating2 379 3.500 1.70245 1.00 5.00 3.0000 4.0000 5.0000 Ranks   N Mean Rank Sum of Ranks rating2 - rating1 Negative Ranks 135a 63.00 7975.00 Positive Ranks 0b .00 .00 Ties 254c   Total 389   a rating2 < rating1   b rating2 > rating1   c rating2 = rating1 Test Statisticsa                                          rating2 - rating1                Z              -10.331b         Asymp. Sig. (2-tailed)     .000      a Wilcoxon Signed Ranks Test   b Based on positive ranks. As far as I can see, the difference is signficant p<0.001 (two-tailed). How about one-tailed? How do I know if it is signficant at one-tailed as my hypothesis is one-tailed...I have predicted that post shaving ratings will be less than pre-shaving.