title: Equivalence tests for STP3 - location of maximum power

text: Lehmann and Romano show in their Th 3.7.1 (Lehmann and Romano 2005, p81) when a UMP equivalence test $\phi_1$ exists and that it has a power that is maximised at some point $\theta^\star$ inside the equivalence region $[\theta_1, \theta_2]$. Now suppose we have a second test $\phi_2$ that only differs in the upper point of the equivalence region $[\theta_1, \theta_3]$, with $\theta_3>\theta_2$. Then I think the power of the second test should be maximized at a point $\theta^\star_2$ that is larger than $\theta^\star$, and that this should follow from the stochastic ordering of the distribution family induced by the test statistic (assumed to be STP3 in Th 3.7.1). Does anyone have any hints or suggestions on how to prove this?


title: Linear regression in matrix notation

text: I was reading Wikipedia's article on linear regression and I realized that I don't quite get the same result for the normal equation: $$\beta = (X^{T}X)^{-1}X^{T}{\bf{y}} = (\frac{1}{n}\sum {\bf{x}}_{i}{\bf{x}}_{i}^{T})^{-1}(\frac{1}{n}\sum {\bf{x}}_{i}{{y}}_{i})$$ where: $${\bf{y}}=\begin{pmatrix}y_{1}\\\y_{2}\\\ ...\\\ y_{n}\end{pmatrix}$$ $${\bf{X}}=\begin{pmatrix}{\bf{x}}_{1}^{T}\\\\{\bf{x}}_{2}^{T}\\\ ...\\\ {\bf{x}}_{n}^{T}\end{pmatrix}$$ $${\bf{\beta}}=\begin{pmatrix}\beta_{1}\\\\\beta_{2}\\\ ...\\\ \beta_{p}\end{pmatrix}$$ Supose we want to minimize an error function $E_{D}$ on the parameters $\beta$: $$E_{D} = \frac{1}{2} \sum_{n=1}^{N} \\{y_{n}-\beta^{T}x_{n}\\}^{2}$$ so we use the gradient on $\beta$ in order to obtain: $$\nabla E_{D} = \sum_{n=1}^{N} \\{y_{n}-\beta^{T}x_{n}\\}x_{n}^{T}=0$$ and solving for $\beta$ $$\beta^{T}=\frac{\sum_{n=1}^{N}y_{n}x_{n}^{T}}{\sum_{n=1}^{N}x_{n}x_{n}^{T}}$$ Taking $\bf{X}$ as described before, we find the correct denominator because $${\bf{X}}^{T}{\bf{X}}=( \bf{x}_{1} \thinspace \bf{x}_{2} ... \bf{x}_{N} ) \begin{pmatrix}{\bf{x}}_{1}^{T}\\\\{\bf{x}}_{2}^{T}\\\ ...\\\ {\bf{x}}_{n}^{T}\end{pmatrix}=\sum_{n=1}^{N}x_{n}x_{n}^{T}$$ However, in the denominator my calculation produces a row vector instead of a column vector. What did I get wrong?


title: multivariate analysis for time/place where animals are sighted

text: I have many marked animals that are sighted occasionally. So I have a table with the following columns:               Animal  DateTime  Lat  Lon      Each animal (about 500) has been sighted many times, so this table has about 60k rows. I wanna apply a multivariate analysis, where each animal will become a point in the plot, and the proximity of two such points will give me an idea of the proximity of these two animals in space and time (i.e. they are moving somewhat together). To achieve this, each row of a new table will represent one animal. For the columns, each one should represent at least a day (or a part of a day). And the values in the table should represent the place, but: is there a way to represent lat/lon in a single value? Or should I put two columns for each time? Or anything else? Thanks in advance!


title: Determining beta distribution parameters $\alpha$ and $\beta$ from two arbitrary points

text: Suppose I have two points $(p_1,x_1)$ and $(p_2,x_2)$ where $p_i$ is a probability on the beta CDF and $x_i$ is a value on that same CDF. How would I go about determining the beta distribution shape parameters $\alpha$ and $\beta$ with only these two points? To clarify, I'm assuming I already know the lower and upper bounds of the distribution L and U.


title: Please help me interpret this goodness of fit test

text: In the game of Bridge, 52 cards are distrubuted randomly among 4 players, each player ending up with 13 cards. I implemented a random number generator using the RNGCryptoServiceProvider form the Microsoft .Net Framerowk to first draw a number from the range 1-52, then from the range 1-51 and so on up to the range 1-2. To see if at a given range of 1 to x, where 2>= x <=52, any number in this range is equally likely to be selected I implemented the following test:   * Draw 52* 10 million times from a range of 1-52 and keep track of how often a number from this range comes up. Apply a goodness of fit test against an expected frequency of 10 million for each number.      * Draw 51* 10 million times from a range of 1-51 and keep track of how often a number from this range comes up. Apply a goodness of fit test against an expected frequency of 10 million for each number. And so forth until we reach the range 1-2. Below you find the p values for each number of cards, when it is lower than 0.05 the observered number of draws for each number in the range is printed. In the sample below for two ranges the p value indeed was lower than 0.05. Should I now conclude that the RNGCryptoServiceProvider is not capable of producing an equal distrbution for all ranges in the range 1-52, 1-51, ... 1-2? Test Name: CardChances Test Outcome: Passed Result StandardOutput:   Debug Trace: P value for 52 cards=0.79474227840984 P value for 51 cards=0.497191145912234 P value for 50 cards=0.391892735550181 P value for 49 cards=0.988683559265785 P value for 48 cards=0.153114936344037 P value for 47 cards=0.927327644547415 P value for 46 cards=0.493708839113774 P value for 45 cards=0.35952495612322 P value for 44 cards=0.979430525949886 P value for 43 cards=0.192281134864909 P value for 42 cards=0.53964228753539 P value for 41 cards=0.175093762735034 P value for 40 cards=0.801225853918414 P value for 39 cards=0.156440626522907 P value for 38 cards=0.186908820122156 P value for 37 cards=0.0828295408060195 P value for 36 cards=0.81408556973856 P value for 35 cards=0.805491467886622 P value for 34 cards=0.336741413336137 P value for 33 cards=0.498471084831715 P value for 32 cards=0.824924933331495 P value for 31 cards=0.163033897868359 P value for 30 cards=0.0967583951496144 P value for 29 cards=0.369476297712753 P value for 28 cards=0.0243850151443583 Card number 1 was drawn 10000359 times. Card number 2 was drawn 9998315 times. Card number 3 was drawn 10004012 times. Card number 4 was drawn 9988701 times. Card number 5 was drawn 9997630 times. Card number 6 was drawn 9996076 times. Card number 7 was drawn 10003135 times. Card number 8 was drawn 9997429 times. Card number 9 was drawn 9999151 times. Card number 10 was drawn 10001832 times. Card number 11 was drawn 10001967 times. Card number 12 was drawn 9994442 times. Card number 13 was drawn 9999789 times. Card number 14 was drawn 9996959 times. Card number 15 was drawn 9998547 times. Card number 16 was drawn 10005698 times. Card number 17 was drawn 10006844 times. Card number 18 was drawn 10001617 times. Card number 19 was drawn 10001233 times. Card number 20 was drawn 9999269 times. Card number 21 was drawn 10005547 times. Card number 22 was drawn 9994631 times. Card number 23 was drawn 10006450 times. Card number 24 was drawn 10001226 times. Card number 25 was drawn 9997165 times. Card number 26 was drawn 10000925 times. Card number 27 was drawn 9999416 times. Card number 28 was drawn 10001635 times. P value for 27 cards=0.614170487404118 P value for 26 cards=0.880271821699275 P value for 25 cards=0.923499170502961 P value for 24 cards=0.624287793606294 P value for 23 cards=0.463782431271392 P value for 22 cards=0.0918275795728352 P value for 21 cards=0.848413464526898 P value for 20 cards=0.437303155192944 P value for 19 cards=0.605302075922096 P value for 18 cards=0.0293628427079152 Card number 1 was drawn 9999153 times. Card number 2 was drawn 10004923 times. Card number 3 was drawn 9993127 times. Card number 4 was drawn 10001637 times. Card number 5 was drawn 10000355 times. Card number 6 was drawn 9999963 times. Card number 7 was drawn 10002201 times. Card number 8 was drawn 9999522 times. Card number 9 was drawn 9997863 times. Card number 10 was drawn 10004625 times. Card number 11 was drawn 10001398 times. Card number 12 was drawn 10005456 times. Card number 13 was drawn 10001133 times. Card number 14 was drawn 9997193 times. Card number 15 was drawn 9992970 times. Card number 16 was drawn 10007465 times. Card number 17 was drawn 9996850 times. Card number 18 was drawn 9994166 times. P value for 17 cards=0.463159193420478 P value for 16 cards=0.901074254855155 P value for 15 cards=0.623196060523292 P value for 14 cards=0.16131153037034 P value for 13 cards=0.39322339086841 P value for 12 cards=0.239425750163672 P value for 11 cards=0.465057310698378 P value for 10 cards=0.147910602845352 P value for 9 cards=0.223291445001854 P value for 8 cards=0.567005883891805 P value for 7 cards=0.660822108546329 P value for 6 cards=0.465832956027658 P value for 5 cards=0.586291324248852 P value for 4 cards=0.537424479438887 P value for 3 cards=0.472376330829859 P value for 2 cards=0.767529320894893


title: How many measurements should be done?

text: I am measuring time of a computer operation. The operation should run roughly same time each time I measure it. How many times should I measure it to get good average and standard deviation? originally posted here: http://physics.stackexchange.com/questions/64917/how- many-measurements-should-be-done , but I think that this forum is more appropriate for my question


title: Query on Quants

text: I am working on decision trees for the first time at job. I have done lot of research on CHAID and CART algorithms but find different answers to a very simple question given below : **What kind of target variables CART can have?** I understand that CART can help both in prediction and classification. I confirm that the target variable for regression tree is continuous. I am getting different answers in various research papers with regards to the target variable of classification tree. Can someone please help me with this? Further, I need to do some analysis on the following : 1)Identifying frauds : I intend to use classification trees of CART/CHAID/logistic regression 2) Forecasting losses : I intend to use linear regression or regression trees of CART 3) Identifying cheque bounce customers : I intend to use classification trees OF CART/CHAID or logistic regression Kindly suggest if this is the right way to go.... thanks quants_mum


title: Can student's t and ANOVA be used to analyze a non scalar invariant interval variable?

text: On Cognitive Sciences I asked the question which statistical test can be used to analyze NASA TLX results. Here I will ask a similar, but more generalized question. Suppose an interval variable which might not be fully scalar invariant:   1. How to verify how scalar invariant results are?   2. Can a student's t test or ANOVA still be be used to analyze whether groups of data are significantly different, or this is entirely out of the question?   3. How/to which extent does scalar invariance influence the usefulness of these tests? If I am asking an altogether wrong question, pointing out why that is so would also constitute a useful answer. ;)


title: linear regression intercept does not match

