title: Longitudinal comparison of two distributions

text: I have the test results of a blood test administered to 2500 people four times at six-month intervals. The results primarily consist of two measures of immune response - one in the presence of certain tuberculosis antigens, one in the absence. Currently, each test evaluates to either positive or negative based on the difference between the antigen response and the nil response (with the idea being that if your immune system responds to TB antigens, you've likely been exposed to the bacterium itself at some point). In essence, the test supposes that a non-exposed individual's distributions of nil and TB responses should be basically identical, whereas a person with TB exposure will have TB responses drawn from a different distribution (of higher values). **Caveat** : the responses are very, very non-normal, and values clump at both the natural floor and the instrument-truncated ceiling. However, it's seems pretty clear in this longitudinal setting that we're getting "false positives" (no actual gold standard for latent tuberculosis, I fear) that are caused by (typically small) fluctuations in the antigen and nil responses. While this might be hard to avoid in some situations (you may only get one chance to test someone), there are many situations in which people are routinely tested for TB every year or so - in the US, this is common for healthcare workers, the military, homeless people staying at shelters, and so on. It seems a shame to ignore prior test results because the extant criteria happen to be cross-sectional. I _think_ that what I'd like to do is what I crudely conceive of as longitudinal mixture analysis. Much like the cross-sectional criteria, I'd like to be able to estimate the probability that an individual's TB and nil responses are drawn from the same distribution - but have that estimate incorporate prior test results, as well as information from the sample as a whole (e.g., can I use the sample-wide distribution of within-individual variabilities to improve my estimates of a specific individual's distribution of nil or TB?). The estimated probability would need to be able to change over time, of course, to account for the possibility of new infection. I've gotten myself totally twisted around trying to think about this in unusual ways, but I feel like this conceptualization is as good as any I'm going to come up with. If something doesn't make sense, please feel free to ask for clarification. If my understanding of the situation seems wrong, please feel free to tell me. Thank you so much for your help. **In response to Srikant:** It's a case of latent classification (TB-infected or not) using the two continuous (but non-normal and truncated) test results. Right now, that classification is done using a cutoff (in its simplified form, TB - nil > .35 -> positive). With test results presented as (nil, TB, result), the basic archetypes* are: Probable Negative: (0.06, 0.15, -) (0.24, 0.23, -) (0.09, 0.11, -) (0.16, 0.15, -)   Probable Positive: (0.05, 3.75, +) (0.05, 1.56, +) (0.06, 5.02, +) (0.08, 4.43, +)   Wobbler: (0.05, 0.29, -) (0.09, 0.68, +) (0.08, 0.31, -) (0.07, 0.28, -) The positive on the second test for the Wobbler is pretty clearly an aberration, but how would you model that? While one line of my thinking is to estimate the "true difference" between TB and nil at each time point using a repeated-measures multilevel model, it occurred to me that what I really want to know is if the person's nil response and TB response are drawn from the same distribution, or if their immune system recognizes the TB antigens and activates, producing an increased response. As for what could cause a positive test other than infection: I'm not sure. I suspect it's typically just within-person variation in results, but there's certainly a possibility of other factors. We do have questionnaires from each time point, but I haven't looked into those too much yet. *Fabricated but illustrative data


title: Among Matlab and Python, which one is good for statistical analysis?

text: I have been using Matlab for my work, but I have started learning Python lately. I employ statistical analysis, more precisely geo-statistics, in my work. I was wanting to ask, from your perspectives, which one among the two languages is good for general statistical/data analysis? What are the pros and cons, other than accessibility, for each? Feel free to edit/correct the question to make it more clear. Thanks.


title: inter vs intra group variance

text: Consider this example:               team <- rep(c("A","B","C"), times=c(7,4,10))     trip <- rep(NA,length(team))     for(i in 1:length(unique(team))){       trip[which(team==unique(team)[i])] <- 1:days[i]     }     obs  < -c(rnorm(days[1],100,30), rnorm(days[1],100,5), rnorm(days[1],100,15))     data <- data.frame(team, trip, obs)          boxplot(obs~team, data)      It is pretty clear that the variance in each team is different, but the mean is similar. How can I infer this statistically? How can I compare intra-group (within- group) variance with the inter-group (between-group) variance?


title: What is a distance or similarity metric that takes into account the improbability of a match?

text: Suppose I want to measure similarity between users. If two users match on an item that is very improbable, I want to give greater weight to that.


title: How can I weight certain factors more heavily when computing distance function?

text: I'm working on a "player phylogeny" to help make comparisons between potential NBA draft picks and current players. An example is shown below: ![enter image description here](http://i.stack.imgur.com/giv5V.png) I'm using the dist() function in R, and am assuming that all factors are equally weighted. It's clear to me by looking at the results that some factors need to carry more weight. So my question is pretty simple. How do I go about doing this?


title: How to equate hit probabilties on 2 different surfaces

text: I'm looking for help in determining how to equate the following:   1. We have a surface of 13 cm X 10 cm.    2. We have a 95% probability of hit on a 16 cm diameter surface. I would like to equate $B$, probability of hit to the 13cm X 10 cm surface. 95% probability of hit on a circle of 16 cm represent what probability of hit on the rectangle of 13cm X 10cm. What is the mathematical steps to arrive to the solution? It is a job related problem not a school one.


title: Is Hansen's skewed-$t$ distribution the same as the skewed-$t$ distribution which is a special case of GH Distribution?

text: I recently studied two asymmetric _t_ distribution both with a name of _skewed-$t$_. I am confused with their differences or are they actually the same? The first one is introduced by Hansen (1994) with pdf: $f(x;\nu,\zeta)=\begin{cases} \begin{array}{cc} bc\left(1+\frac{1}{\nu-2}\left(\frac{bx+a}{1-\zeta}\right)^{2}\right)^{-\frac{\nu+1}{2}} & \quad,\text{if }\: x<-\frac{a}{b}\\ bc\left(1+\frac{1}{\nu-2}\left(\frac{bx+a}{1+\zeta}\right)^{2}\right)^{-\frac{\nu+1}{2}} & \quad,\text{if }\: x\geq-\frac{a}{b} \end{array}\end{cases}$ where $a=4\zeta c\frac{\nu-2}{\nu-1}$, $b^{2}=1+3\zeta^{2}-a^{2}$, $c=\frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\pi(\nu-2)}\Gamma(\frac{\nu}{2})}$ and $\nu$ is DoF while $\zeta$ denotes the skewness. The other one is a special case of Generalized Hyperbolic Distributions $X \thicksim GH_d(\lambda,\chi,\psi,\mu,\Sigma,\zeta)$ when $\lambda=-0.5\nu$, $\chi=\nu$, $\psi=0$ Any suggestion will be appreciated! **Reference** Hansen, B.E. (1994), Autoregressive conditional density estimation, _Intern. Econ. Rev._ , vol. 35, no. 3, 705-730.


title: Estimation/simulation of homing with time effect

text: I have birds translocated from site A (original, capture site) to site B (new, release site) and I want to analyse homing behaviour. Translocation of individuals was performed continuously (several weeks), thus a **bird released at the beginning of the experiment had theoretically highest chance to be recaptured at the original site comparing to a bird which was released before the end of the experiment.** My data represents date of capture(=date of translocation and release) and date of recapture. Is there correct method how to estimate probability of recapture in a population (how to include such time effect)? Would it be possible to simulate homing somehow (how fast they come home)? Thank you for any guidance in advance!


title: Model selection in a paper: what to say about the dropped variables?

text: I have a question (simpler than my previous post today I hope!), which is probably very stupid as nobody has never asked it before. Lets say I am trying to explain the effect of 3 variables (A, B and C) on a dependent one (Y). Biologically speaking, A and B should really have an effect on Y. So I am testing: _Y ~ A + B + C_ But when I use a model selection method (whatever the method is), the 'best' model, the one that fits the data the best, drops A. So I end up with: _Y ~ B + C_ **What can I say about A then?** Can I cite something to justify the dropping (the F statistics, something about the AIC/BIC, etc)? **OR** If I need to show that A has no effect, do I need to use the full model anyway?


title: How to properly set up hypothesis testing in an operations reporting environment

text: I am in charge of quality control of n = 400 reports. Most of the quality checks are automated, but some have to be manual. I need to figure out how many reports my team needs to check manually in order to be fairly confident that all of the reports are error free in the areas which require manual checks. How can I correctly use hypothesis testing to answer this problem? My first instinct is to start with something like this:               Null hypothesis: at least one report contains errors.     Alternate hypothesis: all n reports are error free.      To collect data, we can check k reports manually and find numErrors errors in those k reports. Then,               If numErrors > 0:         Null hypothesis is true with 100% probability     Elif numErrors == 0:         Alternate hypothesis is true with x% probability      What I am struggling with is how to figure out the optimal value of k, and once we have k, how to determine the x% probability the Alternate Hypothesis is true if numErrors == 0. Any ideas? Thanks!


title: How to define train and test sets in financial time series for estimating machine learning parameters

text: After reading some material, I found few options for defining train and test sets:   1. Just splitting with no change.   2. Accumulating/moving window of train set.   3. Leave a relatively small (warming) period between test and train sets, and then use window again (including the warming period). What should be the most accurate way for applying machine learning algorithms and parameters estimation? Thanks in advance!


title: Misunderstanding of Monte Carlo Pi Estimation