text: I have done a linear regression in R, using glm function. The calculated intercept says 0.98, but when I plot it, it does not seem to hit the estimated intercept on Y axis. Its far below. Here are my data and function:               event = c(2.2, 6.4, 3.4, 10.2, 4.45, 2.65, 8.25, 4.65, 3, 6.5, 5.25,      8.65, 7.25, 6.4, 7.75, 7.45)          c(230208, 813178, 316617, 1531919, 576869, 270148, 1090947, 562643,      439885, 745741, 666454, 1078175, 924429, 784333, 1091289, 948062)          fit=glm(event~size)          Call:  glm(formula = chr.co.count.wt ~ size)          Coefficients:     (Intercept)         size         9.783e-01    6.528e-06            Degrees of Freedom: 15 Total (i.e. Null);  14 Residual     Null Deviance:      83.08      Residual Deviance: 2.849    AIC: 23.8                         plot(size,events,col="blue",pch=16,xlab="size",ylab="events",ylim=c(0,12),frame.plot=FALSE,xlim=c(0,2000000),axes = F)          axis(side = 1,at = c(0,0.5e6,1e6,1.5e6),labels =  c(0,0.5e6,1e6,1.5e6))     axis(side = 2,at = seq(from = 0,to = 12,by = 0.5),labels = seq(from = 0,to = 12,by = 0.5))     abline(fit.wt)      ![enter image description here](http://i.stack.imgur.com/iqDK2.png) Why is this discrepancy ? Am i missing something here ? I have also checked the std. err which is 0.27, still higher than what is being observed on plot. Thank you.


title: Independent Variable groups w/ n=1 possible for Kruskal Wallis?

text: I am working with clinical data with 7 groups that constitute my independent variable and a not normally distributed dependent variable. I know I have to use KW test to check for significant differences between my groups. The problem is that 3 of the 7 groups have an n=1. I know KW is not as sensitive to unequal sample sizes as parametric tests but is n=1 even possible? Would I be violating one of the assumptions (same overall dist. of data between groups) if I run KW? I cannot combine the values of the 3 groups as each group is very different from the other. I also don't know if its possible to remove these 3 groups for the KW test as I would be misrepresenting the reality of things Please let me know what you think. Thank you.


title: Bayesian AB testing

text: I am running an AB Test on a page that receives only 5k visits per month. It would take too long to reach traffic levels necessary to measure a +-1% difference between the test and control. I have heard that I can use Bayesian stats to give me a good chance of determining whether the test outperformed. How can I do use Bayesian stats to analyze my current data?                       Visitors    Conversions     Control 1345         1165     Test A  961          298     Test B  1274         438


title: Fitting ARIMA with a drift on R

text: I had some problems fitting an ARIMA - the data are               FTSE <-log(EuStockMarkets[,"FTSE"])       The following link explains the problem and gives the solution http://www.stat.pitt.edu/stoffer/tsa3/Rissues.htm#Issue2   The solution sounded effective, untill I tried some $ARIMA(p,d,q)$ other than $ARIMA(1,1,0)$: with $ARIMA(1,1,1)$, $ARIMA(2,1,1)$ and $ARIMA(1,1,2)$ the procedure (1+) and the procedure (2) give different coefficients, in particular with the procedure (1+) it says that the AR part is not stationary (indeed changing method from 'css' to 'ml' the coeffiecient of $AR$ is 1 or above).   These are the codes for $ARIMA(1,1,1)$:               dFTSE<-diff(FTSE)     model2<-Arima(dFTSE, order = c(1, 0, 1))     model2     model1pluscss<-Arima(FTSE,order=c(1,1,1), xreg=1:length(FTSE))     model1plusml<-Arima(FTSE,order=c(1,1,1), xreg=1:length(FTSE), method='ML')     model1plusml      On the other side the two procedures produce the same output also in $ARIMA(0,1,1)$ and $ARIMA(0,1,2)$. Is there a theoretical explanation that I can't see or something else?


title: SVAR, Cholesky decomposition and impulse-response function in R

text: I have two sets of data from the FRED database: real GDP (y) and GDP deflator (p) and I want to be able to use R in order to estimate a VAR(p) (p determined by AIC) process and generate the sets of impulse-response functions with the short-run assumptions (Sims, 1980) which utilizes the Cholesky decomposition. Since this is a website for learning, this is a detailed explanation of the process so that this post could actually "teach" something to some of you. If you can help me and you already know how to do it, the first paragraph is actually what I am looking for. Impulse-response analysis is the analysis of the dynamic response of an economic variable of interest (e.g. real GDP) to shocks in other economic variables such as demand shocks (e.g. inflation) or supply shocks ( e.g. technology). In order to do that, we may want to use a reduced form vector autoregressive process (RVAR): ![RVAR\(1\)](http://latex.codecogs.com/gif.latex?Y_%7Bt%7D%20%3D%20B_%7B1%7DY_%7Bt-1%7D%20+%20B_%7Bt-2%7D%20+%20...%20+%20B_%7Bt-p%7D%20+%20%5Cvarepsilon%20_%7Bt%7D) Where: ![](http://latex.codecogs.com/gif.latex?Y_%7Bt%7D%20%3D%20%5Cbegin%7Bbmatrix%7D%20Y_%7B1%2Ct%7D%5C%5C%20...%5C%5C%20Y_%7Bk%2Ct%7D%5C%5C%20%5Cend%7Bbmatrix%7D) ![](http://latex.codecogs.com/gif.latex?B%20%3D%20%5Cbegin%7Bbmatrix%7D%20b_%7B11%7D%20%26%20...%20%26%20b_%7B1k%7D%5C%5C%20...%20%26%20...%20%26%20...%20%5C%5C%20b_%7Bk1%7D%20%26%20...%20%26%20b_%7Bkk%7D%20%5Cend%7Bbmatrix%7D) ![](http://latex.codecogs.com/gif.latex?%5Cvarepsilon%20_%7Bt%7D%20%3D%20%5Cbegin%7Bbmatrix%7D%20%5Cvarepsilon%20_%7B1%2Ct%7D%5C%5C%20...%5C%5C%20%5Cvarepsilon%20_%7Bk%2Ct%7D%5C%5C%20%5Cend%7Bbmatrix%7D) The variance-covariance matrix of this process is as follows: ![Variance-covariance matrix](http://latex.codecogs.com/gif.latex?E%5Cvarepsilon%20_%7Bt%7D%5Cvarepsilon%20_%7Bt%7D%5E%7B%27%7D%20%3D%20%5Csum) And it is a symmetric, positive definite matrix whose off-diagonal values are non zero, which means that the error terms are mutually correlated. Consequently, our attempt to trace the dynamic responses of our variable of interest will be hindered. One solution to this problem is the use of a structural vector autoregressive process (SVAR): ![Variance-covariance matrix](http://latex.codecogs.com/gif.latex?A_%7B0%7DY_%7Bt%7D%20%3D%20A_%7B1%7DY_%7Bt-1%7D%20+%20A_%7B2%7DY_%7Bt-2%7D%20+%20...%20+%20A_%7Bp%7DY_%7Bt-p%7D%20+%20u_%7Bt%7D) Where A0 is the contemporaneous relations between the k variables. We can multiply both sides of this equation by the inverse of the contemporaneous effect: ![Structural form VAR](http://latex.codecogs.com/gif.latex?Y_%7Bt%7D%20%3D%20A_%7B0%7D%5E%7B-1%7DA_%7B1%7DY_%7Bt-1%7D%20+%20A_%7B0%7D%5E%7B-1%7DA_%7B2%7DY_%7Bt-2%7D%20+%20...%20+%20A_%7B0%7D%5E%7B-1%7DA_%7Bp%7DY_%7Bt-p%7D%20+%20A_%7B0%7D%5E%7B-1%7Du_%7Bt%7D) Where: ![](http://latex.codecogs.com/gif.latex?%5Cvarepsilon%20_%7Bt%7D%20%3D%20A_%7B0%7D%5E%7B-1%7Du_%7Bt%7D) ![](http://latex.codecogs.com/gif.latex?A_%7Bj%7D%20%3D%20A_%7B0%7DB_%7Bj%7D) ![](http://latex.codecogs.com/gif.latex?E%5Cvarepsilon%20_%7Bt%7D%5Cvarepsilon%20_%7Bt%7D%5E%7B%27%7D%20%3D%20I) So here, we need to estimate A0 (which is assumed to be a lower triangular matrix) in order to fully describe the SVAR. One popular method was proposed by Sims (1980) and involves short-run assumptions using the Cholesky decomposition of the variance-covariance matrix such that: ![](http://latex.codecogs.com/gif.latex?%5Csum%20%3D%20PP%5E%7B%27%7D) Where: ![](http://latex.codecogs.com/gif.latex?P%20%3D%20A_%7B0%7D%5E%7B-1%7D) By recursive substitution of the VAR(1) process: ![](http://latex.codecogs.com/gif.latex?Y_%7Bt+j%7D%20%3D%20B_%7B1%7D%5E%7Bj+1%7D%20+%20Pu_%7Bt+j%7D%20+%20B_%7B1%7DPu_%7Bt+j-1%7D%20+%20...%20+%20B_%7B1%7D%5E%7Bj%7DPu_%7Bt%7D) And finally, the impulse-response function of Y_t+j is: ![](http://latex.codecogs.com/gif.latex?%5Cpsi%20_%7Bj%7D%20%3D%20B_%7B1%7D%5E%7Bj%7DP) * * * As you can see, I understand the process completely and I would like to be able to do it using R. So what I'm trying to do is: (i): Estimate the final VAR(p) process (p determined by AIC) (ii): Generate the impulse-response function Thank you very much.


title: Variance and Asymptotic normality of sample variance of normal distribution

text: Let $X_1,X_2$,....,be a random sample from $N(q,w^2)$; $q,w$ are unknown. Let $S_n$ be the sample standard deviation. i.e $S_n^2=\frac{1}{n-1}\sum(X_i-\bar{X})^2$ What is $Var(S_n^2)$? and how to show that $S_n^2$ is asymptotically normal? I tried to do Variance part using moment generating functions but expressions are getting extremely complicated.


title: Speed of linear modeling

text: Asked already on stackexchange, this place is probably a better fit, if somebody could delete the old post I'd appreciate it. I have a linear regression equation from school , which gives a value between 1 and -1 indicative of whether or not a set of data points are close enough to a linear function ![enter image description here](http://i.stack.imgur.com/sgzA5.gif) and the equation given here http://people.hofstra.edu/stefan_waner/realworld/calctopic1/regression.html under best fit of a line. I would like to use these to do simple gesture detection based on a point in 3-space (x,y,z) - forward, back, left, right, up, down. First I would see if they fall on a line in 2 of the 3 dimensions, then I would see if that line's slope approached zero or infinity. Is this fast enough for functional gesture recognition? If not, could someone propose an alternative algorithm?


title: Variable ordering using PCA

text: I've downloaded a script to draw a correlation matrix using colored circles. This script allows to order variables using PCA, but I'm not sure how it works. The code responsible for ordering is below:               if (order) {         if(!n==m){                 stop("The matrix must be squre if order is TRUE!")         }       x.eigen <- eigen(corr)$vectors[, 1:2]       e1 <- x.eigen[, 1]       e2 <- x.eigen[, 2]       alpha <- ifelse(e1 > 0, atan(e2/e1), atan(e2/e1) + pi)       corr <- corr[order(alpha), order(alpha)]     }      **Question** :   What is the interpretation of such ordering and what theory lies behind it?


title: Whether to average two groups and analyse calculations from averages or should all analyses include all baseline data from each individual?

text: I have followed two groups and measured a number of variables. Some of the variables can be used to calculate new measures. Is it ok to average the two groups and analyse calculations from the averages (ie each group acting as a single entity to provide a new measurement) with separate tests or should all the analyses include all the baseline data from each individual?


title: What is a good distribution to model average sales

text: I would like to run a monte carlo simulation for a financial model what needs to predict sales. I have somewhat successfully done this if I assume the average order size is normally distributed (at say 80) and has a standard deviation (of say 30). However, the problem is that the minimum sale of course can't be less than 0. Can someone suggest a good alternative distribution (or technique for working with a normal distribution) to handle this problem (I am sure it must be a common enough scenario!). Thanks.


title: Summary statistics of the precision-recall curve

text: From what I understand, one can use the AUC of the ROC curve as a summary statistic of the full curve. **Q1.** Are there any similar summary statistics that one can use on a single precision-recall curve? **Q2.** As far as I understand, the $F$-score ($F_1$ or $F_\beta$), is measured at a specific operating precision-recall regime. i.e. one needs to fix a point on the precision-recall curve to obtain a precision rate and it's associated recall rate to obtain a result in the formulas: $F_1 = 2 \cdot \frac{\mathrm{precision} \cdot \mathrm{recall}}{\mathrm{precision} + \mathrm{recall}}$ $F_\beta = (1 + \beta^2) \cdot \frac{\mathrm{precision} \cdot \mathrm{recall}}{(\beta^2 \cdot \mathrm{precision}) + \mathrm{recall}}$ If so, how does one choose such point on the curve?


title: Simulating data for logistic regression with a categorical variable

text: I was trying to create some test data for logistic regression and I found this post How to simulate artificial data for logistic regression? It is a nice answer but it creates only continuous variables. What about a categorical variable x3 with 5 levels (A B C D E) associated with y for the same example as in the link?


title: Cox proportional hazard model fit to complex survey data

text: What type goodness-of-fit statistics may be used to evaluate the fit of a cox proportional hazard model to complex survey data?


title: how does one predict edges / links / connections on a weighted directed graph network?

text: Given an adjacency matrix A for a weighted, directed graph (so matrix elements are not just 0/1 and the matrix is not symmetric), are there any good methods for predicting new edges? I have a VERY large (billions of nodes) dataset with known edges for some connections and NULL values for unobserved connections, and I would like to use the known edges to make predictions of the unobserved ones. It is fine if the method is not exact -- in fact, I would like to keep edges NULL or 0 if the prediction is below some threshold, just to keep the matrix as sparse as possible for data size and processing speed. Any thoughts?


title: Bounds for the population variance?

text: Suppose we have i.i.d. samples $x_1$, $\ldots$, $x_n$ for a (potentially non- normal) random variable $X$ with finite moments. We can use these samples to construct an unbiased estimates of the population mean and population variance $$ \bar{x} = n^{-1} \sum_{i=1}^n x_i \qquad\text{and}\qquad s^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2 \enspace. $$ Without making any assumptions on the distribution of $X$, it is possible to construct probabilistic bounds on the population mean, by using Chebyshev's inequality (see, e.g., wikipedia or the original paper). My question is: do such probabilistic bounds exist for the population variance? In other words, can we say that with probability $\delta$ the population variance $\sigma^2$ will be in some interval $[L(\delta,\\{x_i\\}),U(\delta,\\{x_i\\})]$? And if so, what are the functions $L$ and $U$ that describe the lower and upper bound? For normal distributions the sample variance follows a $\sigma^2 \chi^2_{n-1} (n-1)^{-1}$ distribution. This can be used to construct confidence intervals. However, I am looking for more general bounds that apply also to non-normal settings.


title: What is the best estimate for variability of a measurement

text: I am trying to capture the variability of a specific measurement. I have 9 measurements from each patient – 3 measurements at 1-hour intervals over 3 consecutive days. So it kinda looks like this:               Person |  Day1Trial1 Day1Trial2 Day1Trial3 Day2Trial1 ... Day3Trial3        1                     2        ...      I'm not sure how to proceed from here given very little background in statistics. Do I just calculate the overall variance of all the observations? How do I take into account the error of the device I used, or the variability of a patient's measurements? Can someone please refer me to any literature that I can read that's relevant to this? Are there any specific methods or techniques that seem appropriate in this situation?


title: Time series: correcting the standard errors for autocorrelation

text: I have performed a number of tests to detect any presence of autocorrelation in my monthly return series. The test results confirm that the standard errors are not independent. A Durbin-Watson test result shows an upper bound violation with a d-statistics of 2.16, which implicates (first order) negative autocorrelation. A second, Breusch-Godfrey test, performed to examine higher order correlation points outs that for the first 12 lags the tests fails to reject the null of no serial correlation. Isn’t this strange since the test results are contradicting each other? To get a better understanding of the correlation within the (dependent) variable I also implemented a corrgram. The Q-statistic results in this case shows a significantly autocorrelated data after the first lag (see image). The additional independent variable Q-statistic results show no presence autocorrelation. ![Corrgram_Results](http://i.stack.imgur.com/BXG6B.jpg) What I came across so far while searching on the internet for solutions to solve the autocorrelation are a large number of solutions. In order to obtain meaningful results from my OLS-regression I thought it was best to include a lagged dependent variable in the regression and generate Newey-West standard errors. This is probably a relatively simple but accessible approach for inexperienced statisticians, like myself. **My question is;** Is this a correct and sufficient approach of tackling this particular autocorrelation problem in my data, or should I think about more advanced models? A link to the time series file;Time_Series_Data


title: Maximum Likelihood estimator - confidence interval

text: How can I construct an asymptotic confidence interval for a real parameter, starting from the MLE for that parameter?


title: Applications of Item Response Theory

text: As far as I know, item response theory is mainly used within psychology to analyse survey data, personality questionaires or in educational setting.   * Are there other areas (possibly outside of psychology) where IRT is used?    * And if so, what are they?


title: High-Frequency Time-Series Forecast With A Lower Bound

text: I am helping a friend with a data project. He's interested in building a canary-in-the-coal-mine alert system for his website which tells him when the number of users dips below some critical lower bound. The number of users varies by time of day, and day of the week. I have 550 days of by-the-minute data, though we can probably aggregate that into 5 minute buckets to make it less noisy at the cost of early warnings. Suggestion on how to pick the right size bucket would be helpful. This is what the **smoothed** daily fluctuations look like by DOW:![enter image description here](http://i.stack.imgur.com/HPGeR.png) There appear to be two cycles each day, though there's typically not much of dip in between them. There's also a positive trend over time that I need to take into account:![enter image description here](http://i.stack.imgur.com/XBpT2.png) Obviously, the 4am-9pm period in the first graph should look the same over time, but the rest of the curves have been shifting up and becoming steeper. I am having trouble getting the shape of the daily curve right. I think this calls for some sort spline with 4-6 knots, but I am not sure how to integrate this with the time-series aspect. Approaching it from the other direction, I've been trying to fit this into an unobserved components model (UCM) in Stata, but I am lost as to how to specify the trend, seasonality, and cyclicality, and how to treat the exogenous variables. The few things I'be tried either look off or have not converged in a reasonable amount of time. The other issue is how to translate the RMSE of the UCM forecast into an interval of some sort, or at least the left end of one. This is much more complicated than any time-series model I've built before, so any modeling advice would be appreciated.


title: Who invented k-fold cross-validation?

text: I am looking for a reference to the paper where k-fold cross-validation was introduced (rather than just a good academic reference for the subject). Perhaps it is too far back in the mists of time to unambiguously identify the very first paper, so any early papers where the idea was used would be of interest. The earliest I am aware of are P. A. Lachenbruch and M. R. Mickey, “Estimation of error rates in discriminant analysis,” Technometrics, vol. 10, no. 1, pp. 1–12, Feb. 1968. and A. Luntz and V. Brailovsky, “On estimation of characters obtained in statistical procedure of recognition (in Russian),” Techicheskaya Kibernetica, vol. 3, 1969. but as far as I can tell they only cover leave-one-out cross-validation (my technical Russian isn't all it could be ;o).


title: poor prediction performance from (WEKA) classifiers

text: I have an imbalanced 2-class data set made of two parts - model creation data and prediction data. Model creation data is split into train and test data. Model is trained using the training data and tested on the test data, as usual. I do this several times to make sure that the data space is covered fully. The accuracy I get on each of the test data from the model based on the corresponding train data (produced by the cross validation procedure) is quite satisfactory. However, when I apply any of these models to the prediction data the model prediction performance I get is rather poor. I am unsure as to what is going on. The statistical profile of the prediction data is similar to the training/test data and I expect to get similar prediction accuracy. One difference between the train and prediction data is that they are from different, but contiguous, time periods but there no time-related variables in the model. How do I solve this problem? Any help from the community is greatly appreciated. Thanks.


title: improving my variable names

text: Good variable names are: a) short / easy to type, b) easy to remember, c) understandable / communicative. Am I forgetting anything? Consistency is something to look for. The way I would put it is that consistent naming conventions contribute to the qualities above. Consistency contributes to (b) ease of recall and (c) understandability, though other factors are often more important. There is a clear tradeoff between (a) name length / ease of typing (e.g. all lowercase) and (c) understandability. I'm investing a fair bit of thought in these issues because thousands of people are using the data and I hope many will use my code to prepare the data and facilitate some types of analysis. The data, from the Longitudinal Study of Adolescent Health, is broken down into multiple datasets. My first step was to take the 227 variables in the most commonly used dataset, recode them, give them more meaningful names. Original variable names are things like "aid", "s1", "s2", which I renamed "aid2", "age", and "male.is". There are thousands of other variables in the other datasets which may be merged depending on what the researcher's goals are. As long as I'm renaming variables, I want to make them as useful as possible. Here are some of the issues I've considered. So far, I have only used lower- case and avoided using any dashes or underscores, and I've only used periods for one very specific purpose. This has the virtue of simplicity and consistency, and causes no problems for most variables. But as things get more complex I'm tempted to break my consistency. Take, for example, my variable "talkprobmsum", it would be easier to read as "talkProbMSum" or better still "talk.prob.m.sum", but if I'm going to use capital letters or periods to separate words then shouldn't I do it for all variables? Some variables are recorded at more than one time, e.g. the race variables so I appended .is or .ih to indicate whether they come from the in-school or in- home questionnaire. But there are surely some repeats I'm not aware of yet, would it be better to append a reference to the dataset to the name of every variable? I need to group-center and standardize a lot of variables, the way I've done that is by appending .zms meaning z-score by male and by school. Any general or specific thoughts or resources are greatly appreciated. See this repository for some of my code, and descriptive statistics with a list of variable names. I briefly described the reason for sharing this code here, and it was publicized a bit here, but these last two links aren't really relevant to the issue of variable naming conventions. **Added:** I edited this lightly, mostly just moving a paragraph, to try to avoid some of the confusion evident in the comments. Thanks for thoughts!


title: Getting a data frame of logit probabilities and their confidence intervals

text: I have the following model and have used the effects package to plot the predicted probabilities and the confidence interval lines. However, I was wondering how I'd go about spitting out a data frame in R which has the response value, the low and high ci values, and the predicted values. Kind of like the following               mod1 = glm(won_ping ~ our_bid, data=ndat, family=binomial(link="probit"))     summary(mod1)          library(effects)     plot(effect("our_bid", mod1), rescale.axis=FALSE, multiline=TRUE, xlim=c(0,2000), main="129- AH - Bid model")      Desired output with our response variable and the confidence intervals for the predicted probabilities:               our bid    low_ci     hi_ci      prob     1           0.15       0.21       0.17     2           0.18       0.23       0.20      3           0.20       0.30       0.25        4           0.21       0.25       0.23       ...      How would I go about getting the following result in R. I tried the following but it doesn't work as I want the probabilities and not the log-odds values.               > as.data.frame(effect("our_bid", mod1))        our_bid        fit         se       lower      upper     1       25 -2.3549908 0.04598536 -2.44512045 -2.2648612     2      238 -2.0297771 0.03794491 -2.10414781 -1.9554065     3      451 -1.7045635 0.03724233 -1.77755712 -1.6315699     4      664 -1.3793498 0.04422870 -1.46603649 -1.2926632     5      877 -1.0541362 0.05610148 -1.16409305 -0.9441793     6     1090 -0.7289225 0.07043143 -0.86696557 -0.5908795     7     1303 -0.4037089 0.08599888 -0.57226356 -0.2351541     8     1516 -0.0784952 0.10224011 -0.27888213  0.1218917     9     1729  0.2467185 0.11887929  0.01371934  0.4797176     10    1942  0.5719321 0.13577018  0.30582746  0.8380368


title: Is it ever good to increase significance level?

text: Suppose you are using $\alpha = 0.01$ as the significance level. IF you are not getting significant results, is it better to increase the significance level to $\alpha = 0.05$?


title: Treating a variable as categorical or continuous

text: Suppose a variable $T$ (temperature) is coded as follows: $T = 1$ if temperature < 10, $T = 2$ if 10 < temperature < 40, and $T=3$ if temperature > 40. Is it better to treat $T$ as a categorical variable?


title: Gibbs sampler from conditional distribution

text: I am trying to propose Gibbs sampling with the density below, $$p(y_1,y_2,y_3)\propto \exp [-({{y}_{1}}+{{y}_{2}}+{{y}_{3}}+{{\theta}_{12}{y_1}{y_2}}+{{\theta }_{13}{y_3}{y_1}}+{{\theta }_{23}{y_2}{y_3}})]$$ where, $({{y}_{1}},{{y}_{2}},{{y}_{3}})\in R_{+}^{3}$ and ${{\theta }_{ij}}=i+j$ How do I find the full conditional distribution? And then, I'll generate sample $\\{(y_{1}^{i},y_{2}^{i},y_{3}^{i})\\}$ for $i=1,...n$. I understand Gibbs sampling, sample one variable while keep others fixed.


title: Convergence in mean square - missing step

text: My question regards validating the legitimacy of representing an AR(1) as a MA($\infty $) process. In my notes this is done by showing convergence in mean square, saying that: $\mathbb{E}[(\sum_{i=0}^{n-1}\phi^i\epsilon_{t-i}-Y_t)^2]=\mathbb{E}[\phi^{2n}Y_{t-n}^2]$ And then noting that: $\phi^{2n}\gamma_0 \rightarrow 0$ as $n \rightarrow \infty$ The second step is fine. Could anyone show how the first equality is made? I get the feeling this shouldn't be too hard, but I just don't see it. [For context, see 3.5.1 of these notes, which are similar]


title: Can Principal Component Analysis be used alone to infer major patterns within data instead of k means clustering?

text: Do you think it can be used instead of k means? I obtained a correlation with the first 2 components as they carry over 90% of the weight. Would you agree on the technique?


title: Converting J48 to if-then rules in Weka

text: I have a J48 decision tree model trained with WEKA. I would like to access the rules of the tree in J48 so that I can somehow use them in my code whether with if-else statements or as a decision table I can access in my code. Is this possible? If yes, how?


title: How to compare two datasets using metrics drawn from unknown distributions and with small sample sizes?

text: I have two datasets consisting of metrics from several experiments. Dataset 1 is the collection of results of experiments E performed by user A on product A, repeated N times. Dataset 2 is the collection of results of the same experiments E performed by the same user A on product B, repeated the same N times. N is not large, and cannot be large due to practical limitations (typically around 15-20). The data CANNOT be assumed to be Gaussian. In some cases, it is known to be definitely not normal, and in some cases we do not know for sure. It is just an unknown distribution. We know that the metrics cannot be negative. That is pretty much the only definite information we know. Using this data, how do we compare product A and product B and give a result with some statistical significance? How do we devise a hypothesis test to check if A > B with some x% confidence? Edit


title: Example of a 2nd order stationary, but not strictly stationary process

text: Does anybody have a nice example of a stochastic process that is 2nd-order stationary, but is not strictly stationary?


title: Area under curve of ROC vs. overall accuracy

text: I am a little bit confusing about the Area Under Curve (AUC) of ROC and the overall accuracy.   1. Will the AUC be proportional to the overall accuracy? In other words, when we have a larger overall accuracy will we definitely a get larger AUC? Or are they definitely positively correlated?   2. If they are positively correlated, why we bother report both of the in some publications?   3. In real case, I performed some classification task and got the results like: classifier A got an accuracy 85% and AUC of 0.98 and classifier B got an accuracy of 93% and AUC of 0.92. Then which classifier is better? Or will it be possible to get some results like this (I mean there might be a bug in my implementation)?


title: Post-hoc for 2x2 mixed design ANOVA using SPSS

text: I am analysing an experiment run with 20 participants using a $2\times 2$ mixed design ANOVA. The experiment has   * one within subjects variable _A_ with two levels ( _a1_ , _a2_ )    * one between subjects variable _B_ with two levels ( _b1_ , _b2_ ) The following are the results of the ANOVA analysis: > B F(1,18) = 5,84 p<0,026493 SS=1,44 MSe=0,25 > > A F(1,18) = 11,1 p<0,003734 SS=1,25 MSe=0,11 > > B*A F(1,18) = 7,16 p<0,015412 SS=0,81 MSe=0,11 Is there any reason to run post-hoc test in this context? SPSS complains about variables having less than 3 levels and, for this reason, doesn't allow me to run Post-hoc tests. **Question** Is it possible to run bonferroni post-hoc test in this kind of scenario? If not, why? If yes, how to do it in SPSS?


title: Show that $G(x)$ is a distribution function and find mean

text: > Let $F$ be a distribution function on $\mathbb{R}$ with $F(0)=1$ and $\mu$ > be its mean.Show that $$G(x)=\frac{1}{\mu}\int_{0}^{x}[1-F(t)]dt$$ is a > distribution function. Also find its mean. Trial: For $x_1>x_2$ we have $G(x_1)>G(x_2)$. So, it is incrasing. Again $G(\infty)=\frac{1}{\mu}\int_{0}^{\infty}[1-F(t)]dt=\frac{1}{\mu}\cdot\mu=1$ But how I show $G(-\infty)=\frac{1}{\mu}\int_{0}^{-\infty}[1-F(t)]dt=0$ Lastly help me to find mean. Please help.


title: Variance of a time series fitted to an ARIMA model

text: I think this is a basic question, but maybe I am confusing the concepts. Suppose I fit an ARIMA model to a time series using, for example, the function auto.arima() in the R forecast package. The model assumes constant variance. How do I obtain that variance? Is it the variance of the residuals? If I use the model for forecasting, I know that it gives me the conditional mean. I'd like to know the (constant) variance as well. Thank you. Bruno * * * ## Update 1: I added some code below. The variance given by `sigma2` isn't close to the one calculated from the fitted values. I'm still wondering if `sigma2` is the right option. See figure below for time series plot.               demand.train <- c(10.06286, 9.56286, 10.51914, 12.39571, 14.72857, 15.89429, 15.89429, 17.06143,                    17.72857, 16.56286, 14.23000, 15.39571, 13.06286, 15.39571, 15.39571, 16.56286,                   16.21765, 15.93449, 14.74856, 14.46465, 15.38132)     timePoints.train <- c("Q12006", "Q22006", "Q32006", "Q12007", "Q22007", "Q32007", "Q12008", "Q22008",                           "Q32008", "Q12009", "Q22009", "Q32009", "Q12010", "Q22010", "Q32010", "Q12011",                           "Q22011", "Q32011", "Q12012", "Q22012", "Q32012")          plot(1:length(timePoints.train), demand.train, type="o", xaxt="n", ylim=c(0, max(demand.train) + 2),           ylab="Demand", xlab="Quadrimestre")          title(main="Time Series Demand of Product C", font.main=4)     axis(1, at=1:length(timePoints.train), labels=timePoints.train)     box()          ### ARIMA Fit     library(forecast)          # Time series     demandts.freq <- 3     demandts.train <- ts(demand.train, frequency=demandts.freq, start=c(2006, 1))          # Model fitting     demandts.train.arima <- auto.arima(demandts.train, max.p=10, max.q=10, max.P=10, max.Q=10, max.order=10)     print(demandts.train.arima)     summary(demandts.train.arima)     demandts.train.arima.fit <- fitted(demandts.train.arima)          # Forecast ARIMA (conditional means)     demandts.arima.forecast <- forecast(demandts.train.arima, h = 3, level=95)     print(demandts.arima.forecast)          # Constant variance from ARIMA     demandts.arima.var <- demandts.train.arima$sigma2     print(demandts.arima.var)          # Variance from fitted values     print(var(demandts.train.arima.fit))      ![Time Series Plot](http://i.stack.imgur.com/E5gv0.png)


title: Dealing with Dates of Birth in Predictor

text: I have a predictor variable of type factor which contains the date of birth. From your experience (when dealing with logistic regression), what is the best way of treating date of birth (or similar types)? As a continuos variable or maybe or maybe as a categorical variable which stores different age segments?


title: Derivation of prediction intervals for a normally distributed population with unknown population standard deviation

text: I have via the ISO standard 16269 found the solution to a problem that I've been working on. Based on a couple of independent samples from a normally distributed population, I would like to determine a limit above which all the remaining samples in the population will fall with a given confidence level. The standard does not, however, provide sufficient underlying theory for me to get a satisfactory understanding of the statistics. **From the standard follows:** > A random sample of n observations has been drawn from a normally distributed > population with unknown mean $\mu$ and unknown standard deviation $\sigma$. > The sample mean is $\bar{x}$ and the sample standard deviation is $s$. > > For given values of $n$, $m$ and $\alpha$ the smallest factor $k$ is > required such that one may have at least $100(1-\alpha)%$ confidence that > none of $m$ further observations will exceed $\bar{x}+ks$. > > $k$ is given by > $$\int_0^\infty{g\left(s\right)\int_{-\infty}^\infty{\Phi^m\left(\bar{x}+ks\right)f\left(\bar{x}\right)d\bar{x}ds}}\geq1-\alpha > $$ > > where > $$f(\bar{x})=\sqrt{\frac{n}{2\pi}}\exp\left(-\frac{n}{2}\bar{x}^2\right), > -\infty\lt\bar{x}\lt\infty $$ $$ g(s)= > \frac{\nu^{\nu/2}s^{\nu-1}}{2^{(\nu/2)-1}\Gamma(\frac{\nu}{2})}\exp\left(-\nu > s^2/2\right),s\geq0 $$ $$ > \Phi(t)=\int_{-\infty}^t{\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}u^2\right)du} > $$ and $$ \Gamma\left(\frac{\nu}{2}\right)= > \int_0^\infty{x^{\frac{\nu}{2}-1}\exp\left(-x\right)dx}$$ $$ \nu=n-1 $$ I identify the resemblance with the $\chi^2$, originated from the uncertainty of the variance, as well as the normally distributed contribution from the uncertainty of the mean. But I can't get it together. I am looking for directions/guidance to help me get a better understanding for the expression which defines $k$? If no help is to be found here perhaps somebody can recommend a good reference on the topic? I have searched for the references in the standard, unfortunately none are available via my university.


title: "Least square root" fitting? A fitting method with multiple minima

text: I am looking for the name of a fitting method that will work even if points from multiple dataseries are meshed together. As far as I understand there are two major methods, _least squares_ and _least absolute error_ ; the method I am looking for would be the next step in the direction from LS to LAE:   * **Least Squares** fitting: Minimize the square of the distances of all data points to a function. **The further a point is, the stronger it's deviation is weighted**. May be strongly influenced by outliers. Will always converge to one globally unique solution (at least for simple curves), easy to find with gradient descent methods.  ![enter image description here](http://i.stack.imgur.com/E73dX.png) * * *   * **Least absolute deviations** : Minimize the sum of the distances of all datapoints. **All deviations are weighted equal**. More robust to outliers, may converge to potentially infinite solutions that all have the same total deviation ![enter image description here](http://i.stack.imgur.com/XQBrN.png) (generated by hand; images may not be accurate) * * *   * " **Least Square Root** " ??? (Google does not really have many hits for this name). I am looking for a function that may converge towards _multiple local minima_ , and **weights deviations _stronger the closer a point is_**. Which solution it would converge would _depend on the starting parameters_. I think a least-square-root of all deviations would do this, but I cannot find any reference to this sort of fit. ![enter image description here](http://i.stack.imgur.com/mwPOK.png) Is there an official name for a method like this, and what would it be called? Or is this sort of fit fundamentally not possible with curve fitting methods? I am not looking for a method that would find all local minima, one local minimum depending on starting parameters would suffice. * * * (note that in the example pictures the functions were all lines, but ideally I am looking for a method that is able to iteratively fit points towards arbitrary n-dimensional functions)


title: Linear regression with constrained coefficient

text: I have to perform a linear regression of the form               y=a*x1+b*x2+c      where               b=(1-a)      where I have errors both on the dependent variable y and on the predictors x1,x2. I am not sure how to handle the constraint on the second coefficient and the errors on the dependent variable. Any idea?


title: Regression books for business people

text: I am a student of MBA-Business Intelligence and Data Analytics. I would like to know names of some good books on regression analysis written from the perspective of making business decisions.


title: Questions about specifying linear mixed models in R for repeated measures data with additional nesting structure

text: ## Data Structure               > str(data)      'data.frame':   6138 obs. of  10 variables:      $ RT     : int  484 391 422 516 563 531 406 500 516 578 ...      $ ASCORE : num  5.1 4 3.8 2.6 2.7 6.5 4.9 2.9 2.6 7.2 ...      $ HSCORE : num  6 2.1 7.9 1 6.9 8.9 8.2 3.6 1.7 8.6 ...      $ MVMNT  : Factor w/ 2 levels "_Withd","Appr": 2 2 1 1 2 1 2 1 1 2 ...      $ STIM   : Factor w/ 123 levels " arti"," cele",..: 16 23 82 42 105 4 93 9 34 25 ...      $ DRUG   : Factor w/ 2 levels "Inactive","Pharm": 1 1 1 1 1 1 1 1 1 1 ...      $ FULLNSS: Factor w/ 2 levels "Fasted","Fed": 2 2 2 2 2 2 2 2 2 2 ...      $ PATIENT: Factor w/ 25 levels "Subj01","Subj02",..: 1 1 1 1 1 1 1 1 1 1 ...      $ SESSION: Factor w/ 4 levels "Sess1","Sess2",..: 1 1 1 1 1 1 1 1 1 1 ...      $ TRIAL  : Factor w/ 6138 levels "T0001","T0002",..: 1 2 3 4 5 6 7 8 9 10 ...      ## Full Model Candidate               model.loaded.fit <- lmer(RT ~ ASCORE*HSCORE*MVMNT*DRUG*FULLNSS                                   + (1|PATIENT) + (1|SESSION), data, REML = TRUE)        * Reaction times from trials are clustered within sessions, which in turn are clustered within patients   * Each trial can be characterized by two continuous covariates of ASCORE and HSCORE (ranging between 1-9) and by a movement response (withdraw or approach)   * Sessions are characterized by drug intake (placebo or active pharmacon) and by fullness (fasted or pre-fed) ## Modeling and R Syntax? I'm trying to specify an appropriate _full model_ with a loaded mean structure that can be used as a starting point in a top-down model selection strategy. Specific issues:   * Is the syntax correctly specifying the clustering and random effects?   * Beyond syntax, is this model appropriate for the above within-subject design?   * Should the full model specify all interactions of fixed effects, or only the ones that I am really interested in?   * I have not included the STIM factor in the model, which characterizes the specific stimulus type used in a trial, but which I am not interested to estimate in any way - should I specify that as a random factor given it has 123 levels and very few data points per stimulus type?


title: Significance test for highly skewed Bernoulli distribution

text: I am working with two highly skewed Bernoulli distributions where 96-99+% of the samples are in the "false" category, and the rest are in the "true" category (sort of speak). I am looking for a two-sided test of difference of proportions between the two samples. I can often achieve 500+ "trues" and tens or hundreds of thousands of "falses" in a reasonable time but I'm not sure if approximation to the normal distribution can withstand this extreme skewness. I initially thought I might need something non-parametric, but here, I actually know the distribution. I have been using a student's t-test, while paying attention to sample size estimation, but past experience has led me to be skeptical of its results. Thanks for your help.


title: How can I remove a variable from a linear regression?

text: Let's say I have the following multivariable linear regression: y = A **x** \+ _b_ , **x** being a vector of many variables. Is it possible to remove a specific variable _x_ with non-null coefficient in a way the resulting model will be as it had never learned that variable in the first place? I want to be agnostic in relation to the training procedure, so I don't want to retrain the model with one less variable. My suspicion is that if the mean of the variable is zero, just removing it would work, otherwise the bias _b_ would need to be corrected somehow.


title: Help with hypothesis test and generalised likelihood ratio test

text: I have two groups of people undergoing a medical procedure, a relapse group and a non-relapse group. The liklihood function is: $$f_n({\bf x} ; \theta) = \theta^{n_u} \exp[-\theta(\sum_{i \in C} x_i* + \sum_{i \in C^C} x_i)] $$ C is the non-relapse group, C^C is the compliment of that. N_u is the number of patients in C_C. Data has been collected for two types of performing the medical procedure. > x1 <\- c(28, 32, 49, 84, 357, *933, *1078, *1183, *1560, *2114,*2144) > > x2 <\- c(42, 53, 57, 63, 81, 140, 176, -210, 252, *476, 524, *1037) Where stars indicate the patient is in group C. Both types of medical procedure are modelled using the liklihood given. We wish to compare the parameters in the two types of procedure. I am testing the following hypothesis: $$ H_0: \theta_1 = \theta_2 \text{ vs. } H_1:\theta_1 \neq \theta_2$$ Am I correct in thinking I should use the generalised likelihood ratio test? I am unsure how to run this test as I have only ever tested $$ \lambda = \lambda_0 \text{ vs. } \lambda \neq \lambda_0 $$ where $$\lambda = \lambda_0$$ is a specific value. Any help would be appreciated.


title: Leave One Out Cross Validation

text: I tried to implement the Leave One Out Cross Validation (LOOCV) method to get me a best combination of 4 data points to train my model which is of the form: Y= a + b X1 + c X2. Where a, b and c are the coefficients based on regression. I have a set of 20 data points on the whole to train my model but I want to restrict my model to be trained from a set of 4 data points and be able to predict the other data points fairly well. So I used the LOOCV method to achieve that by starting with 19 data points as my training set and the remaining 1 data point as my test set and calculating the error. Then I eliminated the data point (from both training and test data set) that resulted in the minimum error as this data point has been best predicted by the model. Then I kept on iterating until I came to 4 data points. But the problem is, when I tried to add the data point that has been left out from each iteration and compute the Root Mean Square Error Value by applying the model to the whole set of data points at each iteration, I presume that it should have kind of an asymptotic shape when I plot the Root Mean Square Value Vs the No. of data points i.e. higher RMSE at lower no. of data points and converging to the lowest RMSE values at higher no. of data points. But this is the plot I'm getting for the data set. ![enter image description here](http://i.stack.imgur.com/xDTnJ.jpg) I don't know why is the plot showing this type of a behavior. Also, Is it valid to use the Leave One Out Cross Validation method as I have used it?. Please, clarify.


title: b1 estimator properties in intercept-only model

text: For the model $Y_i=B_1+u_i,i = 1, 2, …, n$ (This model does NOT have an independent variable), how do you:   * Derive the OLS estimator of $B_1$   * Prove that $B_1$ is a linear estimator.   * Prove that $B_1$ is an unbiased estimator.    * Derive an expression for the variance of $B_1$   * Prove that the mean value of the residuals for this model is zero.


title: Can normal GARCH innovations outperform student-t GARCH innovations?

text: I am trying to model gold and oil volatility but I am getting quite confusing results. I am using the rugarch package and using the standard GARCH and EGARCH to model volatility. My results seems to indicate that the standard GARCH model with normal innovations is the best but I am reluctant to believe that as descriptive statistics indeed indicate that the returns are not normally distributed. What I find confusing also is that when using t-distribution instead of normal distribution the p-value for my correlation test decreases. I was expecting to see the EGARCH outperform by a large margin but this doesn't seem to be the case. Hopefully someone will be able to clear that up for me. Here is a summary of my results.![Summary](http://i.stack.imgur.com/XGuY9.png) Can someone explain to me why I'm getting such inconclusive results?


title: Combination of variational methods and empirical Bayes

text: Suppose I have a posterior $p(z, \theta | y, \eta)$ with $y$ observed data, $z$ are hidden variables and $\theta$ are parameters, and $\eta$ is a vector of hyperparameters. I construct a mean field approximation to the posterior using coordinate ascent, i.e. $$ q(z) \gets \exp\left\\{E_q(\log p(z, \theta | y, \eta) | z) + \mbox{const}\right\\} \\\ q(\theta) \gets \exp\left\\{E_q(\log p(z, \theta | y, \eta) | \theta) + \mbox{const}\right\\} $$ where $\mbox{const}$ makes each integrate to $1$. Iterate until convergence. My question is, if I want to do empirical Bayes, and (approximately) profile out $\eta$, is it valid to just augment this with an EM step $$ \eta \gets \arg \max_\eta E_q \log(p(z, \theta | y, \eta))? $$ Based on the evidence lower bound $$ \log p(y | \eta) \ge E_q \log(p(z, \theta | y, \eta)) - E_q \log(q(z, \theta)) $$ it seems like I should be able to get away with this; optimizing the lower bound over $\eta$ is exactly the proposed EM step, so I'm still doing coordinate ascent on the lower bound. I suspect that, if this works, it is the obvious thing to do for people who do variational inference regularly. Since variational methods are fast, my thought is that maybe I could do this to set hyperparameters before launching into a Gibbs sampler for exact (up to MCMC error) inference. **Update:** I've toyed around with this a bit and found that adding the EM step can work okay, but in some situations it seems to increase sensitivity to bad local optima in the variational algorithm. Initializations of variational parameters that normally seem to work okay don't. If I instead only throw in the EM step every few iterations it works better, but then of course I slow down convergence. One approach that works okay is to "burn in" the variational algorithm using a method that I trust, and then add the EM step every iteration thereafter. I'm adding a bounty because I'm still looking for some good general advice.


title: Why is it valid to account for k-1 intercepts w/ only 1 random intercept parameter?

text: I possess a basic understanding of random vs. fixed effects, and how to code random effects models in SAS. However, I'm having trouble wrapping my head around the derivation of random effects terms, and how a random intercept model, for example, can describe the variation in $k$ intercepts with a single parameter ($\sigma^2$ for a normal dbn) rather than $k-1$ parameters, which can amount to a huge saving in degrees of freedom. Isn't that cheating? ;) The old school technique would be to use the maximum likelihood method to solve for $k-1$ parameters for each category.   * How does a random intercept model avoid this by using a single parameter?    * Isn't the end result the same--both random and fixed effect models will estimate $k-1$ intercept terms?


title: PCA demeaning the data

text: What is the motivation for demeaning the data when doing PCA. I've been told to do it, but I've never heard a good and/or intuitive reason for it. Is this a case where doing it just makes the math easier, or is it more like: if you don't do it, you'll get the wrong answer. This is something I've been wondering about for a while. One other question: when you demean the data, do you also need to normalize it (i.e., divide by standard deviation)?


title: Centering variables before running PCA

text: I am learning about PCA, regarding PCA I need to know that given a dataset is it always necessary to use centering? what if I don't center the variables used in PCA?


title: Hidden Markov Model speech recognition

text: What if have N audio signals of different length and for each we know corresponding sentence. How we can train classifer for audio signals->sentence ? And maybe Hidden Markov Model can be used here? More generally my question is about, is there some technique exist that can be trained with audio signal not separated to letters or phonemes by hand but just by providing sentence that contains in it,the same with handwriten strings or it can be analogy that we give classifer image and information what object present on image but not specify rect positions of them


title: Are there statistical techniques that investigate such relationships ...?

text: If we have data set, X and Y variable. Say, we do correlation analysis and get some correlation coefficient. Besides, we find an important fact after observing their relationship: That is, the scatter plot of X and Y has a triangular shape. Which means that (for example) when X values are increasing Y values are increasing for all X values, whereas vice versa is not true; when the Y value is increasing X values are anything for all Y values. What kind of analysis should I do to investigate this? * * * ( _Update in response to_ @Penguin_Knight) Your graph is exactly what I explained. Just take away the negative values with an imaginary y=0 line. As you can see there are many data points the X of which is either 0 or very small, and Y has pretty high value. However there is no data point that have y=0 and x is quite high value. And as you can see it makes the data scatter shape a right angle triangle. Thus we could say that the high X values necessitates high Y values but high Y values does not always have high X values. I find it very interesting. In practice for example I measure Complexity and Size of given entities. And my analysis show that all complex entities have big size but not all big-size entities are complex. Then I conclude that the certain amount of complexity requires defined amount of size. You cannot put more complexity in a given size. This is a bit abstract but you see my point? It is very interesting and I would like to get some help on how this kind of relationships are discussed or described in statistics formally.