text: I am fairly sure that I understand the how Monte Carlo integration works but I am not understanding the formulation of how it is used to estimate Pi. I am going by the procedure outlined in the 5th slide of this presentation http://homepages.inf.ed.ac.uk/imurray2/teaching/09mlss/slides.pdf I understand the preliminary steps. Pi is equal to 4 times the area of a quarter of the unit circle. And the area of the top-right quarter of the unit circle centered at (0,0) is equivalent to the integral of the curve that is the top-right quarter of the unit circle in $0<x<1$ and $0<y<1$. What I don't understand is how this integral is $\iint I((x^2+y^2)<1)P(x,y)dxdy$ where $P(x,y)$ is uniformly distributed in the unit square around the quarter circle (i.e. it is always equal to 1 if $0<x<1$ and $0<y<1$ and 0 otherwise). So this would mean that $I((x^2+y^2)<1)P(x,y)$   is the function that is the top-right quadrant of the unit circle at $0<x<1$ and $0<y<1$ but I do not understand how this is true since the indicator function can only be 1 or 0. I understand that it is probably written in this way to make Monte Carlo sampling easy (i.e. it is an expectation so just sample from $P(x,y)$ and get the average of the samples applied to $I((x^2+y^2)<1)$) but it just does not make intuitive sense to me why that integral represents the area under that curve. Could someone provide an intuitive explanation of this. Maybe show how that integral was derived in a step-by-step way? EDIT: I was able to gain a better understanding by relating the expectation to an area. I will explain it here in case it helps anyone. First start with relating Pi to the area of the top-right quadrant of the unit circle $\pi=4\times A_{tr}$ Then we place the top-right quadrant into the unit square. And under a uniform distribution over the unit square, the area of the circle quadrant is proportional to the probability of obtaining a sample from it. It follows that the following equality holds $P(x^2+y^2<1)=\frac{A_{tr}}{A_{square}}$ and $A_{square}=1$ so $P(x^2+y^2<1)=A_{tr}$ And substituting into the original equation $\pi=4\times P(x^2+y^2<1)$ and it is also true that $P(x^2+y^2<1)=E[I(x^2+y^2<1)]$ which is equal to the original double integral. So I understood it by relating the area to a probability then relating that probability to an expectation that is equivalent to the integral. Let me know if I have made any mistakes.


title: How to apply a Kalman filter to use both previous and future measurements of a random variable?

text: I'm trying to estimate the state of a Gaussian random walk with central tendency based on time series measurements with varying uncertainties. My random variable has the following form: $ \frac{d x}{d t} \equiv F(t) - \alpha x $ Where F is a Gaussian random variable. I've noticed that this problem is analogous to the velocity of a bubble experiencing Brownian motion. (See for example, F. Reif, _Fundamentals of Statistical and Thermal Physics_ , p. 565). As a result of the $ -\alpha x $ term, the position has a central tendency (i.e. the variance does not become infinity as time approaches infinity). Now, like any good physicist, I know that I cannot exactly measure the value $x$. The best I can do is to measure it at time $ t_i $ within some uncertainty, $ \sigma_i $. Using a Kalman filter, I can estimate the value of $x$ from several measurements. Let's call that $ \hat x $. The approach is as follows. For each measurement, we compute: > $ \delta t = t_i - t_{i-1} $ > > $ P(t) = P(t_{i-1}) * e^{-\alpha\, \delta t} + \langle x^2 | \delta t > \rangle $ > > $ K = {{P}\over{P + \sigma_i}} $ > > Our incoming est estimate of $x_i$: > > $ \hat x_{i-} = \hat x_{i-1} \, e^{-\alpha\,\delta t} $ > > $ \hat x_i = \hat x_{i-} + K [x_{obs} - \hat x_{i-}] $ This works great for propagating our estimates forward in time. My question is: If I have measurements at times that span the time at which I want the best estimate, how do I compute an $\hat x(t) $ where $ t_i < t < t_{i+1} $?


title: Efficient calculation of Poisson incidence rates in R?

text: Poisson derived incidence rates can be fairly nice in medical research. I noted this beautiful comparison in The Lancet: ![enter image description here](http://i.stack.imgur.com/VrmQe.png) I believe that kind of calculation, when accompanied by a Cox regression is quiet informative. How would You do something similar to that? I struggle with a similar task, but I am interested in blood pressure (4 categories) as the main exposure of interest. I have info regarding survival time and event for each individual. Is there any efficient way in R, that I could not find via Google nor R-bloggers. There are some tutorials on how to calculate incidence rate ratios, but those are not the same as incidence rates, as far as I know. Grateful for any advice that You might have.


title: Repeatedly measuring accuracy against the hold out set

text: I have an iterative document classification task, corpus size = 300,000 documents. The labels are binary valued (yes/no). I wanted to know whether the following methodology is valid. The assumption is that there is an oracle to label the documents. Randomly select 500 documents and label them using an oracle, lets call this the hold out set H. Set H is never used for training purposes. Select another 500 document set, Tr, label them using an oracle, train a SVM classifier and evaluate the accuracy of the resulting model on set H. Sample another 100 documents, get them labeled from the oracle, add them to Tr, train a SVM classifier and evaluate the accuracy on set H. Repeat the procedure until you reach the desired accuracy level and stop. Personally I feel that a hold-out set so small may not be able to account for all the variation in the corpus in order to determine whether to stop the iterative learning process. But I am not able to quantify this. Would also appreciate some pointers to research papers.


title: Maximum Likelihood Estimation Identity question

text: Is summation of logarithm of a random variable equal to the logarithm of the sum of the random variable.


title: Compute the total variation for multinomial distributions

text: Suppose I have two discrete random variables, $P$ and $Q$, with probability mass functions given by $p(x)$ and $q(x)$. I know that these random variables are multinomials, generated by choosing $K = 5$ English words each. Therefore, a sample of $P$ looks like a very long vector, in which the $i$th entry is the number of times that the $i$th word from the English language appeared among the 5 words in $P$. Represented that way, samples of $P$ and $Q$ look like super-long strings containing thousands of zeroes, with a few non-zero elements (totaling to 5). I have 126 samples each from of $P$ and $Q$, and I would like to try to answer the question _how different are $p(x)$ and $q(x)$_? **Ideally, I would like to estimate the total variation, or $L_1$ norm**. From a couple different sources, I am getting essentially the following advice: (As you read this, recall that, while the $x$ in $p(x)$ is a vector with $L_1$-norm of 5, the distribution $p(x)$, in a multinomial, and so is parametrized by $p_i$-values, probabilities of a given word turning out to be the $i$th English word). **The advice** : Look at the frequency whith which given words appear in your set of samples for $P$, and take that to be the probability of turning up that word. I.e., the probability of the $i$th English word is found by going through each of the $K=5$ words in the the $N=126$ samples, and calculating: $p_i = \frac{1}{N*K}\sum_{n=1}^N \sum_{k=1}^K[P_{n,k}=i]$. Do likewise with the samples of $Q$ to get $q_i$-values. Next, he said, compute the divergence between these empirical distributions: $L_1(p(x),q(x)) = \sum_i|p_i - q_i|$ _This can't be right!_ Imagine that there are thousands of english words (there are), and with only 126 samples, it is very likely that, _even if $P$ and $Q$ came from the same distribution_ , many words would occur only once (in $P$ but not $Q$, or vice versa). This would lead to a greatly overestimated $L_1(p(x),q(x))$. A couple quick simulations show that is true. Am I missing something? And is there a method for calculating $L_1$-norm, using only samples, between two multinomials? Related:   * enlightening, but no computation advice   * implies the dubious (to me) approach that I described   * they address a problem, but not the one I'm worried about. It seems completely unprincipled to calculate KL that way. Besides, I'd prefer a _real metric_   * Now this is more like it, but I can't get the bounds I need. **Edit: More Context** :   Here, $P$ and $Q$ represent 5 labels that people provide to describe images. I have two groups of people, $\mathcal{P}, \mathcal{Q}$ and I want to be able to tell "how different" they are. I would like the answer to be in terms of some metric over the space of probability distributions, so that I can use it to make information theoretic claims. I.e., it isn't useful to me to just be able to say that $P$ is _more_ different from $Q$ than $R$ is from $Q$ -- rather _how_ different? $L_1$ has the desireable property that it is directly bounds the distinguihsability of $P$ and $Q$.


title: Maximum Rank Correlation for panel data

text: Let $Z=(Y, X)$ be an observation from a distribution $P$ where $Y$ is a response variable and $X$ is a vector of regressors. Assuming the following model: $Y = F(X'\beta, u)$ where $X'\beta$ is a linear index, $u$ is a random disturbance ($u \perp X$) and $F$ is a strictly increasing function in each of its arguments. We know that if $Z_i$ is a sample of independent observations from $P$, then we can estimate $\beta$ by simply maximizing the Rank Correlation (Han, 1986) between $Y$ and $X'\beta$. QUESTION: would the same result apply for panel data (e.g., $Y_{ij}$ is the rank of item $j$ for the observation $i$ and there are $j=1,2,..., n$ items for each $i$)?


title: What are the relative merits of testing significance of a regression coefficient

text: If I have several hundred coefficients generated by running multiple variable regression model (keeping it as broad as possible by not specifying the nature of the predictors and outcome variable), it would appear to me that I have two options for assessing the significance of any one result.   1. Identify the largest standardized predictors by calculating z-scores for regression coefficients and assessing how many exceed a fixed threshold (e.g. 2 Z-scores from 0).   2. Use permutation testing to simulate the distribution of covariates under a null hypothesis (of no covariates being associated with the outcome) and calculate a p-value by comparing the actual results, ranked by their order of significance calculated from Z-scores, to the distributions of ranked coefficients in the permutation tests My question is: What are the advantages and disadvantages of each approach?


title: How can I display the formula/expression of a lm result in a plot in R?

text: For example, if I was to plot a trendline in Excel I can choose to include its formula on the chart, it will display something like `myY = 0.5 + 0.5 * myX`. R has quite advanced functionality for creating sophisticated mathematical expressions, so I assume this is fairly straightforward to do. I'm just not sure how to do it.


title: Multlinear regression: analysis of residual of transformed response and predictor variables

text: In the first step of modeling a regression equation I came up with the following model: $T_c = 26.73 + 0.042{\rm Sc} + 0.247{\rm Lc} - 14.709{\rm Lf} + 1.41{\rm Lu} - 0.214{\rm Fc} + 0.041{\rm Ad} - 64.308{\rm Sr} - 20.341{\rm Rc}$ Analyzing the residuals it is obviously noticed it does not follow a normal distribution; the relationship between the response variable and predictor is also found to be non-linear. After the attempt of some transformation in both predictors and response variable, and also stepwise procedures (to minimize the number of variables), I came up with the following model: $T_c^{1.5}= -9.868 + 1.845\ln({\rm Sc}) + 10.23\sqrt{{\rm Lu}} + 6.866 \times 10^{-5} {\rm Ad}^{2.5} + 7.772 \times 10^{-3} {\rm R_c}^{-4.47}$ I am using the method of least square to estimate coefficients and parametric statistic tests based on normality assumptions to test the performance of the model.   1. My question is how the residual must be calculated and analysed: I am using R language and using function `lm`, it returns the residual as residual = $Tc^{1.5}$ predicted - $T_c^{1.5}$ observed. However, I feel the residual should be: residual = original $T_c$ - the result of the right side of the above equation to the power of 2/3. Could you shed a light on this please? Am I missing some simple theory on the model that I am not being able to visualize?   2. The same question is valid to r squared: which one should I choose to calculate r squared, the result of the right side of the above equation to the power of 2/3 (as the predicted variable) or $T_c$ predicted, that is, $T_c^{1.5}$?


title: Accept/reject for beta distribution

text: I am required to show that the following accept/reject algoritmh produces observations from a beta distribution with parameters $\alpha$ and $\beta$. * * *   1. Generate $U_1$ and $U_2$ iid uniform(0,1) random variables. Set $V_1=U^{1/{\alpha}}$ and $V_2=U_2 ^{1/{\beta}}$   2. Set $W=V_1+V_2$. If $W\leq 1 $, set $X=V_1 /W $; else go to step 1.   3. Deliver $X$ * * * Thus, I need to find the conditional distribution of step 2. Now we have: $$P \left[X \leq x \right]= P \left[ V_1 /W \leq x | W \leq 1 \right]= \frac{P \left[V_1 /W \leq x \cap W \leq 1 \right]}{P \left[ W \leq 1 \right]} $$ A simple transformation shows that $V_1 \sim B \left( \alpha ,1 \right) $ and $V_2 \sim B \left(1, \beta \right)$ but other than that I do not know how to proceed as I cannot simplify the fraction. Had I known the distribution of $W=V_1+V_2$ things would have been easier but this is not the case here. Can you please help me compute these probabilities? Thank you.


title: How to calculate confidence interval for difference (or ratio) of two different time series?

text: I have two time series which are sampled at the exact same times. I would like to calculate a confidence interval either for the ratio between the two or the difference between the two. The values tend to be somewhat similar, but both have a large amount of what looks like random noise superimposed on them (much larger than the potential difference, at any particular point in time). What is the appropriate method to calculate a confidence interval in this case? P.S. The data is very different from being normal; it is roughly flat-top, pretty similar shape (but perhaps shifted) for the two series.


title: Books on non-parametric theory with correlated data

text: **Short question:** I am looking for books that deal with correlated data in a systematic and theoretical way. * * * **Long version:** Right now I am developing approximation algorithms for time series data (more concrete I want to approximate a function $f:[a,b]\rightarrow\mathbb{R}^3$ where $[a,b]$ is a time interval). The raw approximation of that data shall not be of concern for this question. However, what has a certain relevance for my works is a proof of correctness of the developed algorithm. That means that I have to take the measurement error of the data that is to be approximated into account. Problem with that is that I do not know a lot on the distribution of the error. Furthermore, I know that high absolute errors are likely to be correlated. Thus, I need to use statistical methods that do not fail completely with correlated data for my theoretical examination of the data. (Intuitively I would guess that quantiles already might be enough for my purpose, if I can justify their use reasonably.) My problem with that is that I do not have a strong knowledge on methods that hold with correlated data. All statistical methods I remember have that small but important _iid_ among the list of their requirements. Thus, the question is: Can you recommend any good books that build up statistical methods for correlated data? As I'm looking for methods and used theory, I do need a book dealing with probability theory rather than just practical methods. My background on statistics is a bit rusty but certainly should be enough for understanding theoretical books as well. For me it is important to gain an understanding on the mathematics used for that kind of problem rather than just to get to know a lot of different tests.


title: Sum of random variables decomposition

text: Let $0< \delta<\frac{1}{2}$. Why do we have for $X_{1},X_{2}$ iid: \begin{equation} \{X_{1}+X_{2}>x\}\subset\{X_{1}>(1-\delta)x\}\cup\{X_{2}>(1-\delta)x\}\cup\{X_{1}> \delta x,\ X_{2}>\delta x\}? \end{equation}


title: Hypothesis testing for a correlation that is zero or negative

text: I would like to test for a correlation that is zero or negative using the following hypothesis test:   $H_0: p>0$ (Null hypothesis: the correlation is positive)   $H_A: p\le0$ (Alternative hypothesis: the correlation is zero or negative) Since this is different from the usual correlation statistical hypothesis tests, $(H_0:p=0, H_A:p\ne0)$, I am having difficulty in figuring out the formulas. How would I perform the calculations for this hypothesis test? (Of course the calculations for $r,\ r^2,\ SS_x,\ SS_y,\ n,$ etc. are all known.) Thank you in advance for your time!


title: Are all attributes/data points inherently nominal?

text: I'm trying to design a statistics package in .NET, and I'd like to apply object-oriented design principles to some of the data structures that I'm developing. As I look through sample data sets for simple classifiers, even though there exist other types of attributes (ordinal being the only other one I've seen, with discrete and continuous being variations on that), ultimately they are all able to be classified by a distinct value in some way. This makes me think that all attribute types are inherently nominal, and that there should be a base type for nominal attributes, as the statement for all other attributes seems to hold true: > <attribute type> **is** a nominal type (Note this doesn't mean it's mutually exclusive) If that's the case, then the class that represents the type of an attribute on a fundamental level would be a nominal attribute type, with the certain methods/properties that are inherent to _all_ attribute types. Off the top of my head, the first method would be to be able to determine whether or not an instance of the variable type is the same as another instance of the same attribute type. Given the above, is every attribute type inherently nominal, even if it might possess other qualities beyond nominal qualities ( _e.g._ , ordinal attributes could have the concept of distance between two instances)?


title: Comparing a sample distribution to an implied distribution

text: I have a task that is specific to inventory management that is currently driving me crazy. To summarize the problem: We regularly must monitor inventory settings to ensure that they represent true demand, to ensure that we are not overstocking nor understocking. The standard procedure is to schedule reviews, where high volume items get reviewed monthly, and low volume items quarterly. But the vast majority of SKUs don't get adjusted; less than 9% are adjusted after each review...meaning that 91% of SKUs do not have a significant change in their demand patterns to warrant any change in their inventory parameters. I'm trying to design a filter so that instead of reviewing every SKU, I only review the SKUs where demand patterns have likely changed from their previously implied demand. To give background into the domain of inventory management: 1) Inventory SKUs have settings that manage their stock levels given a specified depletion rate and variability. There are practical parameters (eg. lead times) as well as calculated parameters (eg. avg daily demand), and often heuristic parameters (eg. targeted in stock rate). 2) Generally speaking, if I know all of the currently used parameters (or even n-1), I know the implied demand distribution that results in those parameters. For example, if I have a safety stock of 96 units, avg daily demand of 10, lead time of 5 days, and a target in-stock rate of 98%, and demand is normally distributed, then I know that the implied demand can be represented by a normal distribution with mean of 10 and standard deviation of 10. I'm trying to find a method that can help me assign probabilities that a sample of demand came from an implied demand distribution with explicit parameters. Is there a way to calculate this probability? It sounds like a conditional probability problem, but I'm not sure how to construct a calculation to determine conditional probabilities using parameterized distributions. To phrase this question as an example: If my last 10 days of demand are c(6,7,7,5,7,8,9,4,4,9), what is the probability that this sample came from a normally distributed population with mean of 5 and stdev of 3?