title: need name, reference, and/or study for the following variable reduction procedure in regression

text: I have seen the following commonly used: 1. fit a model with all variables, 2. in a single reduction step, remove from the model all variables at once that do not fit some criteria (p-value, whatever), 3. calibrate the reduced model, in my case to a new data set, 4. check model results and hopefully everything went well and you can stop. I'm betting this will perform better than stepwise, especially when used with data splitting, but would appreciate a name for this procedure if it exists and perhaps a reference related to it so that I can learn more. Maybe it is so simple/bad/obvious that no one has bothered? -WVG


title: Sample from Wishart distribution with inverse Scale matrix

text: I tried to model precision matrix in a hierarchical Bayesian setup with Wishart prior given d.f. and inverse scale matrix, and matrix normal likelihood, since it's a conjugate prior, my posterior on the precision matrix $K$ ends up in the form: $$ K \sim \text{Wishart}(\text{df}, \Lambda^{-1}) $$ Since the dimension of $\Lambda$ is quite big, I do not wish to take the inverse of the matrix and use the built-in sampler. I looked at the source code for Matlab's `wishrnd`, they used Bartlett decomposition for large dimensions (81+n), which should also work for smaller dimensions albeit its inefficiency. snips of the code:               [n,m] = size(Lambda);     [d,p] = cholcov(Lambda,1);          % Otherwise use the Smith & Hocking procedure     d = eye(size(d)) / d;      % Load diagonal elements with square root of chi-square variates     a = diag(sqrt(chi2rnd(df-(0:n-1))));          % Load upper triangle with independent normal (0, 1) variates     a(itriu(n)) = randn(n*(n-1)/2,1);     % Desired matrix is D'(A'A)D     x = a(:,1:size(d,1))*d;          a = x' * x;          % --------- get indices of upper triangle of p-by-p matrix     function d=itriu(p)          d=ones(p*(p-1)/2,1);     d(1+cumsum(0:p-2))=p+1:-1:3;     d = cumsum(d);      In order to accommodate my needs,i.e. inputing degrees of freedom $df$ and the inverse scale matrix $\Lambda$ I added               d = d \ eye(size(d))      after the second line, so that I would have the inverse of the cholesky decomp of the inverse scale matrix, in others words, the cholesky decomp of the scale matrix. Then everything should be okay from there ( I hoped). I tested this by firstly generating a 3 by 3 positive definite matrix:               >> test = rand(190, 3); Lambda = test'*test;     >> Lambda          Lambda =            62.7642   46.4970   45.6662       46.4970   61.9178   45.4114       45.6662   45.4114   59.1070b      Setting $df = 6$, 100,000 randoms samples were generated and the empirical mean is compare with the first moment of the distribution               >> df = 6;     >> sam = 0; for loop = 1:100000 sam = sam + mywishrnd(Lambda,df); end     >> sam/100000     ans =           0.0956   -0.1069   -0.0689     -0.1069    0.3380   -0.0286     -0.0689   -0.0286    0.3872          >> inv(Lambda)*df          ans =           0.2647   -0.1118   -0.1187     -0.1118    0.2692   -0.1205     -0.1187   -0.1205    0.2857      There is quite a big differences in some entries of the results, but I don't see anything wrong with this theoretically. Is it possibly due to numerical error? Could someone enlighten me on this problem? Thank you!


title: How to perform an Ancova with a covariate with nested factors and unbalanced?Is this analysis appropiated?

text: I'm trying to see differences between years (2008,2009) and between seasons (winter, spring, summer) for a quantitative variable (dependent) in this case is the number of eggs. This variable covariates with a continuous variable, as is the size of the mother (covariate). To view temporal effects during the breeding season, covering 3 seasons, I need to remove the effect of female size. A ANCOVA would be the choice, however, lack a level of factor season, summer in 2009, because it could not perform the sampling. The data are approx. 30 individuals per season. **Sample size by factor level**                             "Year"               "Season"  2008 2009           Winter  31    30           Spring  31    32           Summer  25    NA      There are normality and homogeneity of variance between seasons, years and season: year. I use R for statistical analysis. My issue is that I have applied ANCOVAs, and no differences were found between years, but there were significative differences between seasons and season:year. And I'm not sure if it is correct the type of analysis and as I applied. If I use Type II or III adjusted SS or other type of analysis. Other possibilities? LmFec<-lm(log(eggs)~ log(size) + fseason + fyear+ fseason:fYear,data=b2) (Final model) summary(aov(LmFec))                              Df Sum Sq Mean Sq  F value   Pr(>F)         log(size)       1 17.136  17.136 1217.228  < 2e-16 ***     fseason         2  0.766   0.383   27.196 9.78e-11 ***     fYear           1  0.004   0.004    0.314   0.5764         fseason:fYear   1  0.085   0.085    6.008   0.0154 *       Residuals     143  2.013   0.014                           ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Thanks in advance.