title: Scatter plot of interaction

text: **Is drawing a scatter plot, with fitted regression lines for subgroups, the way to graph a possible interaction in SPSS?** I may be missing something but I could not find a specific way to do it. I have only just migrated to SPSS from Excel! My IV and DV are interval scale measures and the factor is a categorical variable. This is not a homework question! Thanks.


title: Expected hitting time of a level $a$ for Brownian motion

text: Let $\{W_t,t \geq 0\}$ be a standard Brownian motion under $\mathbb P$. Let $T_a$ be the hitting time of level $a$, that is: $$T_a= \text{inf}\{t \geq 0:W_t=a\}.$$ From a proposition, we know that $$\mathbb E[\exp(-\theta T_a)] = \exp(-a \sqrt{2\theta}).$$ How do we make use of the above preposition to calculate $\mathbb E[T_a]$? I am supposed to obtain $\infty$ as the expectation.


Is a normalized version of an exponential family distribution still an exponential family distribution? Here "normalized" means making its mean zero and variance one. According to the following definition of an exponential family distribution from Wikipedia, since both variance and mean are functions of the parameter $\theta$ alone, after normalization wrt variance we still have an exponential family distribution, but after normalization wrt mean $\mu(\theta)$, $$ f_{X'}(x|\theta) = h(x+\mu(\theta))\ \exp[\ \eta(\theta) \cdot T(x+\mu(\theta))\ -\ A(\theta)\ ], $$ I am not sure. > A single-parameter exponential family is a set of probability distributions > whose probability density function (or probability mass function, for the > case of a discrete distribution) can be expressed in the form $$ > f_X(x|\theta) = h(x)\ \exp[\ \eta(\theta) \cdot T(x)\ -\ A(\theta)\ ] $$ > where $T(x), h(x), η(θ)$, and $A(θ)$ are known functions. > > An alternative, equivalent form often given is $$ f_X(x|\theta) = h(x)\ > g(\theta) \exp[\ \eta(\theta) \cdot T(x)\ ]\, $$ or equivalently $$ > f_X(x|\theta) = \exp[\ \eta(\theta) \cdot T(x)\ -\ A(\theta) + B(x)\ ] $$ > The value $θ$ is called the parameter of the family. Thanks!


title: Binary or multinomial logit?

text: I am having some doubts on whether I should use binary or multinomial logit. I am modelling product choice. The dependent variable is the choice of either shampoo or toothpaste. Although my dependent variable has two outcomes, I am wondering whether it is appropriate to use multinomial logit instead of binary. The case is that the customer does not choose among these two products, but among the whole assortment. And I have taken a sample for only 2 products. Do you think I should use the multinomial logit?


title: Territories from observations

text: I have a number of animal observations, and want to deduce the number of territories (i.e. the number of individual animals) from this. More formally, the problem can be stated as follows: Each observation is a triple consisting of two coordinates and a date. Observations on the same date are assumed to be different animals, whereas observations in close proximity but on different dates might come from the same animal. I'm looking for the **minimal number of territories** which can explain all of these observations. Each territory is assumed to be a circle, whose radius is known in advance and the same for all territories. An observation can be assigned to a territory if it falls within that circle, but two observations made on the same date must not be assigned to the same territory. In this model, territories may overlap. An observation may lie within one circle and still be assigned to a different one inside which it lies as well. **Is there an established procedure or algorithm to answer this kind of question?** If so, is it an exact optimal solution or an approximation? Do you know of any R package to perform this kind of computation? At the end of the day, it would also be nice to know which observations have been assigned to the same territory. I believe that with the above definition, this assignment is not uniquely defined in most cases. I'd be also interested in established procedures to make a decision among these alternatives. If you know of techniques to address this general kind of question, i.e. identifying territories, but the assumptions are different from those stated above, I'd still value a comment providing a useful name or reference.


title: ANOVA data prep - sum across

text: I have an experiment with two treatments (3 & 4 factors each) and 4 true replicates on 48 plots. Those experimental plots recruited various quantities of 25 different species so the dataset has 1200 data points. My question is this: if I run a GLM ANOVA on the two expt tmts and their interaction, but not include species in the analysis (I'm interested in totals), do I first need to sum the plots across species? This takes the 1200 point dataset down to 48 data points of course. I was never taught to do this, but that is what someone else did.


title: How to handle missing data in a logistic regression using SPSS?

text: I have a data-set of genetic variants which I'm trying to use as predictors for a simple phenotype, and for starters I use a binary logistic regression in SPSS. I have around 900 individuals, and for each individual around 50 variations and a phenotype. However, I get an unreasonably high amount of removed values when I run the analysis (there is some missing variants data all over the entire table), i.e., only around 50% of my measurements are actually used, and I can't find the exclusion-cut-off SPSS uses for this anywhere. Does anyone know the cut-off that SPSS employs in this? Or does it remove a measurement once it finds a single missing variant?


title: Boruta score goes to minus infinity

text: I'm running the Boruta algorithm with a $179\times 36$ predictor matrix and a numerical response. Most of the variables have a score going to `-Inf`. Should I conclude there are not important or there's a problem ?               bor <- Boruta(X,y)     > dim(X)     [1] 179  36     > bor$ImpHistory[,1]       [1] -0.33842437  0.03724477  0.77207194 -0.15966460  0.73550293 -2.11845770 -0.28029409 -1.73339443  1.80120877 -1.22191877 -0.45272011 -0.71728205      [13] -0.57268506 -0.48498812  0.61710237  0.58776416 -0.38667490  0.51795751  0.32422847 -1.79565022  0.26975995 -1.10919356  0.77731727 -1.10266785      [25]  0.19157178 -1.16694014  0.53795858  0.96632680 -0.99366419 -1.22646545 -1.71817624 -0.68100802 -1.41193813 -0.18065431        -Inf        -Inf      [37]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf      [49]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf      [61]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf      [73]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf      [85]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf      [97]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf     [109]        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf        -Inf     [121]        -Inf


title: Using ANOVA and t-tests with pre-aggregated data

text: I have a set of experiments that I'd like to run some significance tests over. It looks like:               Word Feature Model Test Score Correct Total     allow.v woc hac zellig 0.382353 26 68     allow.v woc hac bellig 1.000000 0 0     allow.v woc kmeans zellig 0.382353 26 68     allow.v woc kmeans bellig 1.000000 0 0     run.v woc eigen zellig 0.308824 21 68     run.v woc eigen bellig 1.000000 0 0     run.v woc agglo zellig 0.323529 22 68      This combines two different types of tests, which i'm calling zellig and bellig for a lack of better names right now. Using this data, I considered treating each score as a response variable that is dependent on the Feature and Model, with each as an observation. With this format, it's pretty easy to do standard evaluations like:               zelligData <- subset(senseData, Test == "zellig")     anova(lm(Score ~ Feature * Model, zelligData))     pairwise.t.test(zelligData$Score, zelligData$Feature, p.adj="none")     pairwise.t.test(zelligData$Score, zelligData$Model, p.adj="none")      That's well and good, however, Score is actually an aggregate value of many observed events, it's the average number of times each combination resulted a correct evaluation. The Correct and Total fields count the two numbers used to generate Score. And now, I'd like to do significance tests for each word with each event counted by Total treated as an observation. Concretely, I'd like to use the above commands on data that looks like:               Word Feature Model Test Score     allow.v woc hac zellig 1     allow.v woc hac zellig 1     allow.v woc hac zellig 1     allow.v woc hac zellig 1     ....     allow.v woc hac zellig 0     allow.v woc hac zellig 0     allow.v woc hac zellig 0     allow.v.pos hac zellig 1     allow.v.pos hac zellig 1     allow.v.pos hac zellig 1     ...     allow.v.pos hac zellig 0     allow.v.pos hac zellig 0     allow.v.pos hac zellig 0     ...      However, I'd like to avoid creating this secondary view of the data, as some words have a very large number of events and it seems wasteful to turn the aggregate values into the expanded form. Is it at all possible to perform anova and paired t-test using my already aggregated counts of events? Would this evaluation be different from just using the aggregated scores as they are? I imagine they would as the aggregated Score value doesn't indicate how many times events there really were in total. Thanks in advance! Edit: After reading how to do ANOVA manually, I should rephrase this question. Is there anyway to pass in the matrix of summary statistics between each factor, as described in Section 2 (Summary Statistics) of the instructions for doing it manually, to a function in R and have it compute the ANOVA result?


title: In PC SAS, how do you copy & paste from the output window?

text: Sometimes, I just want to do a copy & paste from the output window in SAS. I can highlight text with a mouse-drag, but only SOMETIMES does that get copied to the clipboard. It doesn't matter if I use "CTRL-C" or right click -> copy, or edit -> copy Any other SAS users experience this, and do you know a workaround/option/technique that can fix it? Sometimes, I can fix it by clicking in another window, and coming back to the output window, but sometimes I just have to save the output window as a .lst and get the text from another editor.


title: When can you use data-based criteria to specify a regression model?

text: I've heard that when many regression model specifications (say, in OLS) are considered as possibilities for a dataset, this causes multiple comparison problems and the p-values and confidence intervals are no longer reliable. One extreme example of this is stepwise regression. When can I use the data itself to help specify the model, and when is this not a valid approach? Do you always need to have a subject-matter-based theory to form the model?


title: Time series modeling with dynamic regressors in SAS vs. in R

text: I am using both R and SAS for the time series modeling. There is an option in SAS that I could not find so far in any packages developed in R for the time series modeling such as TSA or forecast package, at least to the best of my knowledge! To explain more, if we use the windowing environment in SAS to fit an ARIMA model with a regressor, we basically choose: Solution->Analysis->Time series Forecasting System->Develop Models   Then Fit ARIMA model -> Predictors->Dynamic Regressors If we ask to forecast this model, SAS says “The following regressor(s) do not have any forecasting models. The system will automatically select forecasting models for these regressors”. This means that we have not provided the values of the regressors over the forecasting period, and the system tries to find a model for that. My questions:   1. Is there any package in R with the same capability (explained above) as in SAS to forecast an ARIMA model?    2. How can SAS automatically forecast the regressor(s) and based on what models?


title: Compare the "normality" of two scales

text: I have two scales measuring the same construct, and I would like to know which scale is distributed more closely to a normal distribution. There are various possibilities for comparing a distribution with a normal distribution (e.g., kolmogorov-smirnov test, etc.), but is there any possibility to compare the deviance from normality between two scales directly?


title: NBA Points Per Game

text: I have a database with many thousands of entries of box scores. Each box score has a player name, date played, points scored, minutes played, and some other stats like rebounds, assists, etc. I want to use this information to find out how the frequency of points per game varies by each half a point difference. I'm not sure I'm stating the problem exactly, so here is an example. So let's say we're looking at Kevin Durant, his average and median points per game for this season is 31. So let's say the next game he plays, he's gonna get above 31 points 50% of the time, and below 31 points 50% of the time (and disregard all cases where he hits exactly 31). How often will he get above 30.5 points? How often will he get above 31.5 points? Or rather, how do I use my database to figure this out? If I just check historical games with Kevin Durant and see how often he went over 31.5 as opposed to 31 as opposed to 29.5 the sample size won't be large enough. And if I select 5,000 random box scores with the criteria "select box score where points scored was between 29 and 33", I believe I will seriously bias the distribution and won't be able to learn from it. Thoughts?


title: Analyzing online sales where data are only produced when the sales is made

text: I am trying to model data on the number of online sales are made within a fixed sale period of 3 days. Data are generated only when the sale is made. I think for this kind of data I will be using a type of time series model for count data with inflated zero samples. But this dataset adds one more layer of complexity of the time limit. I don’t know any time series statistical model that limits the length of time. I may consider the stochastic process as a (Poisson like) draw of the limited number of trials: if it’s hourly then 72, if it’s by minutes then 4320 trials. By then, if it’s by minutes, the data ends up with too many zero observations, which ruins the effectiveness of the model. Do you have any better idea? * * * The data is structured like this:                2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-02 19:15:39        45       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-02 22:32:06        46       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-02 23:39:05        47       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-02 23:20:09        48       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 01:32:09        49       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 03:33:11        50       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 02:47:07        51       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 04:05:05        52       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 04:01:08        53       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 06:53:09        55       2011-10-31 07:01:55 2011-11-03 06:59:59  2011-11-03 07:00:54        56       The first and second times (1st and 2nd column) are the start and end times of the deal, which is about 3-day apart as the deal lasts for 3 days, and of course fixed for each deal. The third time (in the 3rd column) is increasing as the deal goes by with the cummulative number of items sold on the last columne as a _positive integer_ in the 4th columne. The data is generated only when the sale is made, so the time between succesive sales is _irregular_ , which could be a key term to search for the appropirate model. After further search for the appropriate model, I found integer autoregressive model (INAR) may be used, but I am not sure... Since the sale is recorded by the second, I may think the data series with T=259200 observations having many (too many!) zero observations. I may transform the data by minutes by summing over the sales numbers withine one minute to reduce to T=4320. This will reduce the ratio of zero-total observations smaller, and apply INAR or some kind of zero- inflated stochastic model (truncated poisson or something like that). Could you give me any good ideas?


title: Dubious use of signal processing principles to identify a trend

text: I am proposing to try and find a trend in some very noisy long term data. The data is basically weekly measurements of something which moved about 5mm over a period of about 8 months. The data is to 1mm accuracey and is very noisy regularly changing +/-1 or 2mm in a week. We only have the data to the nearest mm. We plan to use some basic signal processing with a fast fourier transform to separate out the noise from the raw data. The basic assumption is if we mirror our data set and add it to the end of our existing data set we can create a full wavelength of the data and therefore our data will show up in a fast fourier transform and we can hopefully then separate it out. Given that this sounds a little dubious to me, is this a method worth purusing or is the method of mirroring and appending our data set somehow fundamentally flawed? We are looking at other approaches such as using a low pass filter as well.


title: Can I include nominal predictor variables in ordinal regression?

text: I am in the process of performing ordinal regression on a large data set that measures different types of variables (nominal, ordinal, continuous). My outcome variable is ordinal, hence the ordinal regression. My question is; can I include nominal data (e.g. a binary coded variable representing a 'yes/no' response) as a predictor variable?


title: Gaussian Process goodness of fit

text: Let's say I got a Gaussian Process model $M$ based on some training data. Now I get a stream of sample data of a certain batch size coming in. The GP does not model a time series, but it's trying to regress the value at certain locations $x$, that will be visited mutliple times. I know that at some point there will be an abrupt change in the distribution the data batches are generated from. (At least on certain locations $x$) Now I'm looking for a statistically sound way of detecting this change, that is I want to find the point, when the GP doesn't model the data well enough anymore. I thought I would base this on detecting "big" changes in the likelihood function $L(\theta | x)$. However I'm not sure how to interpret "big" here, as the values of the likelihood function only have a meaning in comparison, but on their own. Note: I'm asking this, because some people said my original question was too abstract. I didn't want to change the whole question, though, as there already were some answers.


title: Can I used a repeated measures ANOVA or do I have to use a multivariate ANOVA?

text: I don't feel very comfortable with stats (still using simple programmes like SPSS rather than R!), and have a question about which is the best test to use with the following design: I have performed an experiment with n = 4 per group (three groups: drug A+ drug B treated, drug A+ vehicle treated, vehicle only). I've then extracted chromatin (= nuclear protein-DNA complexes) from these samples (12 in total), used antibodies to isolate bits of DNA and checked whether particular genes are more highly represented (= enriched) over others. What I want to test primarily is: 1) Is there a difference in gene enrichment between vehicle only and drugA+vehicle treated groups, and if so, which genes? 2) Is there a difference in gene enrichment between drugA+drugB treated and drugA + vehicle treated groups, and if so, which genes? I have tested 5 genes for each sample, and to complicated matters further I have tested two antibody batches (they are meant to contain the same thing, but as the lot differs they might not, and indeed they don't quite give me the same results). Initially, I was going to run a simple multivariate ANOVA with treatment group as the fixed factor, where I just checked for every paramater separately gene 1,antibody batch 1 gene 1, antibody batch 2 gene 2, antibody batch 1 gene 2, antibody batch 2 etc However, I wonder whether I lose power like this, because obviously I expect the result for the two antibody batches to be somewhat related. I also to some extent would expect the effect on the various genes to be somewhat related within a sample, as these genes all belong to a particular biological network. However, I suppose it is wrong to use a repeated measures test (e.g. with treatment as the fixed between-subjects factor, and gene and antibody as within-subjects factors?) Would you recommend some kind of mixed model instead? I hope this is clear enough...it is rather more complicated to explain than I thought. Thanks!


title: Logistic regression across groups

text: When running a logit/probit model for particular sets of outcomes on a set of participants (whether they did a certain behavior at least two times), how can one best control for differences when the individual observations began? In particular, one person may have started to participate 3 years ago, whereas another person may have joined just last year. I'm using STATA for the analysis and do know when the individuals began observations. Related to this question, what is the most effective way to refine the functional form of a the model?


title: Estimating the difference between two functions

text: For each function (f or g) I have some sampled X data (of about 1000 dimensions) and the associated (noisy) y value. From these two data sets I would like to estimate f-g, whose magnitude I expect to be small relative to f and g. The X data for both datasets are drawn from the same distribution. I have tried estimating f and g separately and merely taking the difference between the estimates, but it seems that is leaving some power on the table. After all, if for each data point I had samples from both f and g, it seems reasonable that the best way would be to subtract the respective y's and regress the differences, but in general I don't. Anyone have any advice, references, or just recommendations on what to google for?


title: How does Cornish-Fisher VaR (aka modified VaR) scale with time?

text: I have already posted this question in the quant section, maybe the statistics community is more familiar with the topic: I am thinking about the time-scaling of **Cornish-Fisher VaR** (see e.g.page 130 here for the formula). It involves the **skewness** and the **excess-kurtosis** of returns. The formula is clear and well studied (and criticized) in various papers for a single time period (e.g. daily returns and a VaR with one-day holding period). Does anybody know a reference on how to scale it with time? I would be looking for something like the square-root-of-time rule (e.g. daily returns and a VaR with $d$ days holding period). But scaling skewness, kurtosis and volatility separately and plugging them back does not feel good. Any ideas?


title: Bootstrap PCA with varimax rotation

text: I have microfossil data from across the Cretaceous-Tertiary boundary, 20 species (=variables) and 27 samples (=observations). I have run a PCA with varimax rotation on these data using the rgr package. I now want to determine the significance of the PC loadings using bootstrap BCa confidence intervals from package boot, but I can't get it to work. I have found using Krzanowski cross-validation that 5 PCs appear to be significant (NRETAIN in the script below is thus equal to 5). The script I have been attempting looks like this (the data matrix is called “x”):               BS.vmax <- function(d, f)     {    n <- nrow(d)          unrot <-gx.mva(d)          unrot$eigenvectors              vmax <- gx.rotate(unrot, nrot=NRETAIN)              vload <- vmax$rload[,1:NRETAIN]     }          BS.vmax.boot <- boot(x, BS.vmax, R = 99)     boot.ci(BS.vmax.boot, type = "bca")      This script returns a `BS.vmax.boot` value that looks like this:               ORDINARY NONPARAMETRIC BOOTSTRAP     Call:     boot(data = x, statistic = BS.vmax, R = 99)     Bootstrap Statistics :               original  bias    std. error     t1*   -0.745596283       0           0     t2*   -0.902044503       0           0     ...     lines removed here     ...     t99*  -0.144587913       0           0     t100*  0.076592734       0           0      which doesn't make any sense (because all biases and SEs are zero)! I will later on increase the number of bootstrap replicates to 2000; this is just for testing.


title: PCA loadings and valid principal components

text: I am a beginner and in the process of understanding PCA. I came across the following question > Suppose we have a data set where each data point represents a single > student's scores on a math test, a physics test, a reading comprehension > test, and a vocabulary test. > > We find the first two principal components, which capture 90% of the > variability in the data, and interpret their loadings. We conclude that the > first principal component represents overall academic ability, and the > second represents a contrast between quantitative ability and verbal > ability. and it states that (0.5, 0.5, 0.5, 0.5), PC1, and (0.5, 0.5, -0.5, -0.5), PC2, loadings are valid loading with the following explanation. > the first component is proportional to average score, and the second > component measures the difference between the first pair of scores and the > second pair of scores. I am not able to understand what does this explanation mean. Please help. I am sorry if this is not the right forum for this.


title: Forecasting time series regression in R using lm

text: I would like to predict using a linear model in R. The model that I have is of the form `lm(y~ lag(x))`. It seems like I should be able to predict using the predict function, but this does not look ahead into the future. Here is the lag function that I have been using from within R.               lag1 = function (x) c(NA, x[1:(length(x)-1)])      This lag function adds an NA at the beginning of the data, shifts everything down, and cuts off the last observation. I am mostly interested in the predictions made on the last line of the data. This line of data seems to be ignored by the fitted and predict functions. The best workaround I have found was to multiply each of the model coefficients by each of my dependent variables, and add them up for each prediction. None of the dependent variables are lagged at this point, when I set up this equation. The problem with this is that each time I change the model formula, the names of the variables have to be changed to match the new formula. It seems like there should be a more intuitive solution to this.


title: Statistical Test for Change in Share

text: So I've been dealing with a problem concerning the change in market share. I have the total volume generated (dollars for instance) and number of transactions that those volumes are generated from. As in below:                         Total Volume    Volume of X  Share of X(vol)   N_Total   N_X  Share of X(in N)     Period I       500            50            10%            250     20         8%     Period II      1000           80            8%             450     35        7,8%      And I want to investigate if the change in share of X(vol) from 10% to 8% is significant or not. Normally if I would test the significance of share transactions of X (from 8% to 7,8%), I could use the usual change in proportions formula. But since Volume is not binomial distributed I don't think it would be appropriate in this problem. How should this case be formulated?


title: Multiple Regression - Extreme F-statistic and R-square on 1 - too good to be true?

text: I have used a psychometric survey of 10 items which measure risk perception. I have taken the mean scores of these 10 items from a sample of 355 respondents who rated six hazards via the 10 Psychometric questions. I have performed a multiple regression where my criterion variable is a weighted composite scale from an estimation of risk consequence survey. My results are as follows - "In combination, the neutral Psychometric Dread risk factor accounted for a significant 100% of variability in neutral estimations of the severity of risk consequence, R2 = 1.0, adjusted R2 = 1.0, F (10, 355) = 18521487.85, p = .001" I have never seen an F-statistic this high, and the r-square = 1. My question is - Is this result possible, and if so, how likely it it? Thank you! Mel


title: How do I remove emotional bias in rating survey questions?

text: I want to do a survey about the quality of a slot game. At the end of a game session (when player cashes out), there will be a question on screen: "Please rate the game on a scale of 1 to 5" Player can touch to rate. I can foresee 2 bias sources:   1. Whether player won or lost   2. Whether player encountered unusual events such as hitting a jackpot I assume those 2 are the only major sources of bias; there may be other sources, of course. So in addition to the rating itself, I'll also collect these data:   1. The amount of money before and after the game session   2. The number of jackpots hit **How do I construct a model that will best remove the emotional bias in rating?** Currently my approach is using multipliers. If player wins, multiply the rating score by 0.9 (to offset the positive emotion), if player loses, multiply by 1.1 (to offset the negative emotion), if player hits jackpots, multiply by 0.7. (E.g. if rating is 3/5, player wins, no jackpot hit --> adjusted rating 3*0.9 = 2.7). But a big issue with this approach is determining the multipliers. Thank you very much!


title: Expected value of inverse?

text: If I have a random variable $V$ that is normally distributed with some $\mu$ and $\sigma$, then what is the expected value of $1/V$? I tried doing by delta method, and I get expected value $1/\mu$, but I also read that this is not the correct result. To clarify: My random variable $V$ is a variable that can be measured. The data given back to us are mean and standard deviation of the measurements. (I understand that I should be using the t-distribution to best model it, as we only have an average of 4 reads, but I am trying to work it out for the Normal case first.) The value $V$ is indicative of protein binding to ligand. I wish to understand fractional occupancy ($\alpha = \frac{current signal}{saturation signal}$), in which case I need to know the value of $V$ at saturation, which also can be estimated from the data (take the max value of all values present, this will be a consistent under-estimation). Therefore, if $V \sim N(\mu, \sigma^2)$, I expect $\alpha \sim N(\frac{\mu}{max}, \frac{\sigma^2}{max^2})$. I hope I did the math right. From the Hill equation, $$ \alpha = \frac{\mbox{ligand}^n}{\mbox{ligand}^n + K_d} , $$ therefore $$ K_d = \frac{\mbox{ligand}^n}{\alpha} - \mbox{ligand}^n. $$ Therefore, if $V$ were a random variable, then alpha is a random variable, therefore Kd should be a random variable. I'm trying to find how to model $K_d$ here. My thought was to use the delta method to get an approximation, and that is how I led myself to this question I've posted here.


title: Measuring the standard deviation in Pattern recognition

text: Sometimes in pattern recognition say Character recognition, Hamming distance is used although there are other distance measures. But if the pattern is represented in (1,0,-1) then Hamming distance is easier (Reference : Effect of Hamming distance of patterns on storage capacity of Hopfield Network). In my application, at every time step a character is presented to the network and each character is represented by 20 bits. I am designing a Hopfield network so the network is 20 input and 20 output for auto-associative memory purpose. My concern is along with Hamming Distance I want to use the standard deviation of each input for comparison. I want to record the count of hamming distances between the input pattern and the model output which is less than SD of the input pattern. Should I consider the SD of the entire training data set or 1 character at a time or one varaible at a time in which case there will be 30 SD for each input example. How do I use the standard deviation in this case? Please help Example : `Time Instant 1 : Character A is represented by $[-1 1 1 0 0 \ldots ]`$                       Time Instant 2:           B is represented by $[1 1 -1 1 \ldots]$      and so on. The data set forms a multivariate time series. The model outputs at Time Instant 1 `Ahat`, at time instant 2 `Bhat` and so on. In the general case of 1D time series it is intuitive to use SD but I am not sure how to use SD in case of patterns. Will appreciate ideas and insights. Thank you


title: Likelihood vs. Probability

text: I have difficulties with _Likelihoods_. I do understand Bayes' Theorem $$p(A|B, \mathcal{H}) = \frac{p(B|A, \mathcal{H}) p(A|\mathcal{H})}{p(B|\mathcal{H})}$$ which can be directly deduced from applying $p(A,B) = p(B) \cdot p(A|B) = p (A) p(B|A) = p(B,A)$. Thus in my interpretation, the $p(\cdot)$ functions in Bayes Theorem are somehow all probabibilities, either marginal or conditional. So I have actually thought that Likelihood as a concept was more of a frequentist view of the inverse probability. However, I have now repeatedly seen statements in _Bayesianists'_ books that say that the likelihood is not a probability distribution. Reading MacKay's book yesterday, I stumbled over the following statement > "[...] it is important to note that the terms likelihood and probability are > not synonyms. The quantity $P(n_b|u,N)$ is afunction of both $n_B$ and $u$. > For fixed $u$, $P(n_b|u,N)$ defines a probability over $n_B$, for fixed > $n_B$, $P(n_B|u,N)$ defines the likeihood of $u$."   * I understand this as follows: $p(A|B)$ is a probability of $A$ under given $B$, thus a function $\text{probability} : \mathcal{A}\to [0,1]$. But considering a given value $a \in A$ and evaluating $p(A=a|B)$'s dependency on different $b\in\mathcal{B}$'s we are actually using a different function $L : \mathcal{B}\to[0,1]$.   * Is this interpretation correct?   * Can one then say that maximum likelihood methods could be motivated by the Bayesian theorem, where the prior is chosen to be constant?


title: How do I conduct model selection for logistic regression in SPSS?

text: I am analyzing a set of clinical data where I try to predict an outcome by using certain covariates. I have already done univariate analysis and now am progressing to binary logistic regression, incorporating the covariates that have a p < 0.1 in univariate tests to the model. In doing binary logisitic regression, which method is better enter or one of the forward or backward elimination methods? If I use the enter method, should I manually include and exclude different covariates until all covariates with a final significant contribution to outcome are included in the model? Which ever the method used, what are the factors/statistical measures to consider in selecting the best method?


title: Generate two correlated ARMA(1,1) processes

text: How can I generate two correlated ARMA(1,1) data series where $d_{1,t}=\mu+\phi_{1}d_{1,(t-1)}-\theta_1(e_1(t-1)+e_1(t))$ $d_{2,t}=\mu+\phi_{2}d_{2,(t-1)}-\theta_2(e_2(t-1)+e_2(t))$ and $\rho_{12}$ is desired correlation between $e_1$ and $e_2$?


title: Can p values be used to show impact of treatment

text: Does it make sense to use the differences in p value to show a tendency or the 'importance' of the effect of a treatment. for example, I have treated a contaminated soil and I test the treatments against one control. I can see the high molecular weight contaminants are more affected by the treatment and that the **p values** (testing treatment against control) for these contaminants **are lower** than the ones of the low molecular weight ones, which if im right would mean that the significance is higher. This seems to correlate with the graph, but could this difference in significance seen through the p values be used in itself to show the difference in the impact of the treatment, I mean in the discussion of your results. Thank you in advance.


title: Why don't the results of testing $H_0 : \beta = 0$ and $H_0 : {\rm cor}(X,Y)=0$ agree?

text: I have 4 IVs in my model that directly effect the DV. The results of the correlation & regression analyses showed that: **IV1 &DV:** Pearson Correlation Coefficient: _insignificant_ Regression Beta&t values: _significant_ **IV2 &DV:** Pearson Correlation Coefficient: _significant_ Regression Beta&t values: _insignificant_ **IV3 &DV:** Pearson Correlation Coefficient: _significant_ Regression Beta&t values: _significant_ **IV4 &DV:** Pearson Correlation Coefficient: _significant_ Regression Beta&t values: _significant_ No multicollinearity problem was detected and also all the regression assumptions have been perfectly met! Why the results of the regression analysis turned out to be completely opposite of correlation analysis for **IV1** and **IV2**?!! Why they're being contradictory?! Does it make any sense? Is it acceptable? TQ.


title: Balanced repeated measures design

text: Okay this is my problem:   * I have 12 participants.   * Each participant spent 3 nights in my lab doing a reaction time task at four timepoints during the night(12, 1, 2 and 3 oclock), with one week between each of those nights.   * Each night, each participants was exposed to one of three experimental conditions, such that each participant completed each condition at the end of the three weeks, but with the order of conditions fully balanced over subjects. So in the end I have for each participant, for each condition, four mean reaction times. I wanted to use a linear mixed model with condition and timepoint during the night, and the interaction of condition*timepoint as fixed effects(since I expect reaction time to decrease less in one of the three conditions during the night). As a random effect I wanted to include subject, but I'm not confident of how I would be doing this in SPSS and whether a mixed model is the right way to go. Can someone give me a few hints on how to proceed? I'd really appreciate it! Thanks in advance


title: How to check goodness of fit for a graded response model in R

text: I am using the `grm` function in the `ltm` package in R. There is no function to check the goodness of fit of the output. How can I check if the graded response model is a good fit to the data?


title: Resampling from Monte Carlo simulation runs

text: Say I run a Monte Carlo simulation to generate a regression model of a simulation space. All variable parameters in the original sample have been selected at random from a uniform distribution. Is it then possible to resample from this sample while constraining one or more of the variables to a new distribution (say a normal distribution)? It should be assumed that the new distribution falls within the upper and lower range of the uniform distribution. Alternatively, can this be done analytically based on the regression model? The purpose of this is to continuously update an estimate of a system's performance as influential variables are specified (usually as distributions). The simulation is very time-consuming so if it is possible to resample rather than re-running the Monte Carlo simulation then this would be preferable.


title: Friedman's test for binary data - possible or not?

text: If my dependent variable is a technique scored as 1 for correct step and zero for incorrect can I use Friedmans test to compare technique scores at 3 time points?


title: What is the simplest statistical analysis to compare 2 groups on a simple survey?

text: I designed a simple survey about cancer knowledge among adults. There will be two groups: Those that have first-hand experience with cancer (you or someone in your immediate family with cancer) and those without (negative on the above). There are seven questions related to cancer knowledge, each with one right answer. My question is: What is the simplest and most easily understood statistical analysis for this type of data?


title: Reversing PCA back to original variables

text: I have a set of data that has $n$ samples described by $m$ variables. I do a PCA to reduce it to just 2 dimensions so I can make a nice 2D plot of the data. I understand that the $x,y$ coordinates (i.e., the PCA scores) for the plot are calculated by basically summing the products of the original data (after centering) by the loadings for each variable, so: $PC_1 = X_1L_1 + X_2L_2 + ... + X_nL_n$ My question is, if I pick an arbitrary point in the PCA space (i.e. a value for $PC_1$ and $PC_2$, or $x$ and $y$ in my plot), is there a convenient way to translate that back to a set of the original values (i.e., $X_1,X_2,\dots,X_n$)? Note 100% reversal is obviously not expected (since I'm only using 2 PCs), so a decent approximation is fine.


title: Quantitative results of clustering analysis

text: Currently, I am doing a clustering analysis for two sets of data. One smaller dataset (about 100 data) got ground truth labels, and one larger dataset (about 2000 data) has no ground truth labels. For the smaller dataset, obviously, I can obtain quantitative results like accuracy, sensitivity and specificity. However, for the larger dataset, I have no ground truth and couldn't get any useful quantitative results.   1. The only thing I found useful is the 'mean silhouette value', which can measure the cluster performance. However, it based on some distance measure that can only tell people how separate are the clusters. I am wondering if there are other 'better' or 'more appropriate' quantitative analysis for data without labels.   2. Because the data are without labels, I am also wondering if we can somehow have a 'uncertainty' measure about the clustering results like how confident about the cluster results?   3. For the smaller dataset with labels, except accuracy, sensitivity and specificity, any other quantitative results I can get? For the classification algorithm, we can do a cross-validation, is there any method we can use to do such a cross-validation for clustering? Also, can we get ROC analysis for clustering task? Thanks! A.


title: Which statistical test should be used to test for enrichment of gene lists?

text: I have performed an experiment to test the cellular sensitivity to a certain DNA damage agent. We have found 270 genes that were specifically sensitive to the drug and the total number of genes analyzed was 3668. 38 out of the 270 sensitive genes are classified as "DNA repair genes". If the number of "DNA repair genes" contained in the genome is 112 and the total number of genes in the genome is 3668, are the sensitive genes enrichment in DNA repair genes? Which statistical test should be used? I would appreciate if you could also tell me some tool to calculate the p-value online.


title: Can someone help me understand what type of problem I am looking at? Not sure if this classifies as hypothesis-testing

text: Please pardon me if this question is not clear. I am not sure if I am using the right terminologies. I have conducted an experiment in different environments multiple times. So my data looks something like this:               Environment1  1.2  2.1  1.1  1.5  1.6     Environment2  4.2  2.6  3.5  2.5  2.9     Environment3  7.2  4.6  5.3  4.5  1.6     Environment4  0.0  0.0  1.2  15.0 0.0     Environment5  3.2  2.4  7.2  5.5  6.6     Environment6  23.2  32.1  18.1  1.5  19.6      I can clearly see (or maybe my intuition says) that the experiment was not conducted properly in Environment4 (too low and fluctuating a lot) and Environment5 (way too high) but am not sure how to prove this. Am I supposed to rely on hypothesis-testing with the hypothesis: > The experiment was not conducted properly in Environments 4 and 6. and then use some procedure to prove this? Or is there a standard way of showing this? Can someone please help me how to approach this kind of problems? I am using R.


title: Any online software to calculate pvalue of Fisher exact test?

text: I would like to do a one-sided Fisher´s exact test for an analysis. I have not any statistic software to obtain the pvalues (no SAS, no SPSS). The 2x2 tables are of this type: Do you know any online statistical software to calculate the pvalues? I have tried with some of them but the results indicate pvalue<0.0001 but I need to know the exact number. The 2x2 tables are of this type: ![enter image description here](http://i.stack.imgur.com/BUTLq.png) Thanks a lot in advanced!


title: Bayes Theorem for Continuous Value Attributes

text: I need any solved example/data set which explain how to apply Bayes Theorem for continuous value attributes. I read the book (Machine Learning by Tom Mitchell) and found this equation. But I need some data to verify that I've understood this concept. ![Equation 1](http://i.stack.imgur.com/pmjKA.png) ![Equation 2](http://i.stack.imgur.com/LU8uT.png)


title: Unsupervised Random Forest using Weka

text: I am having some issues understanding how unsupervised Random Forest works according to Breiman. I only have unlabeled data, so the thought arose to use unsupervised Random Forest and use the resulting dissimilarity matrix as input for a cluster algorithm. One "constraint" is that I have to use Weka. Can someone explain to me how exactly I have to understand the usage of original and synthetic data? As I understand it right now I have to:   * artificially add a class attribute with a value of, e.g., "0" for the original data    * make a randomized copy of the original data resulting in the synthetic data and add a class attribute with a value of, e.g., "1" **Questions:**   1. Did I understand it correctly?    2. What am I supposed to do next?    3. Do I have to join these two data sets? Isn't that highly impractical for a, say, very large dataset, as I have to store the whole data while running it with Weka?   4. How do I get the resulting dissimilarity matrix in Weka?  I'm a little confused right now...


title: Poisson distribution help

text: In many cities, neighborhood Crime Watch groups are formed in an attempt to reduce the amount of criminal activity. Suppose one neighborhood that has experienced an average of 10 crimes per year organizes a Crime Watch group. During the first year following the creation of the group, 3 crimes are committed in the neighborhood Use the Poisson distribution to calculate the probability the 3 or fewer crimes are committed in a year, assuming that the average number is still 10 times per year Average = 10 x = 3 Answer would be P(X < 3) = 0.00277? Do you think this event provides some evidence that the Crime watch group has been effective in this neighborhood? Justify your answer below. Im lost at this one, any help would be appreciated. Thank you


title: Constructing one-ninth fraction of the $3^{5}$ design

text: I want to construct a $3^{5-2}$ design with $I=ABC$ and $I=CDE$ as generators. The complete defining relation for this design is: $I=ABC=CDE=ABC^{2}DE=ABD^{2}E^{2}$. This is a resolution III design with $x_4=2x_1+2x_2+x_3$ and $x_5=x_1+x_2+2x_4\:({\rm mod}\:3)$. **Question** I wonder how to get $x_4=2x_1+2x_2+x_3$ and $x5=x_1+x_2+2x_4\:({\rm mod}\:3)$ in this specific problem. **Edited** ![enter image description here](http://i.stack.imgur.com/p9x94.png) This is a question from Design and Analysis of Experiments by Douglas C. Montgomery.


title: When to dichotomize a variable for correlational analysis?

text: This is the question I am trying to answer : I’m working on a thesis project and have a variable that is normally distributed and related to happiness. I am considering to split my sample into depressed and non-depressed and run the correlation and the rest of the analyses that way. I think it will help clarify the differences between happy people and non-happy people. a. Would you advise me to do this? Why/Why not? b. What would happen to the correlation if I left it continuous or dichotomized it?


title: recall calculation for a real world example

text: I have totally `696,876` nodes in set _X_. I have `8,127,040` nodes in set _Y_. I am constructing an edge between set _X_ and set _Y_. After constructing the edges between the 2 sets, I have `546,255` edges between set _X_ and set _Y_. Now, my algorithm has detected `3843` edges out of which `3450` are correct edges. From this link, I see that `recall` is the ratio of number of relevant records retrieved to the total number of relevant records in the database.               recall = A/(A+B) *100      Where, > A = number of relevant records retrieved.   >  B = number of relevant records not retrieved. So based on the above formula, I calculate my `recall` value as,               recall = 3450/546255 *100 = 0.631%      Is my `recall` calculation correct?


title: How to select a model in quasi-poisson GLM with interactions using drop1 command?

text: I want to evaluate the effect of three factors (one categorical, and the other two continuous) on the response variable, which is a count data. I have performed 7 candidate GLM models with quasi-Poisson. The models (hypothetically) are the following:               fit1 <- glm(response ~ factor1 * factor2 * factor3, family=quasipoisson, data=mydata)     fit2 <- glm(response ~ factor1 * factor2 + factor3, family=quasipoisson, data=mydata)     fit3 <- glm(response ~ factor1 + factor2 + factor3, family=quasipoisson, data=mydata)     fit4 <- glm(response ~ factor1 * factor2,           family=quasipoisson, data=mydata)     fit5 <- glm(response ~ factor1 + factor2,           family=quasipoisson, data=mydata)     fit6 <- glm(response ~ factor1,                     family=quasipoisson, data=mydata)     fit7 <- glm(response ~ factor2,                     family=quasipoisson, data=mydata)      With these models (some of them with interactions), can I use the `drop1()` command in each model until find the better model? For example:               drop1(fit1, test="F")     drop1(fit2, test="F")     ...     drop1(fit7, test="F")      Because in `fit1` all the factors have interactions, the output shows only one p value. In contraposition, the output of fit3 model (where all the factors do not have interaction) the output shows the three p values. As a result, I do not know if I am proceeding correctly or not for model selection. I do not know if I can use the `drop1` command when I have a model with interactions.


title: Multicollinearity when individual regressions are significant, but VIFs are low

text: I have 6 variables ($x_{1}...x_{6}$) that I am using to predict $y$. When performing my data analysis, I first tried a multiple linear regression. From this, only two variables were significant. However, when I ran a linear regression comparing each variable individually to $y$, all but one were significant ($p$ anywhere from less than 0.01 to less than 0.001). It was suggested that this was due to multicollinearity. My initial research on this suggests checking for multicollinearity by using VIFs. I downloaded the appropriate package from R, and ended up with the resulting VIFs: 3.35, 3.59, 2.64, 2.24, and 5.56. According to various sources online, the point you should be worried about multicollinearity with your VIFs is either at 4 or 5. I am now stumped about what this means for my data. Do I or do I not have a multicollinearity problem? If I do, then how should I proceed? (I cannot collect more data, and the variables are parts of a model that aren't obviously related) If I don't have this problem, then what should I be taking from my data, particularly the fact that these variables are highly significant individually, but not significant at all when combined. **Edit:** Some questions have been asked regarding the dataset, and so I would like to expand... In this particular case, we're looking to understand how specific social cues (gesture, gaze, etc) affect the likelihood of someone producing some other cue. We would like our model to include all significant attributes, so I am uncomfortable removing some that seem redundant. There are not any hypotheses with this right now. Rather, the problem is unstudied, and we are looking to gain a better understanding of what attributes are important. As far as I can tell, these attributes should be relatively independent of one another (you can't just say gaze and gestures are the same, or one the subset of another). It would be nice to be able to report p values for everything, since we would like other researchers to understand what has been looked at. **Edit 2:** Since it came up somewhere below, my $n$ is 24.


title: Analysing log and square-root transformed variables

text: I'm just wondering if it's possible to analyse and interpret variables that have been transformed using different transformations (eg log and square-root) or if all variables in the analysis have to have been transformed with the same method in order to interpret the results. Thanks, Suzie


title: Detecting deviation of one time-series variable from a set of others

text: I have several streams of time-series data representing sensor readings that are usually closely correlated (at least in the casual usage of the term) though with different scalar multipliers applied. For example, say I have time-series A, B, and C representing historical sensor readings over time. B might be consistently around 75% of A, and C might be consistently around 50% of A. The absolute values change frequently, but the ratios between A / B / C remain fairly constant. My goal is to detect points in time when one of these variables deviates significantly from the others. For example, I want to be able to detect if B's readings drop significantly in relation to A's and C's. The number of variables being compared might be anywhere between 3 and the low hundreds. A sketch of what I'm thinking (which may be totally off-base): build a model based on historical data that allows me to predict the values of the target variable based on the values of the other variables. Use the model to predict what the value of the target variable should be for a certain time period, and compare to the actual value for that time period. If the delta between prediction based on historical data and actual reading is large enough, mark the time range as anomalous. While that sounds reasonable to me, I have no idea how to go about: 1 - Building the predictive model from historical data, or 2 - Comparing the predicted and actual values and applying a test of some kind to determine whether the deviation is significant Can anyone recommend approaches for either of these issues?


title: Drawing floating numbers with [0, 1] from uniform distribution by using numpy

text: I'm currently trying to draw floating numbers from a uniform distribution. The Numpy provides numpy.random.uniform.               import numpy as np     sample = np.random.uniform (0, 1, size = (N,) + (2,) + (2,) * K)      However, this module generates values over the half-open interval [0, 1). How can I draw floating numbers with [0, 1] from a uniform distribution? Thanks.


title: Problem with lagged covariates

text: Suppose we have a binary variable $X$ that indicates whether a person ate pizza during the week. This variable is recorded for an entire year. So we have $X_{1}, \dots, X_{52}$ values (1 or 0). Suppose I am interested in the number of weeks of the last 4 weeks that a person did not eat pizza. I want to use this to predict that hazard of a heart attack this week. Call this variable $Y$. Then would $Y = 4-(X_1+X_2+X_3+X_4)$ where $\text{week} > 4$? Note that if we treat $X$ as an array based on time then $X_1 = X[\text{week}-1]$, $X_2 = X[\text{week}-2]$, $X_3 = X[\text{week}-3]$ and $X_4 = X[\text{week}-4]$ since we are working with lagged covariates.


title: Kullback-Leibler vs Hellinger Distance

text: I am working on this problem in which I have a dataset of n-dimensional examples that come from different and unknown distributions. Given a new sample, I wish to find k examples from the dataset that come from distribution(s) closest to the new sample. Which measure(Kullback-Leibler vs Hellinger Distance) might be more suitable for this and why?


title: Computing type II error $\beta$

text: In hypothesis testing for the value of the population mean, when computing type II error $\beta$, if the alternative hypothesis includes many possibilities, is it always the case that $\beta$ is approximated? In other words, when we compute a numerical value for $\beta$, based on the mean from the sample $\bar{X}$, shouldn't we write $\beta \approx P[\text{do not reject } H_0 \mid \mu = \bar{X} ]$? If there is no other information available, of course, it's the only way to compute $\beta$. It's the best we can do. But since $\bar{X}$ may very well not be $\mu$, shouldn't it be stressed that it's $\approx$ and not $=$? Edit: the original non-answered question was here... http://math.stackexchange.com/questions/298922/computing-type-ii-error-beta


title: Classification techniques for incremental attributes

text: I am really interested in learning techniques which can help me deal with the problem briefly described below. I have a two-class {$0,1$} classification problem. Most of the object attributes are incremental in time: when new object appears we know very little about it; as time passes by its "incremental" dynamic attributes are being updated.   Most of these "dynamic" features are incremental - for example such as   * number of visitors of the web page   * number of times the product was purchased etc. Objects arrive in discrete time. Object's label is quite expensive - in the sense that it is rather difficult for expert to classify newly appeared object. False Positive $(\widehat{y}=1, y=0)$ is VERY BAD, while False Negative $(\widehat{y}=0, y=1)$ is not so much. Main goal is to find some classification model, which achieves good balance between   1. model's high True Positive Rate with really low False Positive Rate and   2. time it takes to accumulate object's dynamic attributes for its confident classification by the model. It is crucial to classify new object $(X,Y)$ as soon as possible if its class label is $1(Y=1)$.


title: Sample size effects on R squared

text: Let's say I generated 2*n values $X_{1},X_{2},...X_{n},Y_{1},Y_{2}...Y_{n}$ from a normal distribution N(10,15) and n = 10,100,1000,20000 and fitted the regression line. Here's some result $$ n=10, R^{2} =0.03919\\ n=100,R^{2} = 0.004381\\ n=1000, R^{2} = 0.0001705\\ n=20000,R^{2} = 2.386e-06 $$ from the result as n (sample size) start to increase, $R^{2}$ start to get smaller. how can i explain this phenomena (a hint will do)? thanks


title: How to predict one time-series from another time-series, if they are related

text: I have been trying to solve this problem for over a year without much progress. It is part of a research project I'm doing, but I will illustrate it with a story example I made up, because the actual domain of the problem is a bit confusing (eye-tracking). You are a plane tracking an enemy ship that travels across the ocean, so you have collected a series of (x,y,time) coordinates of the ship. You know that a hidden submarine travels with the ship to protect it, but while there is a correlation between their positions, the submarine often wanders off from the ship, so while it's often near it, it can also be on the other side of the world occasionally. You want to predict the path of the submarine, but unfortunately it is hidden from you. But one month in April you notice the submarine forgets to hide itself, so you have a series of coordinates for both the submarine and the ship throughout 1,000 trips. Using this data, you'd like to build a model to predict the hidden submarine's path given just the ship's movements. The naive baseline would be to say "submarine position guess = "ship's current position" but from the April data where the submarine was visible, you notice there is a tendency for the submarine to be ahead of the ship a bit, so "submarine position guess = ship's position in 1 minute" is an even better estimate. Furthermore, the April data shows that when the ship pauses in the water for an extended period, the submarine is likely to be far away patrolling the coastal waters. There are other patterns of course. How would you build this model, given the April data as training data, to predict the submarine's path? My current solution is an ad-hoc linear regression where the factors are "trip time", "ship's x coordinate", "was ship idle for 1 day", etc. and then having R figure out the weights and doing a cross-validation. But I would really love a way to generate these factors automatically from the April data. Also, a model that uses sequence or time would be nice, since the linear regression doesn't and I think it's relevant. Thanks for reading through all this and I would be happy to clarify anything.


title: How do I model interaction of three boolean variables for logistic regression when there are empty cells?

text: I collected data to find whether the presence or absence of vision, sound, and touch during a task affected the successful completion of that task. However, there were no samples collected where all three senses were absent. So the dependent variable is boolean success but I have a question about how to model the independent variables in a logistic regression. My initial analysis used a single categorical variable with seven levels representing each combination of senses (seven because there were no cases where all three senses were absent).               summary( glmer( Success ~ Condition + ( 1 | Participant ), family=binomial, data=trials))      When I tried to build a model with the Vision, Sound, and Touch as separate variables, the analysis fails. I believe this is because I have empty cells when including the vision*sound*touch interaction because we did not collect results where all senses were absent.               summary( glmer( Success ~ Vision + Sound + Touch + Vision*Sound + Vision*Touch +                      Sound*Touch + Vision*Sound*Touch + ( 1 | Participant ),                      family=binomial, data=trials))      I followed the suggestion linked above to use the `interaction` function to drop the unused factor (all three senses absent). However, this seems to create a variable that looks like my original single categorical variable.               senses <- interaction( trials$Vision, trials$Sound, trials$Touch, drop=TRUE )     summary( glmer( Success ~ senses + ( 1 | Participant ), family=binomial, data=trials))      As I try to refine this analysis, is there a way to model the senses as separate variables to make the interaction between these variables clearer? That is, to appropriately model the contribution of vision in the `vision`, `vision*sound`, `vision*touch` and `vision*sound*touch` conditions. From the initial analysis, the `vision*sound*touch` interaction is the most interesting.


title: Testing for structural breaks in GARCH Models

text: I am looking for a package in R that can test for structural Breaks in GARCH models. I have estimated my coefficients with `rugarch`, and I am highly suspicious that there might be some structural break happening, however I need to test this more formally. Any help would be appreciated :)


title: Testing usefulness of predictors in MLR

text: How can I test whether 2 predictors/covariates in MLR are really useful in prediction using R? I reckon this is something to do with multicollinearity. So would the CAR library function , `VIF()` be appropriate? or can I just do it using simple `lm()` + `summary()`? Thanks. Update: Is this correct?               > model1 = lm( Y ~ X1)     > model2 = lm( Y ~ X1 + X2)     > anova(model1, model2)     Model 1: Y ~ X1     Model 2: Y ~ X1 + X2       Res.Df     RSS Df Sum of Sq      F  Pr(>F)       1     36 0.31814                                   2     35 0.27624  1  0.041904 5.3093 0.02727 *      ok so how to interpret the result. Pr=0.02 < 0.05 so can i say X2 has an effect and is useful as a predictor ? how about the F value? If i look up F(0.95,1,35) in F table then, F=5.31 > 4.17 = should include X2 ? is this correct? thanks


title: Dealing with project uncertainties: is the sum of the most-likely estimates equal to the sum of the expected times?

text: I found this document on the internet: Dealing with Project Uncertainties In this article I read: > In order for uncertainties to be included in the project estimates it is > necessary to take into account an optimistic estimate ( **o** ), a > pessimistic estimate ( **p** ), and a most-likely estimate ( **m** ). This > approach is referred to as stochastic, or statistical estimating, and is the > correct technique to use when developing what is known today as a PERT > (Program Review and Evaluation Technique) So to obtain the estimate for 1 task one can use this formula: $$ Te = \frac{o+4m+p}{6} $$ It then presents us with this formula (I had already heard of this formula (in the book The black art of Software Estimation)): $$ D = \sum T_{e} + Z\sqrt{\sum}\sigma ^{2} $$ Where $\sum T_{e}$ = the sum of all path expected times and $\sum \sigma ^{2}$ = the sum of all critical path variances and $Z$ is the number of standard deviations of a normal distribution (the standard normal deviate). And D is the **project’s duration** at a desired level of confidence "Z", such as 1.281 for 90% confidence To obtain the variance it presents the formula: $$ Variance = \sigma^{2} = (\frac{p-o}{6})^{2} $$ The article also says: > Regardless of the shape formed by the data values of the estimates, the > shape of the values making up the sampling distribution of the expected > times will approach a normal distribution based on the central limit > theorem, if the sampling size is sufficiently large (usually 30 or greater). > When the sampling size is less than 30, and the Central Limit Theorem cannot > be invoked, a **tdistribution** must be used Now the questions: ¿Is this $D = \sum T_{e} + Z\sqrt{\sum}\sigma ^{2}$ an accurate way to sum estimates? I asked before about file storage estimation, and this formula was not used to provide me with answer... ¿why is that? ¿is it because did not have different estimations with different optimistic, pessimistic and most likely sizes? ¿If sampling size is less than 30... how can I use this "tdistribution" to sum estimates ? ¿what is the formula then?