title: Why does the accuracy of the central limit theorem for a mean depends on the skewness of the random variables being summed?

text: > The accuracy of the central limit theorem for a mean depends on the skewness > of the random variables being summed. Could someone explain that to me? Or could someone propose a book to read? Thank you.


title: T test on Proportions vs Chi Square, 3 groups

text: Can some one help me decide which approach is better? I have three groups, each with responses of Yes or No. I'm trying to pick the BEST out of 3 groups Should I use 3 t tests with a bonferroni adjustment? Or use a chi squared test. How do you pick the best group in a Chi square test of multiple groups Thanks


title: Difficult interpreting linear mixed model result - R lme function

text: I'm fitting an harmonic regression model on data from different plants separately as follows:               d1 <- c(0.03895199, 0.04048451, 0.04514816, 0.03701958, 0.04196877, 0.04325377)     d2 <- c(0.05928356, 0.05418038, 0.05917998, 0.06842356, 0.04053254, 0.05161718)     time=c(1/6,2/6,3/6,4/6,5/6,6/6)     curve <- data.frame(time=c(1/6,2/6,3/6,4/6,5/6,6/6), intensity=d2)     STlm <- lm(intensity ~ sin(2*pi*time)+cos(2*pi*time)   +sin(4*pi*time)+cos(4*pi*time),data=curve)     plot(intensity~time,data=curve)     points(intensity~time,data=curve, pch=20, lwd=5)      lines(intensity~time,data=curve, lty=2)     lines(STlm$fitted~time, col=2, data=curve)      The fitting seems good as can be see by the pictures below ![enter image description here](http://i.stack.imgur.com/0Laf3.png) ![enter image description here](http://i.stack.imgur.com/HOMWH.png) But I want to fit a single model. For what I read until now it seems that I have to consider a random effect sampling from this two plant,and following this reference I tried to fit a mixed model in the following way               data <- data.frame(int= c(d1, d2), time=rep(time,2), idx=rep(1:2,each=6))     fm1.lme <- lme(int ~  sin(2*pi*time)+cos(2*pi*time) +sin(4*pi*time)+cos(4*pi*time), data = data, random = ~ 1 | idx)      # model               0.0483369985+0.0007153365 *sin(2*pi*data$time)+  -0.0031904639*cos(2*pi*data$time)+0.0038267490*sin(4*pi*data$time)+0.0014627727* cos(4*pi*data$time)       I read a result that seems reasonable ![enter image description here](http://i.stack.imgur.com/lSmHI.png) I have some initial doubts: Have I used the correct codification to say that the "int" variable (coming from different plants) is the random effect? In the summary:               Linear mixed-effects model fit by REML      Data: data          AIC       BIC   logLik       -23.00751 -23.38614 18.50376          Random effects:      Formula: ~1 | idx             (Intercept)   Residual     StdDev: 0.009743065 0.00723739          Fixed effects: int ~ sin(2 * pi * time) + cos(2 * pi * time) + sin(4 * pi *      time) + cos(4 * pi * time)                              Value   Std.Error DF   t-value p-value     (Intercept)         0.04833700 0.007199211  6  6.714207  0.0005     sin(2 * pi * time)  0.00071534 0.002954652  6  0.242105  0.8168     cos(2 * pi * time) -0.00319046 0.002954652  6 -1.079810  0.3217     sin(4 * pi * time)  0.00382675 0.002954652  6  1.295161  0.2429     cos(4 * pi * time)  0.00146277 0.002954652  6  0.495074  0.6382      Correlation:                         (Intr) s(2*p*t c(2*p*t s(4*p*t     sin(2 * pi * time) 0                                  cos(2 * pi * time) 0      0                           sin(4 * pi * time) 0      0       0     cos(4 * pi * time) 0      0       0       0          Standardized Within-Group Residuals:             Min          Q1         Med          Q3         Max     -1.14446419 -0.31620917  0.01160475  0.39646718  1.37276892          Number of Observations: 12     Number of Groups: 2       Why do I have intercept and coefficient for the random and fixed models, it is not a single model? What is a good measure of fitting, similar to R2 for a linear model? I also have fitted values for "idx" variable, that I used only to inform that I have two different groups                 fm1.lme$fitted              fixed        idx       1  0.04994394 0.04335105       2  0.04650628 0.03991339       3  0.05299024 0.04639735       4  0.05189541 0.04530252      Are there another way to inform that it does not have meaning for the model? Thanks for the help in advance. Cheers


title: How to obtain the elasticity from a log-level regression

text: I have forgotten my undergrad econometrics, and was hoping that someone could help to refresh my memory. If I have a regression $\ln y=a+bX+e$, and want to evaluate the elasticity of $y$ with respect to $X$ (at the mean of $X$), how do I do this?


title: Book recommendation: sample size determination for hypothesis testing of the mean

text: After a previous attempt to get answers here to what I found out was a nonsensical question I am looking for some good resources I can use to begin to understand the area. The majority of Psychology experiment books in our library focus on statistics, data analysis, qualitative research methodologies and experiment design. But methods and techniques for determining sample size, with practical examples, seem to be few and far between, especially at a beginner level. I have found a few books that I believe would be useful online but they are quite expensive and I would rather get feedback that they were worth the money and that there was nothing better out there.   * Sample Size Determination and Power - Thomas P. Ryan   * Sample Size Calculations: Practical Methods for Engineers and Scientists. Paul Mathews So can anyone recommend books or other resources in the area, specifically for determining the sample size for hypothesis testing of the mean? * * * A brief overview of my research area: I am undertaking research to determine which particular features of a website or software interface affect a users perception using metrics such as trust, fairness, expertise etc. Each participant will rate 9 separate designs, each with a different distortion which I will then compare to the control. Participants will be crowdsourced online. Hope this helps.


title: Why does my bootstrap interval have terrible coverage?

text: I wanted to do a class demonstration where I compare a t-interval to a bootstrap interval and calculate the coverage probability of both. I wanted the data to come from a skewed distribution so I chose to generate the data as `exp(rnorm(10, 0, 2)) + 1`, a sample of size 10 from a shifted lognormal. I wrote a script to draw 1000 samples and, for each sample, calculate both a 95% t-interval and a 95% bootstrap percentile interval based on 1000 replicates. When I run the script, both methods give very similar intervals and both have coverage probability of 50-60%. I was surprised because I thought the bootstrap interval would be better. My question is, have I   * made a mistake in the code?   * made a mistake in calculating the intervals?   * made a mistake by expecting the bootstrap interval to have better coverage properties? Also, is there a way to construct a more reliable CI in this situation?                tCI.total <- 0      bootCI.total <- 0      m <- 10 # sample size      true.mean <- exp(2) + 1          for (i in 1:1000){      samp <- exp(rnorm(m,0,2)) + 1      tCI <- mean(samp) + c(1,-1)*qt(0.025,df=9)*sd(samp)/sqrt(10)           boot.means <- rep(0,1000)      for (j in 1:1000) boot.means[j] <- mean(sample(samp,m,replace=T))      bootCI <- sort(boot.means)[c(0.025*length(boot.means), 0.975*length(boot.means))]           if (true.mean > min(tCI) & true.mean < max(tCI)) tCI.total <- tCI.total + 1      if (true.mean > min(bootCI) & true.mean < max(bootCI)) bootCI.total <- bootCI.total + 1      }     tCI.total/1000     # estimate of t interval coverage probability     bootCI.total/1000  # estimate of bootstrap interval coverage probability


title: Validating sign of mean

text: I have a lot of experiments with the following kind of data. The sample size is $n\geq 1$. For each test $i$ in the sample a cost $c>0$ and a revenue $r_i\geq 0$. The cost $c$ is constant. The revenue $r_i$ assumes a few different values and is must often $0$. The calculate the mean income as $\frac{1}{n}\cdot\sum_{i=1}^{n} (r_i-c)$. The sample size increases over time. At any given time, I want to determine if the sample size is large enough to determine if the mean income in positive or negative.


title: Differntiate between point estimate and interval estimate in terms of assumptions for making such estimates?

text: interval estimate is common. However, point estimate appears to be based on simple theorem of central limit theorem.


title: Pooled regression and random and fixed effects

text: Am researching on poverty and growth. I built a model for poverty which includes 11 variables theoretically affecting it. These are inflation, agricultural growth, gini, GDP growth rate etc. My data had few variables in growth rate and others in percentage and index. Poverty and gini actually are indexed values. I somehow managed to transform all of them in growth rate. Authors suggested these variables that is why I narrowed down my variables to those but my r sq is too low less than 1 for instance when I run the pooled regression, f test proves to be insignicant as well. I am really confused. Am using stata 10 someone please comment on this.


title: Probability that the square of a random integer ends in 1

text: From my intro to probability theory book: What is the probability that a randomly selected positive integer will give a number ending in 1 if it is squared? Hint: It is enough to consider the 1 digit numbers. Answer: 0.2 I do not understand why we only need to consider one digit numbers. I see that 1 & 9 are the one digit numbers that produces a square ending in 1, giving a probability of 2/10 of randomly selecting them, but how can we generalize this to all positive integers?


title: Is there a statistically good way to split a set of data into (generally) uneven groups?

text: So here's the scenario: I have about 100 or so "items", each assigned with a number and total guesses... this number is a guess of what people think the item is worth (in terms of money). This number can be an average of the guesses, or the sum of all the guesses. I want to "split" them into (generally) _n_ uneven groups... so for example, it may turn out that some statistical formula/equation/etc. finds out that item 1, 3, 10, etc. are in group 1, 4, 5, etc. are in group 2, and so on (up to group n). Is there any statistical formula/equation/method that can achieve this? I am not sure, but can Discrete Fourier Transforms do this? I've never learned it before but heard that its used to transform data into sinusoids... Example: You have 1000 balls, each with are colored a unique color. You ask random people to look at a some of the balls, and rate it from a number of 1 (bad) to 1000 (good) how much they like it. At the end of this survey, you know that each ball has their own average rating (not necessarily unique). Now, you want to split up the balls in such a way that the balls are partitioned into, say, 5 groups (1 = bad, ..., 5 = good). But here's the thing: you can't just split them into 5 equal groups; it could be the case that; what if there are balls that have a big number difference (i.e. ball 851 has an average vote of 700, but ball 852 has an average vote of 800?)? The votes are basically denser in some areas than others. I think this is the basic idea... I'm not sure what else to add. The partition of the data **must** be consistent and reproducible.


title: Determining a correlation factor and applying it to a second set of numbers

text: > **Possible Duplicate:**   >  Determining a correction factor and applying it to a second set of numbers 1st- I would like to figure out the correlation/ correction factor between two sets of numbers (time series). 2nd- I want to apply that correction factor to a second set of numbers. As an example- say I have two temperature sensors, and I place them in the same exact environment. Even though they should be reading the same temperature, I can expect a little bit of offset due to error between the sensors. Something like:               Sensor 1 Sensor 2     10            10     10             9     11            11     12            11.8     13            12.9     14            13.9      So start, I'm not sure about the best way to figure out the correction factor. It's pretty easy to run a regression analysis or a correlation analysis by plotting an x y scatter plot and calculating a best fit line. But, I'm not sure if this is the way to go. But the real question is that once I figure out this correction factor, how do I apply it to a second set of numbers in a way that will reflect the error associated with the first set of numbers (basically that sensor 2 is reading slightly lower than sensor 1). To continue the above example, say I take the same two sensors, but now place them in different environments where they will be exposed to different temperatures. Now I have a second set of numbers (below) from the same instruments but no longer reading the same temperature. How do I relate the original correction factor to a set of numbers from different environments.               Sensor 1 Sensor 2     10            13     11            14     12           14     12           14     13           14      Since the sensors are now in two different environments, I can no longer just plug numbers into a y=mx+b type linear regression equation. So how can I account for the known offset between two sensors? For the first set of numbers the regression analysis yields the equation y=1.0925x -1.3125. My thoughts were that I could use the regression equation from the first set of numbers as my correction factor, and then apply that to the second set of numbers in order to adjust them to account for the inherent error of the sensors. I was hoping to adjust the second set of data by adding the y intercept value to sensor 2, but this obviously does not work. So I think I have been barking up the wrong tree


title: Why do you maximize specificity for circumstantial findings?

text: Tests for a medical condition can be classified into two categories:   1. Tests on an asymptomatic (no reason to suspect condition) patient where the result is a circumstantial finding; and   2. Tests on a symptomatic (prior reason to suspect condition) patient where the result is a diagnostic finding. In the first case, the emphasis of any test should be on maximum specificity. For the second case, the emphasis of any test should be on maximum sensitivity. Intuitively, I understand why. However, can somebody give me a formal run- through, preferably from a Bayesian perspective, as to why this is the case? References: Slide presentation, "Categorizing variants after whole genome sequencing", J.S. Berg Lang, E., Naraghi, R. (2005) Neurovascular relationship at the trigeminal root entry zone in persistent idiopathic facial pain: findings from MRI 3D visualisation. _J Neurol Neurosurg Psychiatry 76_ :1506-1509. Full text


title: Is a g-test appropriate for count data?

text: I am analyzing an forest ecology experiment where we counted the number of trees in 5 pairs of plots in a forest. One member of each pair was fenced several years ago to exclude deer (exclusion treatment) and the other left so that deer could access it (control). Within a given pair the 2 plots are adjacent to each other and in similar environmental conditions, but the number, size and species of trees in each plot when the experiment began is slightly different. We are interested in whether excluding deer increases the abundance of trees in the forest. We are particularly interested in whether deer have _decreased_ the number of small trees, either by killing small trees directly or reducing the number of seedlings and sapling that survive to the "small tree" stage. (We consider a tree to be anything taller then 2 meters in height, regardless of how thick in diameter the stem is; a small trees is anything less than 2 cm thick). We have not, however, followed the fates of individual trees since the experiment began. What we have are **counts** of the number of trees in each of the ten plots and the size of each tree. The data, pooled across all plots for each of the two treatments and organized into binned size classes, is below. Because deer only affect the fate of small trees, we would expect that there would be more small trees in our fenced plots. The number of large trees, however, should not be affected because the experiment has not been run long enough to result in a change in the number of large trees. We are therefore hypothesizing that the overall size distribution of trees should be different between the treatments, but that the difference should be driven by smallest size classes. It has been suggested by my collaborators that I analyze the data as a contingency table using a g-test. An undergraduate working on this project had previously applied a g-test to each _row_ of the table and reported a treatment effect for the two smallest size classes but not the others. This is not, however, a valid use of a g-test because, as pointed out when I posted this initially, each tree is not an independent data point. I have the following questions about how to analyze this data correctly: 1\. Would it be correct to generate a table like the one below for _each_ of the five pairs of plots and apply a goodness of fit test like a g-test? Unfortunately, when broken down into the individual plots many cells of the contingency tables are 0 or < 5\. 2\. To do such an analysis, do I need to assume an underlying distribution (eg poisson) of the size classes to generate predicted values? 3\. Would it make more sense to analyze this type of count data with a generalized linear model with a poisson distribution? I believe this would require analyzing each size class separately. My collaborators would like to keep things as simple as possibly, but I am perfectly comfortable with running a GLM. Again, the data is fairly sparse, which might make this type of analysis difficult. **DATA:**                    SIZE_CLASS CONTROL_trt EXCLOSE_trt       1   1_or_less          15          60       2      1_to_2          29          30       3      2_to_3          11          10       4      3_to_5          11          15       5 5_or_larger          45          50


title: Choosing the best model to predict second hand widget prices

text: I am looking for some thoughts on how best to approach a dataset and pick the right model for predicting the second hand price of widgets. ## Some background on the properties of these widgets:   * The widgets come in different models, sizes, colors and a few other properties that are categorizable, along with the level of wear & tear that is unique for each widget (but can be measured accurately).   * The value of widgets depreciates over time (it depreciation profile is similar to 1/x) up until some fixed value which is equal to their scrap value. The scrap value is not constant over time but > 0.    * The widgets’ prices are observed once in a time period and recorded along with other properties for each widget. The price is not updated again.    * There is constant stream of new widgets hitting the market and the data comes pretty uniformly so the number of new widgets coming to the second hand market is pretty stable at say some k widgets per day s.t. the stock in time t+10 would be stock at time t plus 10*k. ## Questions on how to handle the data:   * Because the price is observed only once, what should be the approach to setting upper limit as to how old the price data can be to be taken into consideration in the analysis? What kind of preliminary analysis should be used to determine such cutoffs?    * Also, should there be any weight put on the recency of the price data point in modeling the price? ## How to choose the best model?   * What things should be considered about the underlying data before choosing the right model (treating outliers etc)?   * Is multiple variable regression the right way to approach this data or is there already an establish method of approaching such cases of price prediction?


title: Linear regression vs analysis of variance: how to explain the difference?

text: On the Mathematics site an OP who is just learning statistics gave his description of the difference between linear regression and ANOVA and asked if his interpretation was correct. I responded that linear regression considers how a set of covariates relate to a response in a functional form (could have added "that is linear in the parameters") whereas ANOVA categorizes the response into a class or classes of group(s) and tests for a difference between group means. A member downvoted my answer saying that ANOVA can include continuous predictors as well. His own answer indicated that he was considering the term ANOVA to mean the testing of significant terms from the decomposition of variance in the general linear model. We both gave descriptions of linear regression that agreed. My question is: "What do you think is the best answer? His answer? My answer? An explanation providing the two meanings of ANOVA? Something else? Link: http://math.stackexchange.com/questions/183704


title: What is the difference between linear regression on y with x and x with y?

text: The Pearson correlation coefficient of x and y is the same, whether you compute pearson(x, y) or pearson(y, x). This suggests that doing a linear regression of y given x or x given y should be the same, but I don't think that's the case. Can someone shed light on when the relationship is not symmetric, and how that relates to the Pearson correlation coefficient (which I always think of as summarizing the best fit line)?


title: how to identify the regression equation of Y on X and X on Y?

text: I have one general question. Suppose you are given two regression lines of the variables x and y: Say- 4x-6y+50=0 and 15x-6y-120=0. Then how to identify among these two regression equations, the one which is Y on X and the one which is X on Y?


title: Permutation for large sample sizes (and z-test)

text: I was trying to do a permutation test on a large amount of temperature observations in R. Approximately 1700 temperature observations (`g1`) from one place and 2 times 1700 observations from two other places (`g2`). I would like to test the hull hypothesis $\mu_1-\mu_2 = 0$ against the alternative hypothesis $\mu_1-\mu_2 < 0$. I am doing the following calculations:               k            <– 10000     perms        <- replicate(k, sample(c(g1, g2)))     my.mean.difs <- apply(perms[1:length(g1), ], 2, mean) -                     apply(perms[(length(g1)+1):(length(g1)+length(g2)), ], 2, mean)          p            <- sum(my.mean.difs < (mean(g1) - mean(g2))) / k      The p-value, `p`, gives me 0, straight. I don't understand why I am not getting any difference in the means (`my.mean.difs`) which is smaller than the difference in the actual mean. My immediate guess is that I am using too few resamples, but since the data sets are so big (1700 and 3400 data points) it will be very time consuming to do many more resamples. I could randomly create some smaller data sets, and test them, but I am not sure that is a good procedure. Can anybody explain me why I am getting a p-value of 0, straight, and if it is a bad idea to do permutation on big data sets like this? ## Bonus info: ### Temperatures in group 1 (`g1`)   * $\mu_1 = 24.14$ °C    * $\sigma_1 = 0.96$ ### Temperatures in group 2 (`g2`)   * $\mu_2 = 24.95$ °C    * $\sigma_2 = 1.01$ ## Edit: Based on @whuber's comment and the central limit theorem I tried to make a two-sample Z-test of the means, which should be OK to use here, if the assumption of equal variances is true, right? I did the following calculations:               (mean(g1)-mean(g2)) / sqrt(var(g1)/length(g1) + var(g2)/length(g2))      This gives me an extremely small Z-value of -28.25. Is this magnitude of the Z-value realistic? Sorry for the probably very basic question, but I am quite new to statistics, and from my intro stat course, we worked on much smaller sample sizes, and Z-values much closer to 0. But I guess that the large sample sizes makes it possible to show a clearly significant difference in the means, even if the actual difference not is that big. It that right?


title: How to determine the best split threshold (split point) in Hoeffding tree (for continuous attributes)?

text: In Hoeffding tree (that handles continuous attributes), Hoeffding inequity suggests the minimum amount of instance to be looked at before the best split attribute can be determined with confidence. What I don't get is when the best split attribute is decided, how is it guaranteed that the split threhsold (or split point) for the chosen attribute, based on the limited amount of instance, is also optimal?


title: Random intercepts as response variables: Is there a name for this method?

text: I'm trying to find the name of this method (and ultimately a reference). The approach is as follows: 1) Fit a mixed-effect model with a random intercept $$ E(Y_{ij})= \beta_0+\beta_1x_1+\gamma_{0j} $$ 2) Use the estimated random intercepts as the response variables in a new linear regression $$ E(\gamma_{0j})= \beta_2+\beta_3z_3 $$ Any ideas?


title: GLMNET in R: Difference between DFMAX and PMAX

text: In `glmnet` in R, what is the difference between `dfmax` and `pmax` when creating a model with both ridge and lasso penalties? Is the following interpretation correct? Let's say we have we run a `glmnet` elastic net with 500 predictors and we only want to include in 3 predictors in the final model; i.e. we want the lasso eliminate 497 predictors. In order to do so, do we set `dfmax = 500` and `pmax = 3`?


title: I am unsure how to execute and interpret a simple linear regression model where Y ~ a single X

text: I have the following table of X and Y values:                  X     Y     1    68   55     2    54   38     3    90  95     4    64  63     5    61  58     6    51  40     7    79  74     8    51  32     9    83  84      10  48  45      I have to do the following:   * Make linear regression using the method of least squares   * Draw some conclusions from this model when X = 100   * Calculate the standard error for the model


title: Two-way interaction: interpret output from LMM in lmer

text: I have asked a question in this topic relating interpreting two- and three-way interactions in this link: How to interpret 2-way and 3-way interaction in lmer?               DV is height. Baselines are here (all are categorical):           •animal: lion      •color: white      •sex: female      That is for three-way interaction. However, I wondered if my model has only two levels:               Fixed effects:                                       Estimate Std. Error t value     (Intercept)                       164.6888     7.8180  21.065     rat                               -14.1342     8.2889  -1.705     sexmale                           -16.0883    10.0071  -1.608     colorred                            0.5776     6.2473   0.092     coloryellow                        -14.4048     6.1025  -2.360     rat:sexmale                         15.3645    11.8567   1.296     rat:colorred                        12.5258     4.4028   2.845     rat:coloryellow                     10.3136     4.3196   2.388     sexmale:colorred                     2.0272     5.2773   0.384     sexmale:coloryellow                  5.7643     5.1669   1.116      I wonder if the interpretation of two-level will be the same as in three-way interaction. For example, rat:sexmale = I should interpret as "The difference of height between rat and lion is higher for males than for females." This should mean that when we compare height of male and females, males will have higher significant difference between rat and lion while in females, this differences is not significant. Is my interpretation right?


title: Can I estimate Variance of Gamma from Negative Binomial distribution distributed data, given NB is Gamma-Poisson compound

text: I believe the data I have follows Negative Binomial distribution (over- dispersed Poisson). We know Negative Binomial is a Gamma-Poisson compound distribution. The variance of this Gamma distribution is actually the variance I am interested in. Can I estimate the variance of the Gamma distribution of the compound from just my Negative Binomial data. if $X|\lambda ~ \sim Poisson (\lambda)$ $\lambda \sim Gamma(shape = r, scale = p/(1 − p))$ then, the marginal distribution of X follows $NB(r,p)$, $var(X) = \frac{pr}{(1-p)^2}$ $mean(X) = \frac{pr}{1-p}$ $var(\lambda) = r \cdot \frac{p}{1-p}$ Hence, $var(\lambda) = var(X) - mean(X) = \frac{pr}{(1-p)^2} - \frac{pr}{1-p} = \frac{rp}{1-p}$ Also, we can consider X and $\lambda$ are two random variables, then $Var(X) = E(Var(X|\lambda))+ Var(E(X|\lambda))$ $Var(X) = E(\lambda)+ Var(\lambda)$ $Var(\lambda) = Var(X) - mean(\lambda)$ $E(X) = E(E(X|\lambda)) = E(\lambda)$ From the derivations, it seems like I can estimate the variance of Gamma distribution in the mixture by $Var(X)-mean(X)$. Alternatively, if we parameterize Negative Binomial X using $\mu$ and size, the $var(X) = \mu + \mu^2/size$. So $var(\lambda) = \mu^2/size$.


title: Extract Rules from Matlab Classification Tree

text: I am using Matlab's ClassificationTree.fit to generate a classification tree for my data. Let's say I type:               load fisheriris     t = ClassificationTree.fit(meas,species,'PredictorNames',{'SL' 'SW' 'PL' 'PW'})}      This gives me a ClassificationTree class variable 't'. From this class variable, how do I extract the logical rules that get me to a _particular_ class. For example, the rule to get class 'versicolor' in 'fisheriris' is:               PL >= 2.45 AND PW < 1.75 AND PL < 4.95 AND PW < 1.65.      NOTE:   1. I do not want to show any rules for other classes.   2. I do not want the output of view(t).